||
本节,介绍如何使用R语言的lme4包拟合混合线性模型,计算最佳线性无偏估计(blue)
❝数据来源: Isik F , Holland J , Maltecca C . Genetic Data Analysis for Plant and Animal Breeding. Springer International Publishing, 2017.
❞
「数据及代码下载,请关注公众号:育种数据分析之放飞自我,进入知识星球进行相关下载和学习」
该数据有62个重组自交系(RIL),在4个地点进行试验,随机区组,每个地点2个重复,每个小区种植20株,随机选择5株的表型平均值作为观测值。
library(lme4) library(emmeans) library(data.table) library(tidyverse) library(asreml) dat = fread("MaizeRILs.csv",data.table = F) head(dat) str(dat) col = 1:5 dat[,col] = dat %>% select(all_of(col)) %>% map_df(as.factor) str(dat)
之前,我批量转化为因子,都是用for循环,这次用的是map函数,可读性更强。
这里,使用lme4
包进行blue值计算,然后使用emmeans
包进行预测君子(predict means)的计算,这样就可以将predict means作为表型值进行GWAS分析了。
# lme4 m1 = lmer(height ~ RIL + (1|location:RIL) + (1|location) + (1|location:rep), data=dat) summary(m1) re1 = emmeans(m1,"RIL") %>% as.data.frame() head(re1)
这里,
然后通过emmeans
计算RIL的预测均值。
「注意,lme4直接计算的固定因子(RIL)的效应值(BLUE值),不是我们最终的目的,因为它是效应值,有正有负,我们需要用预测均值将其变为与表型数据尺度一样的水平。」
emmeans这一列就是预测均值了。
library(asreml) m2 = asreml(height ~ RIL, random = ~ location + location:RIL + location:rep,data=dat) summary(m2)$varcomp re2 = predict(m2,classify = "RIL")$pval %>% as.data.frame() head(re2)
这里,用predict进行预测均值的计算,predicted.value这一列就是blue值,两者结果一致。
95%的同学,在计算GWAS分析表型值计算时,都是用上面的模型计算出blue值,然后直接进行计算,其实还有更好的模型。
比如设置每个地点的残差异质,然后和残差同质的模型进行LRT检验,选择最优的模型。
比如设置每个地点与品种的互作的方差异质,比较方差同质的模型,选择最优的模型。
下节见。
「数据及代码下载,请关注公众号:育种数据分析之放飞自我,进入知识星球进行相关下载和学习」
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-9-1 14:40
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社