多年多年数据如何计算BLUE值

有老师写信给我，询问我如何计算BLUE值，问的人多了，就写一篇博客解释一下。

其实大家来写信，主要是问代码如何写，而我写博客，也是讲代码如何写。

如果对你有帮助，还请多多点赞，转发，十分感谢。

1. 为何要计算BLUE值？

一年多点或者多年多点的植物数据中，一个基因型（品种）往往有多个表型数据，但只有一个基因型，在GWAS关联分析中，就需要一个基因型对应一个表型数据。

之所以有多个表型数据的原因：

• 或者是多个重复
• 或者是多个地点的数据
• 或者是多个年份的数据

问题：如何计算得到一个表型数据呢？

解答：可以使用多个表型值的平均值，作为品种的表型值，现在有更好的方法：BLUE值。

2. 为何使用BLUE值？

一般，有两个选择，BLUE值或者BLUP值，在GWAS中大都使用的BLUE值。

BLUE和BLUP的区别：

• BLUE值是混合线性模型中固定因子的估计效应值
• BLUP值是混合线性模型中随机因子的估计效应值

BLUE和BLUP的代表：

• BLUE值着重在于评估品种现在的表现
• BLUP值着重在于预测品种将来的表现

BLUE和BLUP的方差变化

• BLUE只是对表型值根据地点，年份进行矫正，得到的数据和原来数据尺度一样
• BLUP值会对表型数据进行压缩

3. 示例数据

数据为learnasreml中的MET数据集。数据包括2年，5个地点，每个地点4个重复，共有10品种，观测值为产量（yield）

代码

library(learnasreml)
library(magrittr)
data("MET")
head(MET)

数据

4. lme4包如何分析

模型：

• 固定因子：Cul
• 随机因子：Year + Location + Location:Rep

代码：

dat = MET
m1 = lmer(Yield ~ Cul  + (1|Location) + (1|Location:Rep) + (1|Year), data=dat)
as.data.frame(fixef(m1))

注意：植物中，一般的BLUE值需要加上截距（Intercept）。因为BLUE值中，第一个水平会当做0，其它为相对值，可以手动进行相加，也可以使用lsmeans包中的lsmeans。

library(lsmeans)
re = lsmeans(m1,"Cul") 
re

数据中的lsmeans即为品种的BLUE值，可以作为GWAS或者GS的表型值进行后续的计算。

5. asreml对比结果

众所周知，asreml是一个非常强大的商业软件，如果用asreml进行结果对比，可以判断lme4计算是否正确。

代码

library(asreml)
m2 = asreml(Yield ~ Cul , random = ~Location/Rep + Year, data=dat)
pre = as.data.frame(predict(m2,classify = "Cul")$pvals)
pre

结果

结果中的predicted.value即为品种的BLUE值。

两者结果对比

library(tidyverse)
ls_value = as.data.frame(re)
as_value  = pre

head(ls_value)
head(as_value)

merge(ls_value,as_value,by = "Cul") %>% select(lsmean ,predicted.value) %>% cor

相关性分析结果：

散点图：

可以看出，两者结果完全一致。