表型分析之最优无偏预测
最佳线性无偏预测(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.007620.585182.53832#### Random effects:
## Groups Name Variance Std.Dev.
## lines (Intercept)5.508592.3470
## env (Intercept)0.090910.3015
## Residual 3.481511.8659## Number of obs:575, groups: lines,209; env,3
#### Fixed effects:## Estimate Std.Error t value
##(Intercept)60.84650.2511242.3
我们可以得到遗传方差(即lines的方差)
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)
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.39990.00710.38502.8693
#### Random effects:## Groups Name Variance Std.Dev.## env:lines (Intercept)940997.00
## lines (Intercept)815490.30
## env (Intercept)129802360.28
## rep %in% env (Intercept)30449174.50
## Residual 612178.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.0272.22.939