遺傳力又稱遺傳率，指遺傳方差在總方差（表型方差）中所占的比值吊骤，可以作為雜種后代進(jìn)行選擇的一個(gè)指標(biāo)缎岗。遺傳力分為單株遺傳力、家系遺傳力白粉、小區(qū)遺傳力传泊、個(gè)體遺傳力。動(dòng)物中一般用個(gè)體遺傳力蜗元，植物中一般用家系遺傳力或渤。
遺傳力介紹詳細(xì)介紹見鄧飛老師博客https://zhuanlan.zhihu.com/p/368057210?ivk_sa=1024320u
https://cloud.tencent.com/developer/article/1445670
對于不同的數(shù)據(jù)，遺傳力計(jì)算方法有所不同奕扣，本篇文章是對多年單點(diǎn)有重復(fù)數(shù)據(jù)進(jìn)行遺傳力計(jì)算薪鹦。
讀取數(shù)據(jù)
#設(shè)置工作目錄
> setwd("D:/GWAS_phe")
#調(diào)用R包
> library('Matrix')
> library('lme4')
#讀取表型數(shù)據(jù)（這里需要原始數(shù)據(jù)）
> dat <- read.table("TL.txt", header = T, check.names = F, sep = "\t")
> head(dat)
Cul Blk Year TL
1 1 1 2017 40.37
2 1 1 2017 62.99
3 1 1 2017 90.68
4 1 1 2017 42.09
5 1 1 2017 57.25
6 1 2 2017 25.30
#第一列為品種Cul（188個(gè)品種），第二列為區(qū)組Blk（三個(gè)區(qū)組惯豆、每個(gè)區(qū)組5個(gè)單株重復(fù)）池磁、第三列為年份Year（兩年），第4列為性狀楷兽。
在計(jì)算前地熄，需要將考慮的因素變?yōu)橐蜃?Factor)
> for(i in 1:3) dat[,i] = as.factor(dat[,i]) #前三列
> str(dat)
'data.frame': 5640 obs. of 4 variables:
$ Cul : Factor w/ 188 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
$ Blk : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 2 2 2 2 2 ...
$ Year: Factor w/ 2 levels "2017","2020": 1 1 1 1 1 1 1 1 1 1 ...
$ TL : num 40.4 63 90.7 42.1 57.2 ...
計(jì)算
參照某平臺的課程，到這里芯杀，一切都是熟悉的樣子端考，天真的我開始和大多數(shù)時(shí)候一樣的操作，復(fù)制粘貼——改數(shù)據(jù)名稱揭厚，內(nèi)心毫無波瀾却特，甚至有些急迫地等待結(jié)果好繼續(xù)下面的分析，然而……
> options(lmerControl=list(check.nobs.vs.rankZ = "warning",
+ check.nobs.vs.nlev = "warning",
+ check.nobs.vs.nRE = "warning",
+ check.nlev.gtreq.5 = "warning",
+ check.nlev.gtr.1 = "warning"))
> m1 =lmer(TL~(1|Blk%in%Year)+(1|Year)+(1|Cul)+(1|Year:Cul),data=TL)
Warning messages:
1: grouping factors must have > 1 sampled level
2: grouping factors with < 5 sampled levels may give unreliable estimates
3: In checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
Model failed to converge with max|grad| = 1.55852 (tol = 0.002, component 1)
這時(shí)出現(xiàn)了一長串warning筛圆，警告lme4:模型無法與max | grad |收斂裂明，但是有輸出結(jié)果，所以只是一晃而過太援，也沒有太在意闽晦，然后查看一下：
> summary(m1)
Linear mixed model fit by REML ['lmerMod']
Formula: TL ~ (1 | Blk %in% Year) + (1 | Year) + (1 | Cul) + (1 | Year:Cul)
Data: TL
REML criterion at convergence: 38021.9
Scaled residuals:
Min 1Q Median 3Q Max
-4.3273 -0.5796 -0.0600 0.5496 5.6062
Random effects:
Groups Name Variance Std.Dev.
Year:Cul (Intercept) 143.441 11.977
Cul (Intercept) 143.441 11.977
Year (Intercept) 143.441 11.977
Blk %in% Year (Intercept) 1.434 1.198
Residual 143.441 11.977
Number of obs: 4733, groups: Year:Cul, 333; Cul, 188; Year, 2; Blk %in% Year, 1
Fixed effects:
Estimate Std. Error t value
(Intercept) 37.658 8.626 4.365
optimizer (nloptwrap) convergence code: 0 (OK)
Model failed to converge with max|grad| = 1.55852 (tol = 0.002, component 1)
哎？提岔？仙蛉？品種、年份唧垦、殘差捅儒、品種和年份的互作，方差竟然一模一樣！這巧还，這就不對了吧鞭莽！為了找到問題所在，不再忽略警告麸祷，重新運(yùn)行一次澎怒。這里建議大家，不要一開始就選擇忽略所有warning阶牍，這里的數(shù)據(jù)是恰好算出來一模一樣喷面，如果不是這樣的話，很容易被誤導(dǎo)走孽，拿到錯(cuò)誤結(jié)果惧辈。
> m1 =lmer(TL~(1|Blk%in%Year)+(1|Year)+(1|Cul)+(1|Year:Cul),data=TL)
錯(cuò)誤: grouping factors must have > 1 sampled level
結(jié)果依然報(bào)錯(cuò)：grouping factors must have > 1 sampled level，并且沒有輸出結(jié)果磕瓷。
繼續(xù)搜帖子盒齿，然后發(fā)現(xiàn)，教程里要不然就是單年多點(diǎn)的數(shù)據(jù)困食，要不然就是多年多點(diǎn)有重復(fù)的數(shù)據(jù)边翁，為什么沒有多年單點(diǎn)有重復(fù)的呢？多年多點(diǎn)的模型如下：
m =lmer(Trait~(1|Line)+(1|Year)+(1|Loc)+(1|Line:Loc) +(1|Line:Year),data=data)
emm……為什么這么模型里完全沒出現(xiàn)區(qū)組和重復(fù)呢硕盹？
參者這個(gè)模型符匾，如果把Loc全部刪掉，那我的重復(fù)就沒有任何意義了瘩例；如果把這個(gè)模型里的Loc都替換成Blk啊胶，也不對，Blk不是一個(gè)獨(dú)立的因子垛贤，不能單獨(dú)存在于函數(shù)里创淡。最符合我的理解的，Blk%in%Year南吮，重復(fù)嵌套在年份里，也不對誊酌，想到前面Blue值計(jì)算中部凑，Blk和Year的互作效應(yīng)，繼續(xù)嘗試：
> m2 =lmer(TL~(1|Year)+(1|Cul)+(1|Year:Cul) +(1|Blk:Year),data=TL)
Warning message:
grouping factors with < 5 sampled levels may give unreliable estimates
> summary(m2)
Linear mixed model fit by REML ['lmerMod']
Formula: TL ~ (1 | Year) + (1 | Cul) + (1 | Year:Cul) + (1 | Blk:Year)
Data: TL
REML criterion at convergence: 37999.2
Scaled residuals:
Min 1Q Median 3Q Max
-4.2590 -0.5803 -0.0636 0.5481 5.4665
Random effects:
Groups Name Variance Std.Dev.
Year:Cul (Intercept) 126.755 11.259
Cul (Intercept) 131.938 11.486
Blk:Year (Intercept) 1.452 1.205
Year (Intercept) 123.803 11.127
Residual 143.540 11.981
Number of obs: 4733, groups: Year:Cul, 333; Cul, 188; Blk:Year, 6; Year, 2
Fixed effects:
Estimate Std. Error t value
(Intercept) 37.668 7.956 4.735
m2雖然出現(xiàn)了warning碧浊，但是有運(yùn)算結(jié)果涂邀。可是我也不能確定箱锐，這個(gè)是不是準(zhǔn)確的結(jié)果比勉。
SAS遺傳力計(jì)算
為驗(yàn)證R中結(jié)果的可靠性，利用SAS進(jìn)行了計(jì)算驗(yàn)證：
proc mixed data=dat;
class Year Blk Cul;
model TL = / solution;
random Year Blk(Year) Cul Cul*Year / solution;
run;
從上到下依次為環(huán)境方差、區(qū)組方差浩聋、基因型方差观蜗、基因型與環(huán)境互作方差、誤差方差衣洁。
數(shù)字上來看墓捻，SAS與m2的結(jié)果基本一致。
我的疑問在于坊夫，SAS中砖第，寫法為Blk(Year)和Cul*Year，分別是嵌套和互作环凿，但是為什么在lme4中梧兼，都是(1|Year:Cul) 和(1|Blk:Year)交互的寫法？而且這樣得到的結(jié)果竟然是一致的智听。如果有大佬理解其中的原理羽杰，還煩請浪費(fèi)幾分鐘，告訴我為什么瞭稼，不勝感激忽洛！
到這里為止，各組分的方差終于可以確定环肘，剩下的部分就是套公式了欲虚，公式如下：
Vg：遺傳方差（Cul,131.91)
Vge：基因與環(huán)境的互作方差 (Year：Cul, 126.75)
l：環(huán)境個(gè)數(shù) (年份：2）
VΣ：殘差 （Residual: 143.54)
r：區(qū)組個(gè)數(shù) （3）
> result <- summary(m2)
> var <- as.data.frame(result$varcor)
> var
grp var1 var2 vcov sdcor
1 Year:Cul (Intercept) <NA> 126.754983 11.258552
2 Cul (Intercept) <NA> 131.938053 11.486429
3 Blk:Year (Intercept) <NA> 1.451874 1.204937
4 Year (Intercept) <NA> 123.803121 11.126685
5 Residual <NA> <NA> 143.540460 11.980837
> H2 <- var[2,4]/(var[2,4]+var[1,4]/2+var[5,4]/(2*3))
> H2
[1] 0.6018002
當(dāng)然，在上述模型中沒有考慮另外一些互作悔雹，比如Cul:Blk复哆，Cul:Blk:Year等等，是因?yàn)榛プ骺紤]的太多腌零，遺傳力計(jì)算會很復(fù)雜梯找，所以這樣設(shè)置模型主要是便于計(jì)算。
引用轉(zhuǎn)載請注明出處，如有錯(cuò)誤敬請指出。
