遗传力又称遗传率，指遗传方差在总方差（表型方差）中所占的比值，可以作为杂种后代进行选择的一个指标。遗传力分为单株遗传力、家系遗传力、小区遗传力、个体遗传力。动物中一般用个体遗传力，植物中一般用家系遗传力。
遗传力介绍详细介绍见邓飞老师博客https://zhuanlan.zhihu.com/p/368057210?ivk_sa=1024320u
对于不同的数据，遗传力计算方法有所不同，本篇文章是对多年单点有重复数据进行遗传力计算。
读取数据
#设置工作目录
> setwd("D:/GWAS_phe")
#调用R包
> library('Matrix')
> library('lme4')
#读取表型数据（这里需要原始数据）
> 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个品种），第二列为区组Blk（三个区组、每个区组5个单株重复）、第三列为年份Year（两年），第4列为性状。
在计算前，需要将考虑的因素变为因子(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 ...
计算
参照某平台的课程，到这里，一切都是熟悉的样子，天真的我开始和大多数时候一样的操作，复制粘贴——改数据名称，内心毫无波澜，甚至有些急迫地等待结果好继续下面的分析，然而……
> 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)
这时出现了一长串warning，警告lme4:模型无法与max | grad |收敛，但是有输出结果，所以只是一晃而过，也没有太在意，然后查看一下：
> 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)
哎？？？品种、年份、残差、品种和年份的互作，方差竟然一模一样！这，这就不对了吧！为了找到问题所在，不再忽略警告，重新运行一次。这里建议大家，不要一开始就选择忽略所有warning，这里的数据是恰好算出来一模一样，如果不是这样的话，很容易被误导，拿到错误结果。
> m1 =lmer(TL~(1|Blk%in%Year)+(1|Year)+(1|Cul)+(1|Year:Cul),data=TL)
错误: grouping factors must have > 1 sampled level
结果依然报错：grouping factors must have > 1 sampled level，并且没有输出结果。
继续搜帖子，然后发现，教程里要不然就是单年多点的数据，要不然就是多年多点有重复的数据，为什么没有多年单点有重复的呢？多年多点的模型如下：
m =lmer(Trait~(1|Line)+(1|Year)+(1|Loc)+(1|Line:Loc) +(1|Line:Year),data=data)
emm……为什么这么模型里完全没出现区组和重复呢？
参者这个模型，如果把Loc全部删掉，那我的重复就没有任何意义了；如果把这个模型里的Loc都替换成Blk，也不对，Blk不是一个独立的因子，不能单独存在于函数里。最符合我的理解的，Blk%in%Year，重复嵌套在年份里，也不对，想到前面Blue值计算中，Blk和Year的互作效应，继续尝试：
> 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虽然出现了warning，但是有运算结果。可是我也不能确定，这个是不是准确的结果。
SAS遗传力计算
为验证R中结果的可靠性，利用SAS进行了计算验证：
proc mixed data=dat;
class Year Blk Cul;
model TL = / solution;
random Year Blk(Year) Cul Cul*Year / solution;
run;
从上到下依次为环境方差、区组方差、基因型方差、基因型与环境互作方差、误差方差。
数字上来看，SAS与m2的结果基本一致。
我的疑问在于，SAS中，写法为Blk(Year)和Cul*Year，分别是嵌套和互作，但是为什么在lme4中，都是(1|Year:Cul) 和(1|Blk:Year)交互的写法？而且这样得到的结果竟然是一致的。如果有大佬理解其中的原理，还烦请浪费几分钟，告诉我为什么，不胜感激！
到这里为止，各组分的方差终于可以确定，剩下的部分就是套公式了，公式如下：
Vg：遗传方差（Cul,131.91)
Vge：基因与环境的互作方差 (Year：Cul, 126.75)
l：环境个数 (年份：2）
VΣ：残差 （Residual: 143.54)
r：区组个数 （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
当然，在上述模型中没有考虑另外一些互作，比如Cul:Blk，Cul:Blk:Year等等，是因为互作考虑的太多，遗传力计算会很复杂，所以这样设置模型主要是便于计算。
引用转载请注明出处，如有错误敬请指出。
网友评论