用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模型构建单性状GBLUPST-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/(VGVE) # 狭义遗传力验证集准确性评估cor(gebv[vv,], DT[vv,X1], use complete) # 预测准确性模型输出关键元素解析$sigma - 方差组分(原始尺度) $sigma_scaled - 标准化方差组分 $U - 随机效应预测值(BLUP) $VarU - BLUP方差 $PevU - 预测误差方差 $fitted - 拟合值 $residuals - 残差3. 多性状GBLUP进阶应用当同时分析多个相关性状时多性状GBLUPMT-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[,(half1):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参数限制方差为正大规模数据优化设置getPEVFALSE节省内存模型比较利用AIC/BIC选择最佳模型并行计算对多性状模型可启用多核加速常见错误及解决方案# 错误1Gu矩阵维度不匹配 # 解决检查dim(K)与length(unique(data$id)) # 错误2模型不收敛 # 解决调整init参数提供初始值或检查数据异常值 # 错误3内存不足 # 解决对大数据集使用sparseTRUE选项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(xGEBV, fillSet)) geom_density(alpha0.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折交叉验证更可靠