根据固定因子和协变量矫正表型值做GWAS和GS

邓飞

发布于 2026-06-03 13:01:44

540

大家好，我是邓飞。

我们在做GWAS和GS的时候，有时候会考虑固定因子的影响，比如植物中的年份、地点，比如动物中的场年季性别效应，对于这些固定因子，常用的做法是放到一起建模：

GWAS：y ~ Fixed + covar +SNP

GS: y ~ Fixed + covar + ADD

固定因子和协变量，固定因子是Factor，协变量是数字变量，比如年龄，PCA结果等，有些软件区分，比如GCTA的dcovar和qcovar，有些软件只支持数字协变量，比如plink和GEMMA（根红苗正的GWAS软件：GEMMA），这时候就需要用plink的--dummy将固定因子变为哑变量（数字协变量）才能分析（plink软件cookbook）。

为了对数据进行统一操作，我一般是将所有的固定因子都变为数字协变量，这样所有的软件都支持操作了。GWAS分析没什么问题。

问题是GS分析，有时候会用到交叉验证，有固定因子在时，表型数据里面有固定效应，直接用于5层抽样结果不准确，正确的做法是计算固定因子的矫正值（残差值+平均值）（用固定因素矫正表型值，要用残差值？）

所以，考虑固定因子和协变量的GWAS和GS分析过时了，先矫正，最后模型中不考虑固定因子是省事的做法。

理论上是这样讲，但是需要拿一套数据比较一下，心里才有谱：

> mod1 = asreml(BWT ~ BYEAR + SEX, random = ~ vm(ANIMAL,ainv), residual = ~ idv(units),data=d1)
ASReml Version 4.2 02/06/2026 17:27:30
          LogLik        Sigma2     DF     wall
 1     -2707.455           1.0    819   17:27:30
 2     -2207.352           1.0    819   17:27:30
 3     -1665.874           1.0    819   17:27:30
 4     -1307.258           1.0    819   17:27:30
 5     -1133.910           1.0    819   17:27:30
 6     -1096.411           1.0    819   17:27:30
 7     -1093.234           1.0    819   17:27:31
 8     -1093.197           1.0    819   17:27:31
 9     -1093.197           1.0    819   17:27:31
> summary(mod1)$varcomp
                 component std.error  z.ratio bound %ch
vm(ANIMAL, ainv)  2.685783 0.4437252 6.052807     P   0
units!units       2.324005 0.3475859 6.686132     P   0
units!R           1.000000        NA       NA     F   0

> vpredict(mod1,h2 ~ V1/(V1+V2))
    Estimate         SE
h2 0.5361071 0.0743135

这个是两个固定因子的动物模型，加性方差组分为2.68，残差方差组分为2.32，遗传力为0.53。

我们下面先对BWT进行固定因子的矫正，就是表型数据中减去这两个固定因子的效应值，简单来说就是：effect + mean

> temp = asreml(BWT ~ BYEAR + SEX, residual = ~ idv(units),data=d1)
ASReml Version 4.2 02/06/2026 17:29:29
          LogLik        Sigma2     DF     wall
 1     -5362.537           1.0    819   17:29:29
 2     -4169.383           1.0    819   17:29:29
 3     -2843.121           1.0    819   17:29:29
 4     -1897.806           1.0    819   17:29:29
 5     -1356.519           1.0    819   17:29:29
 6     -1172.507           1.0    819   17:29:29
 7     -1131.533           1.0    819   17:29:29
 8     -1127.845           1.0    819   17:29:29
 9     -1127.796           1.0    819   17:29:29
10     -1127.796           1.0    819   17:29:29
> d1$BWT_correct = resid(temp) + mean(d1$BWT,na.rm=T)
> head(d1)
  ANIMAL MOTHER BYEAR SEX   BWT TARSUS BWT_correct
1   1029   1145   968   1 10.77  24.77    8.100141
2   1299    811   968   1  9.30  22.46    6.630141
3    643    642   970   2  3.98  12.89    3.541682
4   1183   1186   970   1  5.39  20.47    6.872058
5   1238   1237   970   2 12.12     NA   11.681682
6    891    895   970   1    NA     NA          NA

可以看到，新的BWT_correct就是根据BYEAR和SEX校正后的BWT值，用它进行动物模型建模，比较一下两者的结果：

> mod2 = asreml(BWT_correct ~ 1, random = ~ vm(ANIMAL,ainv), residual = ~ idv(units),data=d1)
ASReml Version 4.2 02/06/2026 17:30:26
          LogLik        Sigma2     DF     wall
 1     -3880.040           1.0    853   17:30:26
 2     -3054.503           1.0    853   17:30:26
 3     -2146.341           1.0    853   17:30:26
 4     -1517.504           1.0    853   17:30:26
 5     -1180.624           1.0    853   17:30:26
 6     -1084.535           1.0    853   17:30:26
 7     -1069.854           1.0    853   17:30:26
 8     -1069.258           1.0    853   17:30:26
 9     -1069.257           1.0    853   17:30:26
> summary(mod2)$varcomp
                 component std.error  z.ratio bound %ch
vm(ANIMAL, ainv)  2.501618 0.4153511 6.022898     P   0
units!units       2.315922 0.3275531 7.070373     P   0
units!R           1.000000        NA       NA     F   0
> blup1 = coef(mod1)$random %>% as.data.frame()
> blup2 = coef(mod2)$random %>% as.data.frame()
> cor(blup1$effect,blup2$effect)
[1] 0.9977903

方差组分有点变化，blup相关系数为0.997，基本一致。

为何相关系数不是1？为何方差组分不完全一样？

我问了一下豆包结合自己的理解：reml方法是迭代的，另外，固定效应本身有误差，所以会有一些出入，但是相关系数0.99，可以认为两者一致了。

比较一下原始的BWT和BWT_correct值：

> cor(d1$BWT,d1$BWT_correct,use="complete.obs")
[1] 0.8191949

可以看到，矫正后的值相关系数0.819，这个变化还是比较大的。

看一下校正后的值不同固定因子的分布：

> d1 %>% group_by(SEX) %>% summarise(xx = mean(BWT,na.rm=T), yy = mean(BWT_correct,na.rm=T))
# A tibble: 2 × 3
  SEX      xx    yy
  <fct> <dbl> <dbl>
1 1      6.10  7.37
2 2      8.33  7.37
> d1 %>% group_by(BYEAR) %>% summarise(xx = mean(BWT,na.rm=T), yy = mean(BWT_correct,na.rm=T))
# A tibble: 34 × 3
   BYEAR    xx    yy
   <fct> <dbl> <dbl>
 1 1000   4.78  7.37
 2 1001   5.08  7.37
 3 1002   6.94  7.37
 4 968   10.0   7.37
 5 970    7.16  7.37
 6 971    6.75  7.37
 7 972    7.37  7.37
 8 973    7.58  7.37
 9 974    6.84  7.37
10 975    6.95  7.37

可以看到，校正后的值对不同固定因子不同水平结果是一样的，也就是说固定因子的效应已经去掉了，用校正后的表型值去分析不用考虑固定因子了。

这里有一个坑：动物不同场年季，ID不会有重复，植物中如果不同年份或者不同地点重复品种比较多，这种方法不是最优的，最优的是计算品种的BLUE值，作为矫正值进行后续的GWAS和GS分析（GWAS多环境表型数据用BLUE值还是BLUP值？）

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2026-06-02，如有侵权请联系 cloudcommunity@tencent.com 删除