表型分析之最優(yōu)無偏預測
最佳線性無偏預測（Best Linear Unbiased Prediction汛兜，簡稱BLUP）可以對多環(huán)境數(shù)據(jù)進行整合，去除環(huán)境效應膳帕，得到個體穩(wěn)定遺傳的表型辆沦。BLUP是表型處理的常用做法涂邀。R包lme4中l(wèi)mer函數(shù)是BLUP分析常用的方法，在很多NG文章都引用了該方法剿牺。
下面將用實際數(shù)據(jù)演示多環(huán)境無重復數(shù)據(jù)和多環(huán)境有重復數(shù)據(jù)的過程裳擎。首先安裝lme4包。
install.packages("lme4")
多環(huán)境無重復BLUP
數(shù)據(jù)格式如下斯碌，數(shù)據(jù)是每個環(huán)境疊加的一死。 有人喜歡用數(shù)字表示系名或環(huán)境，這樣應該把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的方差）和殘差方差
伪煤。遺傳力
.
是遺傳方差，
是殘差方差凛辣，
是環(huán)境個數(shù)抱既。
我們用ranef函數(shù)獲取隨機效應值,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="")
多環(huán)境有重復BLUP
數(shù)據(jù)格式也是疊加的蝙砌，多了一列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)
多環(huán)境有重復的分析中跋理，重復rep是嵌套在環(huán)境里择克，rep%in%env 表示嵌套。有重復的數(shù)據(jù)還可以分析基因型和環(huán)境的互作效應前普。
遺傳力公式為肚邢。
為基因與環(huán)境互作方差，
一個環(huán)境的重復數(shù)，
為環(huán)境個數(shù)骡湖。
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
