GWAS计算BLUE值1--计算最小二乘均值（lsmeans） #2021.12.11
上一次，我计划写个系列，为何？要用BLUE值作表型进行GWAS分析，GWAS分析多年多点或者一年多点的数据时，如何计算矫正后的均值（BLUE值），肝了一上午，写了四篇，从原理到计算方法到代码展示，后面四天的素材就有了，总结一些东西，总能理解更深，保持输出才能不断输入，加油！
本节，介绍如何使用R语言的lm拟合一般线性模型，计算最小二乘均值（lsmeans）
1. 试验数据
❝数据来源：Isik F , Holland J , Maltecca C . Genetic Data Analysis for Plant and Animal Breeding. Springer International Publishing, 2017.❞
该数据有62个重组自交系（RIL），在4个地点进行试验，随机区组，每个地点2个重复，每个小区种植20株，随机选择5株的表型平均值作为观测值。
2. 用R语言进行一般线性模型拟合
拟合模型：
代码语言：javascript复制m1 = lm(height ~ location*RIL location:rep,data=dat)
summary(m1)
anova(m1)
上面是方差分析的结果。
系数的结果是：
注意，这里的值是系数，不是最小二乘均值。这里，如果我们要计算第一个品种RIL1的lsmeans（最小二乘均值），我们需要：
即我们需要整体均值 品种RIL1的回归系数 地点的效应平均值 地点内区组效应品均值 品种RIL1和地点互作的效应品均值
3. 手动计算RIL11的最小二乘均值
这里我们要计算RIL-11这个品种，
「整体均值」
代码语言：javascript复制coef = summary(m1)$coefficients
coef["(Intercept)",1]
可以看到整体均值为179.7973
「地点效应平均值」注意，这里共四个地点，但是只有三个效应值，因为有一个强制为0了，我们在计算平均值时，需要3个地点效应值的和除以4才可以。
代码语言：javascript复制coef[c("locationCLY","locationPPAC","locationTPAC"),1]
sum(coef[c("locationCLY","locationPPAC","locationTPAC"),1])/4
可以看到，效应值为2.466935
「地点内区组效应值」
注意，这里共有4个效应值（4*（2-1）），但是分母应该是4*2=8.
代码语言：javascript复制coef[c("locationARC:rep2","locationPPAC:rep2","locationCLY:rep2","locationTPAC:rep2"),1]
sum(coef[c("locationARC:rep2","locationPPAC:rep2","locationCLY:rep2","locationTPAC:rep2"),1])/8
可以看到，效应值为-1642473
「品种与地点互作的效应」
代码语言：javascript复制coef[c("locationCLY:RILRIL-11","locationPPAC:RILRIL-11","locationTPAC:RILRIL-11"),1]
sum(coef[c("locationCLY:RILRIL-11","locationPPAC:RILRIL-11","locationTPAC:RILRIL-11"),1])/4
loc_ril11 = sum(coef[c("locationCLY:RILRIL-11","locationPPAC:RILRIL-11","locationTPAC:RILRIL-11"),1])/4
可以看到，效应值为-3.425
「RIL-11的效应值为：」
代码语言：javascript复制coef["RILRIL-11",1]
ril11 = coef["RILRIL-11",1]
RIL-11的效应值为4.2
「最后，RIL-11的最小二乘均值为：182.875」
代码语言：javascript复制mu loc loc_rep loc_ril11 ril11
4. 使用函数计算最小二乘均值
之前都是用lsmeans这个包，现在用emmeans，可以看作是lsmeans的升级包。
但是，数据量大时，这个包也是巨慢。
代码语言：javascript复制library(emmeans)
re1 = emmeans(m1,"RIL") %>% as.data.frame()
head(re1,10)
结果是一致的。
5. 用asreml进行估计
代码语言：javascript复制library(asreml)
m2 = asreml(height ~ location*RIL location:rep,data=dat)
summary(m2)$varcomp
re2 = predict(m2,classify = "RIL")$pval %>% as.data.frame()
head(re2)
这里用predict函数，计算RIL的lsmeans：
可以看到结果是一致的。
6. 总结
一般，很少用一般线性模型去估算最小二乘均值，都是用混合线性模型，将品种作为固定因子，去估计BLUE值（最佳线性无偏估计）。
用一般线性模型，演示一下如何计算lsmeans，通过手动计算和函数计算两种形式，理解计算方法。
另外，lsmeans和整体平均值不一样，它比平均值更能代表表型值。所以，如果不使用混合线性模型，使用lsmeans作为表型值，也要比平均值更好。
