【原】根據(jù)多年多點(diǎn)數(shù)據(jù)使用R語言計(jì)算遺傳力和育種值

育種數(shù)據(jù)分析 2021-11-18

展開全文

1. 背景

之前寫過一篇博客, 介紹領(lǐng)導(dǎo)安利我嗶哩嗶哩的故事, 介紹到我將我從YouTube上收集的關(guān)于混合線性模型, 關(guān)于GWAS, 關(guān)于GS, 關(guān)于農(nóng)業(yè)數(shù)據(jù)分析相關(guān)的視頻, 上傳到了嗶哩嗶哩上面. 今天我們看一下介紹多年多點(diǎn)遺傳力及BLUP值計(jì)算的視頻內(nèi)容. 閱讀原文可以查看視頻, 這里我用文字和代碼進(jìn)行重演.

2. 本次微信文的目標(biāo)

獲得一個(gè)多年多點(diǎn)的數(shù)據(jù)
計(jì)算品種性狀的遺傳力
計(jì)算每個(gè)品種的育種值(BLUP)

3. 試驗(yàn)設(shè)計(jì)

兩年: 2009和2010年
兩點(diǎn): CA和OH
每個(gè)年, 每個(gè)地點(diǎn), 2個(gè)重復(fù)(其中OH在2010年僅有1個(gè)重復(fù))
試驗(yàn)設(shè)計(jì)是完全隨機(jī)區(qū)組(RCBD)
共有143個(gè)品種
按照小區(qū)(plot)進(jìn)行考種

4. 數(shù)據(jù)探索性分析

「預(yù)覽數(shù)據(jù)」數(shù)據(jù)包括品種(Line), 重復(fù)(Rep), 年份(Year), 地點(diǎn)(Loc), 收獲日期(Harvest), 產(chǎn)量(Yield), Brix, PH, TA這三個(gè)也是觀測(cè)值, 主要是品種性狀.

需要分析的觀測(cè)值: 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在不同地點(diǎn)的箱線圖:

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

5. 重新轉(zhuǎn)化數(shù)據(jù)

這里建模之前, 需要對(duì)數(shù)據(jù)進(jìn)行轉(zhuǎn)化, 將需要考慮的因素變?yōu)橐蜃?Factor), 將需要分析的性狀變?yōu)閿?shù)值(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 ...

可以看到, 轉(zhuǎn)化之后, 前五列變?yōu)榱薋actor, 后四列為num

6. 建模

這里使用混合線性模型, 使用R包lme4

觀測(cè)值: Brix 固定因子: 無隨機(jī)因子: 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)

結(jié)果:

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. 遺傳力的計(jì)算

這里 Vg 是Line的方差組分: 0.191933 VGL 是品種與地點(diǎn)的交互方差組分: 0.009483 VGY 為品種與年份的交互方差組分: 0.010209 Ve 為殘差方差組分: 0.257820

這里, 品種與地點(diǎn)以及年份互作, 沒有考慮, 因此不做計(jì)算.

地點(diǎn)為2, 年份為2, 重復(fù)為2

> 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. 計(jì)算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. 對(duì)比BLUP值和平均值

可以看出, BLUP值和平均值趨勢(shì)基本一致, 但是有個(gè)別品種, BLUP值和平均值變化較大.

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

10. 不足

這篇無疑是開山之作, 但是也有一些不足:

一般來說, 多年多點(diǎn)分析中, 我們將地點(diǎn), 年份, 地點(diǎn):年份, 地點(diǎn):年份:重復(fù)作為規(guī)定因子, 品種, 品種與地點(diǎn), 品種與年份, 品種與地點(diǎn)與年份作為隨機(jī)因子. 文章中將所有因子都作為隨機(jī)因子, 太過武斷.
遺傳力公式有錯(cuò)誤:這里沒有考慮: 品種:年份:地點(diǎn), 殘差的分母應(yīng)該是2*2*2 = 8 , 而不是4.
因素沒有考慮完整, 可能是數(shù)據(jù)量有限, 沒有考慮地點(diǎn):年份:重復(fù), 沒有考慮地點(diǎn):年份:品種
計(jì)算遺傳力沒有標(biāo)準(zhǔn)誤, 標(biāo)準(zhǔn)誤可以反映出計(jì)算的好壞.
沒有給出隨機(jī)因子的顯著性, 可以使用LRT檢驗(yàn)

11. 商業(yè)軟件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