如图

邓飞

发布于 2020-11-26 11:07:36

5460

发布于 2020-11-26 11:07:36

1. 背景

很多多年多点，或者一年多点的数据，要计算遗传力，关于数据如何整理，关于模型如何定义，遗传力公式如何写，如何考虑地点、年份、重复的个数，这些问题难倒了很多想要作此分析的小伙伴。这里，根据一个我之前上传到B站的视频，把里面的数据和代码进行演示如何计算相关的参数。另外，视频中也有一些错误或者不足的地方，我做了说明，后面我用红色字体标识了一下。

另外，我本次推送时，也把我之前录制的几期视频也传到了公众号上，主要有：

下载安装最新版的R和RStudio

安装git和plink软件

RStudio使用的10个技巧

如何安装GWAS软件包：GAPIT

如果我有新电脑一定要安装这四款软件

感兴趣的同学可以关注本公众号好，查看推送内容。

2. 本次微信文的目标

获得一个多年多点的数据
计算品种性状的遗传力
计算每个品种的育种值(BLUP)

3. 试验设计

两年: 2009和2010年
两点: CA和OH
每个年, 每个地点, 2个重复(其中OH在2010年仅有1个重复)
试验设计是完全随机区组(RCBD)
共有143个品种
按照小区(plot)进行考种

4. 数据探索性分析

「预览数据」数据包括品种(Line), 重复(Rep), 年份(Year), 地点(Loc), 收获日期(Harvest), 产量(Yield), Brix, PH, TA这三个也是观测值, 主要是品种性状.

需要分析的观测值: Yield, Brix, pH, TA 需要考虑的因子: Line, Rep, Year, Loc

> head(dat)
      Line Rep Year Loc Harvest    Yield Brix   pH   TA
1 SCT_0001   1 2009  OH  9/7/13       NA 5.70 4.42 0.39
2 SCT_0001   1 2009  CA 9/23/10 38789.80 4.75 4.47 0.34
3 SCT_0001   1 2010  CA    <NA>       NA 4.34 4.44 0.29
4 SCT_0001   2 2009  OH 9/14/13       NA 5.40 4.48 0.38
5 SCT_0001   2 2009  CA 9/14/10       NA 4.84 4.67 0.25
6 SCT_0001   2 2010  OH 9/30/14 93832.21 5.40 4.58 0.31

「查看Brix的直方图」

hist(dat$Brix,col="gold")

查看Brix在不同地点的箱线图:

boxplot(dat$Brix~dat$Loc, xlab="Location", ylab="Degrees Brix", main="Degrees Brix by Location", col="pink")

5. 重新转化数据

这里建模之前, 需要对数据进行转化, 将需要考虑的因素变为因子(Factor), 将需要分析的性状变为数值(number)

> str(dat)
'data.frame': 986 obs. of  9 variables:
 $ Line   : Factor w/ 143 levels "SCT_0001","SCT_0002",..: 1 1 1 1 1 1 1 2 2 2 ...
 $ Rep    : int  1 1 1 2 2 2 2 1 1 1 ...
 $ Year   : int  2009 2009 2010 2009 2009 2010 2010 2009 2009 2010 ...
 $ Loc    : Factor w/ 2 levels "CA","OH": 2 1 1 2 1 2 1 2 1 1 ...
 $ Harvest: Factor w/ 26 levels "10/10/10","10/2/10",..: 25 19 NA 11 10 24 NA 25 9 NA ...
 $ Yield  : num  NA 38790 NA NA NA ...
 $ Brix   : num  5.7 4.75 4.34 5.4 4.84 5.4 5.16 5.5 5.65 4.54 ...
 $ pH     : num  4.42 4.47 4.44 4.48 4.67 4.58 4.47 4.44 4.47 4.52 ...
 $ TA     : num  0.39 0.34 0.29 0.38 0.25 0.31 0.29 0.44 0.33 0.23 ...
> for(i in 1:5) dat[,i] = as.factor(dat[,i])
> str(dat)
'data.frame': 986 obs. of  9 variables:
 $ Line   : Factor w/ 143 levels "SCT_0001","SCT_0002",..: 1 1 1 1 1 1 1 2 2 2 ...
 $ Rep    : Factor w/ 2 levels "1","2": 1 1 1 2 2 2 2 1 1 1 ...
 $ Year   : Factor w/ 2 levels "2009","2010": 1 1 2 1 1 2 2 1 1 2 ...
 $ Loc    : Factor w/ 2 levels "CA","OH": 2 1 1 2 1 2 1 2 1 1 ...
 $ Harvest: Factor w/ 26 levels "10/10/10","10/2/10",..: 25 19 NA 11 10 24 NA 25 9 NA ...
 $ Yield  : num  NA 38790 NA NA NA ...
 $ Brix   : num  5.7 4.75 4.34 5.4 4.84 5.4 5.16 5.5 5.65 4.54 ...
 $ pH     : num  4.42 4.47 4.44 4.48 4.67 4.58 4.47 4.44 4.47 4.52 ...
 $ TA     : num  0.39 0.34 0.29 0.38 0.25 0.31 0.29 0.44 0.33 0.23 ...

可以看到, 转化之后, 前五列变为了Factor, 后四列为num

6. 建模

这里使用混合线性模型, 使用R包lme4

观测值: Brix 固定因子: 无；随机因子: Line + Loc + Year + Loc.Year.Rep + Line:Loc + Line:Year + Line:Year:Loc

library(lme4)
library(lmerTest)
mod1 = lmer(Brix ~ (1|Line) + (1|Loc) + (1|Year) + 
              (1|Line:Loc) + (1|Line:Year) , data=dat)
summary(mod1)

