用R包sommer做基因组选择：从单性状到多性状GBLUP，一份给育种新手的保姆级代码指南
用R包sommer实现基因组选择：从单性状到多性状GBLUP实战指南
当我在研究生阶段第一次接触基因组选择时，面对复杂的统计模型和编程实现，曾一度感到无从下手。直到发现了R语言的sommer包，这个强大的混合线性模型工具彻底改变了我的研究效率。本文将分享如何用sommer包从零开始构建GBLUP模型，特别适合刚入门基因组选择的育种工作者和农业数据分析师。
1. 环境准备与数据预处理
在开始建模前，我们需要确保环境配置正确并理解数据的基本结构。sommer包的优势在于它能处理复杂的遗传评估模型，但前提是数据必须规范整理。
首先安装并加载必要的包：
install.packages("sommer") library(sommer)sommer包自带的小麦数据集非常适合教学演示：
data(DT_wheat) # 表型数据 data(GT_wheat) # 基因型数据 DT <- DT_wheat GT <- GT_wheat关键预处理步骤：
- 检查数据维度：
dim(DT)显示599个样本的4个性状，dim(GT)显示599×1279的SNP矩阵 - 统一标识符：确保表型和基因型数据的样本ID完全一致
rownames(GT) <- rownames(DT) # 强制统一行名 DT$id <- as.factor(rownames(DT)) # 创建因子型ID列常见问题排查表：
| 问题现象 | 可能原因 | 解决方案 |
|---|---|---|
| 模型报错"维度不匹配" | 样本ID未对齐 | 检查rownames(DT)和rownames(GT) |
| 结果出现NA值 | 基因型数据含缺失 | 用A.mat()前先impute缺失基因型 |
| 模型不收敛 | 数据尺度差异大 | 对表型数据进行标准化 |
提示：始终使用set.seed()固定随机数种子，确保结果可重复。例如：
set.seed(2023)
2. 单性状GBLUP模型构建
单性状GBLUP（ST-GBLUP）是基因组选择的基础模型，我们先从最简单的形式开始。以X1性状为例，随机选取20%个体作为验证集：
vv <- sample(rownames(DT), round(nrow(DT)*0.2)) # 20%验证集 y.trn <- DT y.trn[vv, "X1"] <- NA # 屏蔽验证集表型构建加性关系矩阵和GBLUP模型：
K <- A.mat(GT) # 计算基因组关系矩阵 st_model <- mmer( X1 ~ 1, random = ~ vs(id, Gu = K), rcov = ~ units, data = y.trn, verbose = FALSE )结果提取与解读：
- 育种值提取：
gebv <- st_model$U$u:id$X1 - 遗传力计算：
VG <- st_model$sigma$`u:id`[1,1] # 遗传方差 VE <- st_model$sigma$units[1,1] # 残差方差 h2 <- VG/(VG+VE) # 狭义遗传力- 验证集准确性评估：
cor(gebv[vv,], DT[vv,"X1"], use = "complete") # 预测准确性模型输出关键元素解析：
$sigma - 方差组分(原始尺度) $sigma_scaled - 标准化方差组分 $U - 随机效应预测值(BLUP) $VarU - BLUP方差 $PevU - 预测误差方差 $fitted - 拟合值 $residuals - 残差3. 多性状GBLUP进阶应用
当同时分析多个相关性状时，多性状GBLUP（MT-GBLUP）可以利用性状间的遗传相关提高预测精度。我们联合分析X1和X2性状：
y.trn[vv, c("X1","X2")] <- NA # 同时屏蔽两个性状 mt_model <- mmer( cbind(X1, X2) ~ 1, random = ~ vs(id, Gu = K), rcov = ~ units, data = y.trn, verbose = FALSE )多性状模型特有输出解析：
# 遗传协方差矩阵 G <- mt_model$sigma$`u:id` # 遗传相关 rg <- G[1,2]/sqrt(G[1,1]*G[2,2]) # 表型相关 P <- G + mt_model$sigma$units rp <- P[1,2]/sqrt(P[1,1]*P[2,2])多性状分析优势对比：
- 效率提升：单次分析获得所有性状参数
- 精度提高：利用性状相关增强预测
- 信息整合：直接估计遗传相关
注意：当性状间遗传相关接近0时，多性状模型可能不会带来明显改善
4. 复杂模型扩展与实战技巧
对于更复杂的遗传架构，sommer支持多K矩阵模型。例如将基因组分成两部分分别构建关系矩阵：
# 将SNP分成两部分 half <- floor(ncol(GT)/2) K1 <- A.mat(GT[,1:half]) K2 <- A.mat(GT[,(half+1):ncol(GT)]) # 双K矩阵模型 multiK_model <- mmer( cbind(X1, X2) ~ 1, random = ~ vs(id, Gu = K1) + vs(id, Gu = K2), rcov = ~ units, data = y.trn, verbose = FALSE )高级应用技巧：
- 方差组分约束：通过
constraint参数限制方差为正 - 大规模数据优化：设置
getPEV=FALSE节省内存 - 模型比较：利用AIC/BIC选择最佳模型
- 并行计算：对多性状模型可启用多核加速
常见错误及解决方案：
# 错误1：Gu矩阵维度不匹配 # 解决：检查dim(K)与length(unique(data$id)) # 错误2：模型不收敛 # 解决：调整init参数提供初始值，或检查数据异常值 # 错误3：内存不足 # 解决：对大数据集使用sparse=TRUE选项5. 结果可视化与报告生成
优秀的分析需要专业的可视化呈现。以下是推荐的图形化方法：
育种值分布图：
library(ggplot2) gebv_df <- data.frame( ID = rownames(DT), GEBV = st_model$U$`u:id`$X1, Set = ifelse(rownames(DT) %in% vv, "Validation", "Training") ) ggplot(gebv_df, aes(x=GEBV, fill=Set)) + geom_density(alpha=0.5) + ggtitle("GEBV Distribution")遗传相关矩阵热图：
library(pheatmap) pheatmap(G, display_numbers = TRUE, cluster_rows = FALSE, cluster_cols = FALSE, main = "Genetic Covariance Matrix")实战中发现，将关键结果整理成表格最受育种家欢迎。使用knitr::kable()生成出版级表格：
results <- data.frame( Trait = c("X1", "X2"), h2 = c(h2_X1, h2_X2), Accuracy = c(acc_X1, acc_X2) ) knitr::kable(results, digits = 3, caption = "Genomic Prediction Results")最后提醒，实际应用中要特别注意：
- 基因型质量控制：先进行MAF过滤和缺失率筛选
- 表型数据调整：考虑固定效应（如年份、地点）的影响
- 交叉验证设计：采用k折交叉验证更可靠