GWAS计算BLUE值2--LMM计算BLUE值

2024-06-07 15:54:24

本节，介绍如何使用R语言的lme4包拟合混合线性模型，计算最佳线性无偏估计（blue）

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. 读取数据及转换为因子

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函数，可读性更强。

3. 使用lme4包进行blue值计算

这里，使用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)

这里，

RIL作为固定因子
地点和品种互作，作为随机因子
地点内区组，作为随机因子

然后通过emmeans计算RIL的预测均值。

注意，lme4直接计算的固定因子（RIL）的效应值（BLUE值），不是我们最终的目的，因为它是效应值，有正有负，我们需要用预测均值将其变为与表型数据尺度一样的水平。

emmeans这一列就是预测均值了。

4. 使用asreml包进行blue值计算

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值，两者结果一致。

5. 你以为这样就结束了？

95%的同学，在计算GWAS分析表型值计算时，都是用上面的模型计算出blue值，然后直接进行计算，其实还有更好的模型。

比如设置每个地点的残差异质，然后和残差同质的模型进行LRT检验，选择最优的模型。

比如设置每个地点与品种的互作的方差异质，比较方差同质的模型，选择最优的模型。

下节见。

数据及代码下载，请关注公众号：育种数据分析之放飞自我，进入知识星球进行相关下载和学习

GWAS计算BLUE值2--LMM计算BLUE值相关推荐

pandas使用nunique函数计算dataframe每个数据列的独特值的个数（count number of unique values in each column of dataframe）
pandas使用nunique函数计算dataframe每个数据列的独特值的个数(count number of unique values in each column of dataframe) ...
pandas使用groupby函数和cumsum函数计算每个分组内的数值累加值、并生成新的dataframe数据列（ cumulative sum of each group in dataframe
pandas使用groupby函数和cumsum函数计算每个分组内的数值累加值.并生成新的dataframe数据列( cumulative sum of each group in dataframe ...
python尝试不同的随机数进行数据划分、使用卡方检验依次计算不同随机数划分下训练接和测试集所有分类特征的卡方检验的p值，如果所有p值都大于0.05则训练集和测试集都具有统计显著性、数据划分合理
python尝试不同的随机数进行数据划分.使用卡方检验依次计算不同随机数划分下训练接和测试集所有分类特征(categorical)的卡方检验的p值,如果所有p值都大于0.05则退出循环.则训练集和测试 ...
R语言ggplot2可视化：使用dplyr包计算每个分组个数的比例(对计算获得的百分比进行近似，值保留整数部分)、使用ggplot2可视化条形图(bar plot)、并在条形图上添加百分比标签
R语言ggplot2可视化:使用dplyr包计算每个分组个数的比例(对计算获得的百分比进行近似,值保留整数部分).使用ggplot2可视化条形图(bar plot).并在条形图上添加百分比标签目录
pandas使用cut函数基于分位数进行连续值分箱（手动计算分位数）处理后出现NaN值原因及解决
pandas使用cut函数基于分位数进行连续值分箱(手动计算分位数)处理后出现NaN值原因及解决目录 pandas使用cut函数基于分位数进行连续值分箱(手动计算分位数)处理后出现NaN值原因及解决 ...
熵权法中计算的熵值与决策树的熵值完全不一样之谜
熵权法中的熵值计算公式如图所示: 比如说某个评价的指标完全一样,1,1,1,1,1,1 那么m=6,p1到p6的概率均等于1/6:这个时候的熵值是最大的:所以在计算指标权重时,用这种方法反而是数据越小 ...
处理时间_3_计算两个时间列工作日差值
计算两个时间列工作日差值需求描述需求:对EMP表里员工KING和SMITH的hiredate入职时间差,这里单位是天且是工作日时间,即周末不计算在内. 解决方法:通过DATEDIFF函数来完成. ...
处理时间_2_计算两个时间列的差值
计算两个时间列的差值需求描述需求:对EMP表里员工KING和SMITH的hiredate入职时间差,这里单位是分钟.小时.天.周.月.年. 解决方法:通过DATEDIFF函数来完成. 注: 数据库 ...
c语言调用函数计算分段函数值,输入x,计算并输出下列分段函数f(x)的值(保留2位小数) c语言...
计算分段函数输入 x ,计算并输出 y 的值.公式如下当x >= 0时,f(x) = x^0.5,当x小于0时,f(x #include#include//[1{intmain()//你所写的 ...
python中的content方法_content最新:python计算Content-MD5并获取文件的Content-MD5值方式_爱安网 LoveAn.com...
关于"content"的最新内容聚合阅读这篇文章主要介绍了python计算Content-MD5并获取文件的Content-MD5值方式,具有很好的参考价值,希望对大家有所帮助 ...

最新文章

热门文章