结果:

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: Brix ~ (1 | Line) + (1 | Loc) + (1 | Year) + (1 | Line:Loc) +      (1 | Line:Year)
   Data: dat

REML criterion at convergence: 1742.1

Scaled residuals: 
    Min      1Q  Median      3Q     Max 
-3.0258 -0.6133 -0.0533  0.5293  5.0653 

Random effects:
 Groups    Name        Variance Std.Dev.
 Line:Year (Intercept) 0.010209 0.10104 
 Line:Loc  (Intercept) 0.009483 0.09738 
 Line      (Intercept) 0.191933 0.43810 
 Year      (Intercept) 0.031648 0.17790 
 Loc       (Intercept) 0.188486 0.43415 
 Residual              0.257820 0.50776 
Number of obs: 967, groups:  Line:Year, 284; Line:Loc, 281; Line, 143; Year, 2; Loc, 2

Fixed effects:
            Estimate Std. Error     df t value Pr(>|t|)  
(Intercept)   5.2817     0.3343 1.3583    15.8   0.0167 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

7. 遗传力的计算

这里 Vg 是Line的方差组分: 0.191933；VGL 是品种与地点的交互方差组分: 0.009483；VGY 为品种与年份的交互方差组分: 0.010209；Ve 为残差方差组分: 0.257820

这里, 品种与地点以及年份互作, 没有考虑, 因此不做计算.

地点为2, 年份为2, 重复为2

h^2 = \frac{ 0.191933}{ 0.191933 + 0.009483/2 + 0.010209/2 +0.257820/8 }

> vg = 0.191933
> vgl = 0.009483
> ve = 0.257820
> # 遗传力
> vg = 0.191933
> vgl = 0.009483
> vgy = 0.010209
> ve = 0.257820
> h2 = vg/(vg + vgl/2 + vgy/2 + ve/8)
> h2
[1] 0.8202037

8. 计算BLUP值

> brixblup = ranef(mod1)$Line
> head(brixblup)
         (Intercept)
SCT_0001 -0.13628589
SCT_0002 -0.14944411
SCT_0003 -0.05347823
SCT_0004 -0.03878384
SCT_0005  1.13568896
SCT_0006  0.49028454
> write.csv(brixblup,"brixblup.csv")

9. 对比BLUP值和平均值

可以看出, BLUP值和平均值趋势基本一致, 但是有个别品种, BLUP值和平均值变化较大.

mm = as.data.frame(tapply(dat$Brix, dat$Line, na.rm=T, mean))
names(mm) = "mean"
plot(brixblup$blup,mm$mean)

10. 不足

这篇无疑是开山之作, 但是也有一些不足:

一般来说, 多年多点分析中, 我们将地点, 年份, 地点:年份, 地点:年份:重复作为固定因子, 品种, 品种与地点, 品种与年份, 品种与地点与年份作为随机因子. 文章中将所有因子都作为随机因子, 太过武断.
遗传力公式有错误:

这里没有考虑: 品种:年份:地点, 残差的分母应该是2*2*2 = 8 , 而不是4.

因素没有考虑完整, 可能是数据量有限, 没有考虑地点:年份:重复, 没有考虑地点:年份:品种
计算遗传力没有标准误, 标准误可以反映出计算的好坏.
没有给出随机因子的显著性, 可以使用LRT检验

11. 商业软件asreml的解决方案

# asrmel的解决方案
library(asreml)
library(learnasreml)
head(dat)
mod2 = asreml(Brix ~ 1 , random = ~Loc + Year + Line+Line:Year + Line:Loc, data=dat)
summary(mod2)$varcomp
pin(mod2, h2 ~ V3/(V3+V4/2+V5/2 + V6/8))

> mod2 = asreml(Brix ~ 1 , random = ~Loc + Year + Line+Line:Year + Line:Loc, data=dat)
ASReml: Tue May 14 21:03:02 2019

     LogLik         S2      DF      wall     cpu
    -27.8710      0.3114   966  21:03:02     0.0
     -8.1437      0.2968   966  21:03:02     0.0
      8.3477      0.2802   966  21:03:02     0.0
     16.0150      0.2642   966  21:03:02     0.0
     16.6323      0.2585   966  21:03:02     0.0
     16.6563      0.2579   966  21:03:02     0.0
     16.6571      0.2578   966  21:03:02     0.0
     16.6571      0.2578   966  21:03:02     0.0

Finished on: Tue May 14 21:03:02 2019
 
LogLikelihood Converged 
> summary(mod2)$varcomp
                        gamma   component  std.error    z.ratio constraint
Loc!Loc.var        0.73097982 0.188461143 0.26741513  0.7047512   Positive
Year!Year.var      0.12272634 0.031641291 0.04564871  0.6931475   Positive
Line!Line.var      0.74444586 0.191932955 0.02948641  6.5092010   Positive
Line:Year!Line.var 0.03960966 0.010212159 0.01139269  0.8963779   Positive
Line:Loc!Line.var  0.03677090 0.009480269 0.01125209  0.8425343   Positive
R!variance         1.00000000 0.257819899 0.01557455 16.5539268   Positive
> pin(mod2, h2 ~ V3/(V3+V4/2+V5/2 + V6/8))
   Estimate        SE
h2 0.820203 0.0394735

「更简洁, 更友好, 更方便.」

公众号回复: blup 下载文档(PPT), 数据, 代码:

「练习」

计算其它三个性状的遗传力
比较性状BLUP值和平均值的异同, 考虑为何要用BLUP作为选择标准
考虑还有没有其它分析的切入点

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

原始发表：2020-11-24，如有侵权请联系 cloudcommunity@tencent.com 删除

python

本文分享自育种数据分析之放飞自我微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

python

登录后参与评论

0 条评论

热度

如图

如图