多年多点的表型处理表型值处理BLUP和BLUE值

2023-10-24 本文已影响0人 wo_monic

当有多个年份和多个地点的表型值之后，基因型只有一套，这时如果用来做GWAS的时候，就需要对表型值进行重新计算。
多年多点的表型值用于GWAS分析前，一般有以下三种方式供预处理表型：

求多年多点表型的均值
使用blup值
使用blue值

BLUP值和BLUE值的区别：

估计随机效应（BLUP）：用于估计随机效应的最佳线性无偏预测值（Best Linear Unbiased Predictor）
估计固定效应（BLUE）：用于估计固定效应的最佳线性无偏估计值（Best Linear Unbiased Estimator）
BLUP值预测品种将来的表现
BLUE值预测品种现在的表现。
BLUP值会对数据进行压缩。
BLUE值会对年份、地点进行矫正，得到的结果和原数据尺度一样。

多年多点无重复的模型

lmer构建的是线性回归模型，我们需要提前知道固定因子和随机因子。
lmer(FL~Taxa+(1|location)+(1|year)+(1|location:year),data=All_fiber_noNA)

上面的公式中，

FL是指定性状的列名，是要预测的性状
Taxa是品种的列名，是固定因子
location是位置的列名，是随机因子
year是年份的列名，是随机因子
All_fiber_noNA是数据框的名称
主的效应是品种、年份、地点、还需要考虑的是品种x年份，品种x地点，年份x地点，品种x年份x地点

多年多点有重复的模型

FL是指定性状的列名，是要预测的性状
Taxa是品种的列名，是固定因子
location是位置的列名，是随机因子
year是年份的列名，是随机因子
Rep是重复，是随机因子
All_fiber_noNA是数据框的名称

使用lme4计算blup值和blue值,这里是多年多点无重复实验

#lme4提供三个函数用于建立模型：线性 ( lmer)、广义线性 ( glmer) 和非线性 ( nlmer.)
mod1 <- lmer(FL~(1|Taxa)+(1|location)+(1|year)+ (1|Taxa:location) + (1|Taxa:year),data=dataframe)
mod2 <- lmer(FL~(1|Taxa)+(1|location)+(1|year)+ (1|Taxa:location) + (1|Taxa:year)+(1|location:year),data=dataframe)
mod3 <- lmer(FL~Taxa+(1|location)+(1|year)+ (1|year:location),data=dataframe)
mod4 <- lmer(FL~Taxa+(1|location)+(1|year)+ (1|year:location)+(1|Taxa:location),data=dataframe)
mod5 <- lmer(FL~Taxa+(1|location)+(1|year)+ (1|year:location)+(1|Taxa:location)+(1|Taxa:year),data=dataframe)
mod6 <- lmer(FL~Taxa+(1|location)+(1|year)+ (1|year:location)+(1|Taxa:location)+(1|Taxa:year)+(1:Taxa:year:location),data=dataframe)
summary(mod1)
summary(mod2)
#计算模型的方差
anova(mod1,mod2,mod3,mod4,mod5,mod6)
# refitting model(s) with ML (instead of REML)
# Data: dataframe
# Models:
#   mod1: FL ~ (1 | Taxa) + (1 | location) + (1 | year) + (1 | Taxa:location) + (1 | Taxa:year)
# mod2: FL ~ (1 | Taxa) + (1 | location) + (1 | year) + (1 | Taxa:location) + (1 | Taxa:year) + (1 | location:year)
# mod3: FL ~ Taxa + (1 | location) + (1 | year) + (1 | year:location)
# mod4: FL ~ Taxa + (1 | location) + (1 | year) + (1 | year:location) + (1 | Taxa:location)
# mod5: FL ~ Taxa + (1 | location) + (1 | year) + (1 | year:location) + (1 | Taxa:location) + (1 | Taxa:year)
# mod6: FL ~ Taxa + (1 | location) + (1 | year) + (1 | year:location) + (1 | Taxa:location) + (1 | Taxa:year) + (1:Taxa:year:location)
# npar   AIC   BIC  logLik deviance  Chisq  Df Pr(>Chisq)    
# mod1    7 17576 17623 -8781.1    17562                          
# mod2    8 16248 16302 -8116.1    16232 1329.9   1     <2e-16 ***
#   mod3  423 15090 17914 -7122.0    14244 1988.2 415     <2e-16 ***
#   mod4  424 15092 17923 -7122.0    14244    0.0   1          1    
# mod5  425 15094 17932 -7122.0    14244    0.0   1          1    
# mod6  425 15094 17932 -7122.0    14244    0.0   0               
# ---
#   Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#显著的是mod2和mod3,所以选择这2个里AIC和BIC值最小的一个，

mod1和2是把Taxa当作随机因子来计算，
mod3,4,5,6是把Taxa当作固定因子来计算。
看到最终是mod2和mod3是达到了极显著水平。所以使用这2个模型进行计算blup和blue值。

评价模型优劣的指标

BIC（Bayesian Information Criterion）和AIC（Akaike Information Criterion）是模型选择中常用的两个信息准则。它们用来衡量统计模型的拟合优度和复杂度，并在比较不同模型时提供一种相对优势判断的标准。
在具体应用中，BIC和AIC的值越小越好。

计算blup值和blue值

使用mod2可以计算blup值

FL_blup <- ranef(mod2)$Taxa
colnames(FL_blup) <- "BLUP"
FL_blup <- FL_blup%>%rownames_to_column("Taxa") #这个是blup值

使用mod3计算blue值

##blue值的计算
FL_blue  <-  as.data.frame(fixef(mod3))
head(FL_blue)
#install.packages("lsmeans")
#植物中，需要使用lsmeans添加截距
library('lsmeans')
re <- lsmeans(mod3,"Taxa")
Blue <- data.frame(re)[,1:2] #这个是Blue值
colnames(Blue) <- c("Taxa","BLUE")
Blup_Blue <- inner_join(FL_blup,Blue,by="Taxa")

使用mod2无法计算blue值，mod3无法计算blup值。
mod2把所有的变量都当作随机因子，就没有固定因子，所以没法算固定效应值（BLUE）。
因为mod3里我把Taxa（品种）当作固定因子，所以无法算品种的随机效应值（BLUP）。

可视化一下最终的数据的分布

hist(Blup_Blue$BLUP,col = "pink")
hist(Blup_Blue$BLUE,col = "orange")

BLUP值的分布

BLUE值的分布

可以看出虽然mod2和mod3的模型不一样，算出的BLUP和BLUE值的分布还是一样的。

参考地址1 https://www.jianshu.com/p/f4f1b5b75830
参考地址2 https://zhuanlan.zhihu.com/p/93495710