表型分析之最优无偏预测
最佳线性无偏预测（Best Linear Unbiased Prediction，简称BLUP）可以对多环境数据进行整合，去除环境效应，得到个体稳定遗传的表型。BLUP是表型处理的常用做法。R包lme4中lmer函数是BLUP分析常用的方法，在很多NG文章都引用了该方法。
下面将用实际数据演示多环境无重复数据和多环境有重复数据的过程。首先安装lme4包。
install.packages("lme4")
多环境无重复BLUP
数据格式如下，数据是每个环境叠加的。 有人喜欢用数字表示系名或环境，这样应该把lines和env转换为因子。缺失值用NA表示。
| lines | env | y |
|---|---|---|
| L1 | env1 | 66.72533 |
| L2 | env1 | 53.82899 |
| L3 | env1 | 58.04559 |
| L4 | env1 | 63.09452 |
| L5 | env1 | 57.59054 |
| L6 | env1 | 61.37506 |
library(lme4)
data=read.table("data_blup.txt",header = T)
head(data)
data$lines=factor(data$lines)
data$env=factor(data$env)
接下我们用lmer进行BLUP分析,在lmer中 1|env 表示把env当作随机效应，我们把env和lines当作随机效应。
blp=lmer(y~(1|env)+(1|lines),data=data)
H=5.509/(5.509+3.481/3)
summary(blp)
## Linear mixed model fit by REML ['lmerMod']
## Formula: y ~ (1 | env) + (1 | lines)
## Data: data
##
## REML criterion at convergence: 2700.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.69680 -0.52821 -0.00762 0.58518 2.53832
##
## Random effects:
## Groups Name Variance Std.Dev.
## lines (Intercept) 5.50859 2.3470
## env (Intercept) 0.09091 0.3015
## Residual 3.48151 1.8659
## Number of obs: 575, groups: lines, 209; env, 3
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 60.8465 0.2511 242.3
我们可以得到遗传方差（即lines的方差）和残差方差
。遗传力
.
是遗传方差，
是残差方差，
是环境个数。
我们用ranef函数获取随机效应值,blups返回一个list,包含env和lines的随机效应值即BLUP。blp@beta为整体均值。
blups= ranef(blp)
names(blups)
lines=blups$lines+blp@beta
res=data.frame(id=rownames(lines),blup=lines)
write.table(res,file="data_blup_result.txt",row.names = F,quote = F,sep="\t")
hist(lines[,1],col="#0AB3CA",border="white",xlab="BLUP of lines",main="")
多环境有重复BLUP
数据格式也是叠加的，多了一列rep表示重复。
| env | lines | rep | y |
|---|---|---|---|
| env1 | L1 | rep1 | 373.6640 |
| env1 | L2 | rep1 | 526.4561 |
| env1 | L3 | rep1 | 544.7073 |
| env1 | L4 | rep1 | 602.2171 |
| env1 | L5 | rep1 | 573.5111 |
| env1 | L6 | rep1 | 415.2294 |
library(lme4)
data=read.table("data_rep_blup.txt",header = T)
head(data)
data$lines=factor(data$lines)
data$env=factor(data$env)
data$rep=factor(data$rep)
多环境有重复的分析中，重复rep是嵌套在环境里，rep%in%env 表示嵌套。有重复的数据还可以分析基因型和环境的互作效应。
遗传力公式为。
为基因与环境互作方差，
一个环境的重复数，
为环境个数。
blp=lmer(y~(1|rep%in%env)+(1|env)+(1|lines)+(1|env:lines),data=data,
control=lmerControl(check.nobs.vs.nlev = "ignore",
check.nobs.vs.rankZ = "ignore",
check.nlev.gtr.1 = "ignore",
check.nobs.vs.nRE="ignore"))
H=8154/(8154+9409/3+6121/3/3)
lines=blups$lines+blp@beta
blp.out=data.frame(id=rownames(lines),blp=lines)
write.table(res,file="data_blup_rep_result.txt",row.names = F,quote = F,sep="\t")
summary(blp)
hist(lines[,1],col="#0AB3CA",border="white",xlab="BLUP of lines",main="")
## Linear mixed model fit by REML ['lmerMod']
## Formula:
## y ~ (1 | rep %in% env) + (1 | env) + (1 | lines) + (1 | env:lines)
## Data: data
## Control:
## lmerControl(check.nobs.vs.nlev = "ignore", check.nobs.vs.rankZ = "ignore",
## check.nlev.gtr.1 = "ignore", check.nobs.vs.nRE = "ignore")
##
## REML criterion at convergence: 3754.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6599 -0.3999 0.0071 0.3850 2.8693
##
## Random effects:
## Groups Name Variance Std.Dev.
## env:lines (Intercept) 9409 97.00
## lines (Intercept) 8154 90.30
## env (Intercept) 129802 360.28
## rep %in% env (Intercept) 30449 174.50
## Residual 6121 78.24
## Number of obs: 306, groups:
## env:lines, 102; lines, 34; env, 3; rep %in% env, 1
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 800.0 272.2 2.939
网友评论