混合线性模型中BLUE值 VS BLUP值 计算育种值与GWAS表型性状数据整合
混合线性模型中BLUE值 VS BLUP值
GWAS关联分析课程推荐：https://bdtcd.xet.tech/s/2KgXQq
Best Liner Unbiased Prediction 最佳线性无偏估计
育种值既然重要，计算育种值的方法要配套，最好的方法目前就是BLUP。
BLUP是一种数据分析方法，用的是线性模型方程组，采用的是最小二乘法，目标是离差平方和最小。
你有一些试验数据，很想对这个数据进行分析，找出里面的规律，对数据的分析，有一个办法一直使用的很好，就是回归的思想，或者说成是数据拟合，因为数据间的规律不能确定，对育种这件事来说，我们可以对育种数据直接下部分结论，比如，品种基因型不同，而这种不同会表现在我采集的观测数据中，于是我们可以对育种数据进行效应剖分。
Y=μ + 效应 + σ
这就出来一个线性模型，而这个线性模型是由我的试验数据决定的。
这个模型的含义就是，试验数据里面有效应，有误差，每个观测值里面都包含了效应的误差，当然还有别的。
这时候我们再看原先的试验数据，就可以深入到数据内部了。
如果我们没有模型，这时候进行数据回归，数据拟合是盲目的，没有什么生物学意义支撑，但有了一个准确剖分试验数据的模型后，基于这个模型进行数据拟合，就有意义了。
之所以用线性模型，是因为非线性模型通常有办法化为线性模型，线性模型计算容易，意义好理解，更符合生物学意义，比如，一个品种加入一个产量高的基因，这个品种的产量就高上去了，观测值里面就加进了一部分这个新基因的效应。再加一个新基因，产量就又高一截，产量就层层加上去了，这就是效应可加性。
数据拟合一个巧妙的思路就是，最小二乘，就是离差平方和最小，针对一组试验数据，要用一个模型去进行数学描述，怎样描述才是最完美的呢？就是这个模型可以描述这组试验数据的每个观测值，这有多好啊，任何一个观测值都可以用同一个模型表达出来，可惜实现不了，但是对平均数可以实现这个效果，这也足够完美了。
你的试验数据中某些数据跟一个效应有关，而由于这个效应的存在，导致了观测数据的大小不同，这特定部分的数据可以计算一个平均数出来，确定了这个平均数，就可以确定每个观测值跟这个平均数的离差，离差平方后累加，累加和÷（观测值数目-1），就可以获得方差，这个离差，就是我们要关心的效应对应部分，这是方差分析的思路，可惜我们现在是用线性模型来做，手里只有一大把的未知数，每个观测值就是一个方程，好了，给这个超大方程组，加一个有着明确生物学意义的约束条件，凡是跟这个效应有关的观测值的x, ∑xi=0。
什么意思？跟这个效应有关的所有观测值，这个效应的累加和=0，这样，方程组获得解后，用方程组的解去计算跟这个效应有关的每个观测值，之后累加起来，累积和÷观测值数目，就是一个平均数，就是提到的这个效应的相关观测值的平均数。而每一个xi 就是离差，就是效应，比如，育种值，育种值就是基因效应引发的离差。
BLUP，离差平方和最小，就是用回归的思想，根据线性模型的效应，分解出数据中的各个平均数和离差。
BLUP是个线性方程组，它有多个方程组成，一个观测值就是一个方程。
除了这些方程以外，还有一些方程，比如，你的试验在两个不同的试验地点进行的，所以数据需要区分一个地点效应，因为这是很明确，很有必要的，同一材料在不同地点数据有不同，除了误差外，两者间的一个差异就是地点差异导致的，需要单独分离出来，不要影响数据分析。
Y=μ + 地点效应 b + 育种值效应a + σ
现在你希望的模型是这样的了。
如果一个观测值采集自地点1，我就在方程组里地点1对应的x那里写个1，如果没有就写个0，这样就标示出了观测值跟不同地点间的联系。
这个地点的效应与育种值效应不同，地点效应属于固定效应，育种值效应为随机效应。
怎么理解呢？
试验种在地点1，地点2，我只需要将地点的影响提取出来，对试验分析不要造成干扰，就达到目的了。
而育种值效应，获得试验结果，是为了估计这些试验品种的好坏，这次试验，仅是品种的样本在参加试验，需要通过样本推断总体，本身就具有概率风险，统计上习惯称为随机效应。
线性方程组又要加一个约束条件，就是∑b=0。
跟地点1有关的观测值都包含一个地点1的固定效应b，所有这些b累加和=0。
这样线性方程组就可以计算出地点1的平均值了。
BLUP是最佳线性无偏预测，BLUE是最佳线性无偏估计，
BLUP 对应随机效应，BLUE处理固定效应
BLUP，BLUE分别是两个线性方程组，将他们叠放在一起，组成一个全新的方程组。
由于又有固定效应，又有随机效应，所以就是“混合”了，
因为都是线性方程组，都来自同一个线性模型，于是就叫“线性混合模型”
准确的说法是混合模型方程组（mixed modal equations MME）
这就是Henderson 1948年提出的著名BLUP方法。因为随机效应育种值的预测更有意义，于是BLUP就更出名了，BLUE很少去提。
BLUP方法就是为计算动物育种值才出现的，最初计算机跟不上，解不了大的方程组，使用受限，70年代，BLUP改变了动物育种，如今成为世界公认的最好的动物育种值分析方法。
BLUP表型分析之最优无偏预测R语言分析
最佳线性无偏预测（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
DOI 10.1007/s00122-013-2158-x参考：
http://blog.sina.com.cn/s/blog_13da6d8a60102wcfa.html
https://www.jianshu.com/p/668e3be27530
https://github.com/mighster/Basic_Stats_Scripts/blob/211cb6d47b7308c23e2f84d0dcc42a4a414a4086/IDC_demo.R
- 发表于 2020-11-27 14:54
- 阅读 ( 10249 )
- 分类：GWAS
