GWAS | 遗传力计算

遗传力又称遗传率，指遗传方差在总方差（表型方差）中所占的比值，可以作为杂种后代进行选择的一个指标。遗传力分为单株遗传力、家系遗传力、小区遗传力、个体遗传力。动物中一般用个体遗传力，植物中一般用家系遗传力。

遗传力介绍详细介绍见邓飞老师博客https://zhuanlan.zhihu.com/p/368057210?ivk_sa=1024320u

https://cloud.tencent.com/developer/article/1445670

对于不同的数据，遗传力计算方法有所不同，本篇文章是对多年单点有重复数据进行遗传力计算。

读取数据

#设置工作目录
> 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;

image

从上到下依次为环境方差、区组方差、基因型方差、基因型与环境互作方差、误差方差。
数字上来看，SAS与m2的结果基本一致。
我的疑问在于，SAS中，写法为Blk(Year)和Cul*Year，分别是嵌套和互作，但是为什么在lme4中，都是(1|Year:Cul) 和(1|Blk:Year)交互的写法？而且这样得到的结果竟然是一致的。如果有大佬理解其中的原理，还烦请浪费几分钟，告诉我为什么，不胜感激！

到这里为止，各组分的方差终于可以确定，剩下的部分就是套公式了，公式如下：

image

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等等，是因为互作考虑的太多，遗传力计算会很复杂，所以这样设置模型主要是便于计算。

转自：https://www.jianshu.com/p/d116831897a2