做基因组选择Genomic Selection, GS的朋友应该都经历过这种场面基因型数据拿到了十几个性状跑一轮GBLUP一出结果就停在那里想试试其他算法又搞不清它们之间的差别最后还是只写上“使用GBLUP进行预测”。尤其在准备文章的时候审稿人一句“为什么只比较一种算法”就会让整个模型部分显得单薄。这篇内容想把这几年我用过的7个主流GS算法彻底讲透——包括它们各自的统计假设、适合什么样的遗传架构、数据量不一样的时候怎么选以及最关键的一步直接套用一份自己的数据从基因型矩阵跑到交叉验证精度再把方法和结果写成可以直接放进论文里的表述。读完你会发现算法选型并没有那么玄但确实需要一点对底层逻辑的理解。1. 选对算法之前先看懂GS预测精度波动的根源很多人有一个误解觉得“算法越复杂预测精度越高”。这句话在GS领域很大概率是错的。预测精度高不高核心取决于你的统计模型对标记效应的假设跟你做的是不是“Bayes”开头没有必然关系。基因组选择的基本逻辑是用覆盖全基因组的SNP标记去捕捉影响性状的QTL。因为标记足够密任何一个QTL都和某些标记存在连锁不平衡LD所以不需要预先知道哪些标记重要只要把全部标记放进模型就能间接估计出每个个体的基因组育种值GEBV。问题就在于“把全部标记放进模型”这件事不同算法对标记效应的先验分布有完全不同的假设。而这个假设跟真实遗传架构匹不匹配直接决定了预测效果。举个例子。产量、生育期这类性状通常由成千上万个微效基因控制效应分布比较均匀。这种时候RR-BLUP和GBLUP那种“所有标记效应方差相等”的假设就很稳。但如果你研究抗病性性状背后可能只有两三个主效基因加一堆微效修饰基因。要是还用“所有标记一视同仁”的模型主效信号就会被稀释在上万个标记里反而拉低预测精度。这时候BayesB、BayesCπ这些能做变量选择的算法就有优势——它们主动让一部分标记效应归零把信号留给真正相关的位点。另一个影响精度的重要变量是样本量。GS模型的标记数常常几万甚至几十万而样本量可能只有几百。高维小样本局面下复杂模型容易过拟合训练集精度很好验证集上一塌糊涂。所以当你手里只有150个个体的时候上来就跑BayesB大概率得不到什么惊喜反而GBLUP这种正则化十足的方法最稳。这不是算法不行是数据量喂不饱复杂模型的参数空间。理解了这层关系再去看算法之间的差异条理就清楚多了。2. 7种主流GS算法拆解先验假设决定行为2.1 RR-BLUP标记效应全相同方差的岭回归RR-BLUP全称Ridge Regression Best Linear Unbiased Prediction本质就是全标记岭回归。它假设每一个SNP的效应都服从一个均值为0、方差相同的正态分布换句话说所有标记的地位完全平等谁也别想搞特殊。这个假设对应的是“微效多基因”模型。优点是正则化强、计算非常快即使标记有几万几十万个也一样跑得动而且不需要做变量筛选把所有标记扔进去就行。在样本量和标记数比例悬殊的情况下RR-BLUP通常不会翻车。它估计出来的标记效应可以直接用来算个体GEBV也可以顺便看一下哪些标记效应比较大——虽然在这种情况下每个标记的效应都被压得很均匀定位精度一般。不过RR-BLUP对“少数基因控制大效应”的情况不友好。假如某个主效QTL的效应是100模型会硬生生把100拆成1万份平均分到1万个跟它存在LD的标记上。虽然预测GEBV时总和还是对的但单个标记的效应大小没有生物学上的直观意义。2.2 GBLUP从标记直接到亲缘矩阵GBLUPGenomic Best Linear Unbiased Prediction是动物育种上用得最广泛的模型。它的思路是绕开“标记效应”直接用标记构建个体之间的基因组亲缘关系矩阵G矩阵然后估计每个个体的随机育种值。模型写成 y Xb Zu e其中u的方差结构里放G矩阵而不是传统A矩阵。GBLUP和RR-BLUP在数学上是等价的。只要G矩阵的定义合理两者预测出的GEBV结果几乎一致。区别在于GBLUP不输出标记效应所以如果你想顺便找一找“哪些SNP贡献大”它帮不上忙。但它的计算效率非常高因为只需要处理一个个体的n×n矩阵而不用管几十万个标记。对于样本量几千、标记上百万的数据GBLUP是首选。日常操作中我用rrBLUP包的A.mat()函数从基因型矩阵构建G矩阵再用mixed.solve()做一次混合模型求解非常简单。2.3 BayesA给每个标记单独开一个方差账户BayesA开始进入“贝叶斯家族”。它的核心变化是每个SNP效应不再共享同一个方差而是各自拥有独立的方差每个方差的先验服从尺度化逆卡方分布。这意味着模型允许部分标记效应比较大部分标记效应比较小——相当于给每个标记“单独开账户”。这个设计比RR-BLUP灵活得多对混合遗传架构少数大效应多数微效应的适应能力更强。但BayesA没有真正的变量选择机制所有标记效应都不会归零只是大的大、小的小。它也是最基础的贝叶斯GS模型BayesB、BayesC都是对它做加法得到的。实现上BayesA在R的BGLR包里一行代码就能跑但MCMC的收敛速度比RR-BLUP慢得多。运行时间在样本量和标记数都大的时候会明显增加这是你要提前做好心理准备的地方。2.4 BayesB让一部分标记效应“强制归零”BayesB在BayesA的基础上加了一个关键的概率参数π表示每个SNP有π的概率效应为0有1-π的概率服从一个厚尾分布。这个设定天然带有变量选择功能模型会主动把一大批标记的效应压缩成0只保留少数“有关系”的标记。对有主效QTL的性状BayesB通常比RR-BLUP和GBLUP表现出更高的预测精度。但代价也很明显MCMC采样极慢在标记数量特别大时拟合一条链可能要跑几小时甚至一整天。早期实现中π是提前给定的一般设成0.95或0.99如果你的数据里QTL个数没那么多这个默认值还算合理但如果QTL比例很高指定一个过大的π反而会损失精度。2.5 BayesCπ共同方差加上动态πBayesCπ可以看成BayesB的“稳定性改良版”。BayesB允许每个非零效应标记有自己的方差信息会被高频标记分走BayesCπ则假设所有非零标记效应共享同一个共同方差只保留π决定哪些标记进入模型。更友好的一点是π不再是你拍脑袋给的常数而是作为未知参数由MCMC从数据中估计。实际使用中BayesCπ在我遇到的多数中等到高遗传力性状上都表现不错。它兼顾了变量选择和参数估计的稳定性而且对π先验不敏感是贝叶斯方法里比较省心的一款。在BGLR包里model参数直接填BayesC即可它实现的就是包含π估计的BayesCπ模型。2.6 Bayesian LASSOL1惩罚的贝叶斯版本Bayesian LASSOBL本质上是L1正则化的贝叶斯版本。它给每个SNP效应套了一个拉普拉斯先验效果类似把大部分效应往0方向压缩。和BayesB“硬归零”不同BL把效应压缩到接近0但不会精确成为0属于软收缩。这个特性让BL对异常值比较稳健在中等遗传力的复杂性状上往往表现不错收敛性也比BayesB顺畅。不过它没有真正的变量选择机制标记效应不会变成严格的0所以在“找出重要标记”这类需求上帮助有限。如果目的只是预测GEBVBL是个值得加入比较阵营的模型。2.7 BayesRSNP效应分成几个档次的混合模型BayesR是这几年比较热门的模型。它不是逐個标记估计效应的连续分布而是假设每个SNP效应来自几个不同方差的正态分布混合。以最常见的四分量版为例一部分SNP效应为0一部分效应小一部分中一部分大混合比例由数据自适应估计。BayesR的特点是灵活它能自动适配从“完全微效多基因”到“存在大效应QTL”的各种遗传架构。因为它把效应分档既不会像RR-BLUP那样把大效应均匀稀释也不会像BayesB那样对π过度敏感。像GCTB软件里实现的BayesR在人与动物基因组选择研究中应用越来越多。代价同样是计算量大不过GCTB用了一些加速策略实际体验比纯R跑BayesB还要快一点。2.8 七个算法的选择速查表算法标记效应假设是否做变量选择速度最适场景RR-BLUP所有标记同方差正态否快微效多基因、样本量小、快速基线GBLUP由G矩阵间接建模个体效应否很快大规模个体、不做标记定位BayesA每个标记独立方差否中混合遗传架构、标记数适中BayesB非零效应各带方差π归零是慢存在主效QTL的性状BayesCπ非零效应共享方差π数据估计是中慢主效QTL微效基因混合、稳健首选BL拉普拉斯先验软压缩软压缩中中等遗传力复杂性状BayesR效应分多个正态分布档次是分档归零中慢各种遗传架构、大样本这张表只是推荐起点具体数据里谁最优一定要跑交叉验证说话。跑之前记住一句话先快后慢。先用RR-BLUP和GBLUP摸清基线再上贝叶斯家族能省掉很多不必要的等待时间。3. 直接用你的数据跑通全流程R代码实战3.1 基因型数据和表型的准备进入代码前数据格式必须统一。我通常准备两个文件一个是基因型矩阵X行是样本列是SNP标记编码采用0/1/2方式0为AA纯合、1为AB杂合、2为BB纯合缺失值最好先做填补另一个是表型向量y建议用经过固定效应校正后的BLUE或BLUP值而不是原始观测值。如果直接把年份、地点这些环境因素丢进模型不是不行但BGLR里要为固定效应单独建ETA列表写起来会多一层。基因型质控在正式建模前完成最基础的三道门槛MAF小于0.05的标记删掉这些低频标记的效应极不稳定缺失率大于0.1的样和标记删掉或先填补严重偏离哈迪温伯格平衡的标记删掉通常是分型错误的信号。质控可以放在PLINK里做也可以用R完成。如果数据行数比较多建议用PLINK速度快很多。质控后剩下的标记数直接决定后续模型运行时间尤其是贝叶斯方法这一步省下来的时间非常可观。3.2 RR-BLUP和GBLUP的快速实现RR-BLUP和GBLUP使用rrBLUP包就能一次跑完。library(rrBLUP) # X: n行m列基因型矩阵0/1/2编码 # y: n长度表型向量 X - as.matrix(read.table(genotype.txt, header TRUE, row.names 1)) y - read.table(phenotype.txt, header TRUE)$blup # RR-BLUP: 直接估计标记效应 ans_rr - mixed.solve(y y, X X, K NULL) # ans_rr$u 为标记效应ans_rr$beta 为截距 # 训练群体个体GEBV X %*% ans_rr$u ans_rr$beta # GBLUP: 用A.mat构建亲缘关系矩阵再求解BLUP K - A.mat(X) # 默认按VanRaden方法计算 ans_gb - mixed.solve(y y, K K) # ans_gb$u 为每个个体的GEBV对这组代码解释几句。mixed.solve估计的是方差分量和BLUP解整个计算过程不需要做交叉验证时跑得飞快。A.mat()默认使用VanRaden方法构建G矩阵它的对角线和离对角都表示个体之间的基因组亲缘关系。GBLUP和RR-BLUP在这个框架下会给出几乎一致的GEBV排序如果你只想在论文里报告其中一个建议优先写GBLUP因为它更贴近常规BLUP育种模型语言审稿人更熟悉。3.3 用BGLR跑BayesA、BayesB、BayesCπ和BL贝叶斯家族我统一用BGLR包。所有模型共用同一个BGLR()函数只改ETA列表里的model参数。library(BGLR) nIter - 12000 burnIn - 2000 # 设置模型列表 models - c(BayesA, BayesB, BayesC, BL) results - list() for (m in models) { ETA - list(list(X X, model m)) fm - BGLR(y y, ETA ETA, nIter nIter, burnIn burnIn, thin 5) results[[m]] - fm # fm$yHat 是拟合值fm$ETA[[1]]$b 是标记效应后验均值 }跑完之后最常用的检查项是fm$fit的逐迭代数值看看MCMC链的RMSE是否在后期趋于稳定。如果burnIn之后还要很久才能稳定就把nIter加大到20000甚至30000。BGLR里BayesC对应的就是BayesCππ由数据估计不需要手动设。注意BayesB在这个代码下同样完成但它是最慢的一个。如果你的标记有10万以上样本有几百建议先在全部标记里跳过BayesB或先对标记做LD修剪再跑。否则一条链两三个小时很正常而五折交叉验证要跑五遍时间成本直接翻五倍。3.4 跑BayesR的GCTB命令行BayesR在R里没有特别统一易用的实现我通常直接用GCTB软件跑它在处理大数据时效率明显更高。gctb --bfile mydata --bayesR --pheno pheno.txt --chain-length 10000 --burn-in 2000 --out myrun # 输出文件: myrun.mcmc.snp标记效应, myrun.mcmc.samples方差分量采样输入基因型用PLINK二进制格式.bed/.bim/.fam表型文件是两列的文本文件第一列FID第二列IID第三列表型值缺测用-9表示。跑完后从.mcmc.snp文件里取每个标记的后验均值再和X矩阵相乘就能得到个体育种值。GCTB的BayesR还会输出每个SNP属于哪个方差分量的后验概率这非常有用可以帮你识别“关键标记”。3.5 一个完整的5折交叉验证代码框架交叉验证才是GS算法比较的重头戏。BGLR有个很省事的特性y向量里允许出现NA模型会基于后验分布自动预测缺失位置的样本值。这意味着可以把验证集的表型设成NA拟合完成后直接读预测值完全不需要手动实现BLUP预测公式。set.seed(123) n - length(y) nfolds - 5 folds - sample(rep(1:nfolds, length.out n)) models - c(BRR, BayesA, BayesB, BayesC, BL) # 说明BRR是RR-BLUP的贝叶斯等价版本用它做交叉验证等价于RR-BLUP cor_pred - matrix(NA, nrow nfolds, ncol length(models)) colnames(cor_pred) - models for (i in 1:nfolds) { test_idx - which(folds i) y_na - y y_na[test_idx] - NA for (m in models) { ETA - list(list(X X, model m)) fm - BGLR(y y_na, ETA ETA, nIter 6000, burnIn 1000, verbose FALSE) cor_pred[i, m] - cor(fm$yHat[test_idx], y[test_idx]) } } apply(cor_pred, 2, mean) # 每个算法的平均预测能力这个框架有几处需要注意。第一BRR模型和RR-BLUP预测结果高度一致写论文时可以直接说“使用RR-BLUP贝叶斯岭回归形式实现”。第二nIter设成6000只是为了快速看结果正式跑建议提到12000以上。第三cor()要求预测值和观测值都没有缺失如果y里有极端的确定性NA先处理好再跑。第四五折交叉验证的fold划分直接决定结果务必设定随机种子并在论文里写出种子号以便复现。GBLUP的交叉验证也可以放进同一个框架把model换成RKHSeta里提供线性核矩阵K即可。实际操作中我更倾向单独用rrBLUP做GBLUP验证因为A.mat()构建的G矩阵更标准而RKHS在BGLR里需要自己传核函数多一道手续。4. 预测精度计算和交叉验证设计里的坑4.1 表型用什么值BLUE还是BLUP这是很多新手翻车的第一站。如果你拿原始观测表型直接跑GS环境效应没扣除预测精度会被环境噪声拉低而且不同地点、年份的影响会让模型去拟合环境而不是基因型。标准做法是先用混合线性模型校正固定效应得到每个个体的BLUE值或者用单步法从动物模型里提取BLUP值作为响应变量。比如我的一个玉米数据集有3个地点2个年份我习惯先跑一个lme4混合模型library(lme4) fit - lmer(yield ~ (1|genotype) location year (1|location:year), data pheno) blup_vals - ranef(fit)$genotype fixef(fit)[1]然后把blup_vals当作y放进GS模型。注意如果后续要在论文里计算预测精度最好把遗传力估计也用同套数据报告出来因为预测精度的校正公式依赖遗传力。4.2 训练集和验证集的拆分逻辑随机5折交叉验证是最常见的做法但它默认假设验证集个体和训练集个体来自同一个随机交配群体亲缘关系比较接近。如果验证集里有训练集的子女或全同胞预测会变得异常容易得到的精度其实是一种“上限乐观值”。这就是为什么很多高水平文章会强调“预测场景”。要么做随机划分对应“群体内部未知个体的预测”要么做家系水平验证将整个家系同时归入验证集对应“新家系的早期选择”。这两种场景给出的精度含义完全不同。如果审稿人问你“验证集和训练集是否存在亲缘重叠”你必须有意识地回答。我通常建议论文正文用随机交叉验证展示算法比较的结果补充材料里放一个家系水平验证证明模型在新材料中仍然稳健。4.3 预测能力与预测精度怎么报告这里有一个术语陷阱predictive ability和predictive accuracy在不少文献中混用但它们其实不一样。预测能力predictive ability直接就是验证集预测GEBV和校正表型之间的Pearson相关系数通常写为r(y_pred, y_obs)。预测精度predictive accuracy则用r除以遗传力平方根得到acc r / sqrt(h²)。后者在理论上更接近于“真实育种值和估计值之间的相关”数值会偏高一些。如果文章里只能报告一个指标我建议两个都列出来并在表格里明确标注哪个是ability哪个是accuracy。不过要注意如果遗传力估计偏低校正后的精度可能大于1这时就尴尬了。现实中很多实证研究里这个情况并不少见妥帖的做法是在方法部分写清楚遗传力来源并在结果里注释“校正后精度在某些模型中接近或超过1可能是遗传力低估导致”这样就不会被当成错误。4.4 贝叶斯算法的MCMC收敛确认贝叶斯模型跑完了不等于结果可信。BGLR默认设置下nIter6000、burnIn1000这种参数组合对复杂模型是偏短的尤其BayesB这种变量选择模型标记效应的后验分布需要比较长的链才能稳定。判断收敛的土办法是看fm$fit每轮的值在控制台画出迭代轨迹如果前20%还在往下掉说明burnIn不够如果后段仍然大幅波动说明链太短需要加长。更严谨的做法是把同样的模型用不同随机种子跑两条独立链比较两条链在验证集上的预测相关性如果差异在0.01以内基本可以认为MCMC结果稳定了。论文里写MCMC设置时参考这种句式“使用12000次迭代前2000次作为预热保存间隔为5次并用两条独立链验证收敛性”就足够向审稿人交代。5. 发表导向论文里怎么描述GS算法才算专业5.1 材料与方法的可直接套用模板这部分是我觉得最实用的地方。拿到一个新的GS数据做完算法比较后论文的“基因组选择预测”小节照着这个框架改姓名即可采用基因组选择Genomic Selection, GS方法对目标性状的个体育种值进行预测。将全基因组SNP标记纳入统计模型并用以下7种模型进行比较RR-BLUPRidge Regression Best Linear Unbiased Prediction、GBLUPGenomic Best Linear Unbiased Prediction、BayesA、BayesB、BayesCπ、Bayesian LASSO和BayesR。RR-BLUP和GBLUP均假设所有标记效应服从同方差正态分布其中GBLUP通过基因组亲缘关系矩阵G矩阵直接估计个体效应。BayesA假设每个标记具有各自独立的效应方差BayesB在BayesA的基础上引入变量选择概率π以π的概率令标记效应为零BayesCπ将非零标记效应统一分配共同方差并将π作为未知参数由数据估计Bayesian LASSO基于拉普拉斯先验对标记效应进行压缩BayesR假设标记效应来自4个不同方差的正态分布混合比例由数据驱动估计。模型性能采用5折交叉验证重复5次进行评估。每次将全部个体随机划分为5个子集依次以4个子集作为训练集、1个子集作为验证集。预测能力定义为验证集个体预测GEBV与校正表型之间的Pearson相关系数预测精度由预测能力除以性状遗传力平方根获得。所有统计分析在R 4.3.1中完成RR-BLUP和GBLUP使用rrBLUP包实现贝叶斯模型使用BGLR包实现BayesR使用GCTB软件实现。MCMC参数统一设为12000次迭代、前2000次作为burn-in。这段文字把模型、版本、参数、评估指标全说清了审稿人无从挑刺。5.2 结果图和表的推荐呈现方式结果部分没有比“一个大表格一个箱线图”更高效的组合。表格列出每个算法的预测能力均值和标准差再加一列预测精度最后一列放平均运行时间。运行时间虽然看起来不核心但能在无形中告诉审稿人“我确实跑完了这些模型”。箱线图用ggplot2画横轴为算法名称纵轴为预测能力把各个fold的点画上去。这种图可以直观展现哪些算法稳定、哪些算法的fold间波动大。如果做了多个性状建议用分面图把不同性状放在同一行一眼就能看出不同遗传架构下算法排名的变化。5.3 审稿人常问的几个GS算法问题问题一“为什么用5折而不是10折交叉验证”回答要点5折在训练集比例80%和10折比例90%之间预测能力差异通常很小但5折计算量只有10折的一半。在贝叶斯方法耗时较多时5折是效率和可靠性之间的合理折中。问题二“训练集和验证集之间亲缘关系如何”回答要点明确写出随机划分下的个体可能具有不同程度亲缘关系报告训练集与验证集平均亲缘关系值有条件的再配上家系水平验证的结果。问题三“比较不同模型时是否确保标记数、质控标准一致”回答要点必须一致。所有模型使用同一套质控后的标记和同一套交叉验证划分而且随机种子完全相同这样可以保证算法间差异来自模型假设而不是数据处理流程。问题四“遗传力估计使用了什么方法是否影响预测精度校正”回答要点写明遗传力来自全基因组标记基于AI-REML或贝叶斯方差分量的估计并说明校正公式的局限性。这些问题提前在论文里主动交代清楚被审稿人追问的概率会小很多。这几年我自己跑GS数据的一个习惯是拿到数据先做一套完整的质控然后直接上RR-BLUP和GBLUP看基线如果时间充裕再挑2到3个遗传架构差异大的性状跑BayesCπ和BayesR。新增的算法如果只比GBLUP高0.01的精度我也会放进去但会在讨论里说明这个提升的实际意义要结合选择强度来看。真正的价值不在于用多复杂的模型而在于把数据质量、交叉验证设计和生物学解释串起来形成一个读者能信服的故事。算法比较本身只是工具箱里的几个选项别让它们互相打架先想清楚你的性状长什么样再决定让哪个模型上场。 
长河君 · 2026/9/7 23:38:34 · 阅读 1.2k 基因组选择算法全解析：7种主流模型对比与R语言实战
做基因组选择Genomic Selection, GS的朋友应该都经历过这种场面基因型数据拿到了十几个性状跑一轮GBLUP一出结果就停在那里想试试其他算法又搞不清它们之间的差别最后还是只写上“使用GBLUP进行预测”。尤其在准备文章的时候审稿人一句“为什么只比较一种算法”就会让整个模型部分显得单薄。这篇内容想把这几年我用过的7个主流GS算法彻底讲透——包括它们各自的统计假设、适合什么样的遗传架构、数据量不一样的时候怎么选以及最关键的一步直接套用一份自己的数据从基因型矩阵跑到交叉验证精度再把方法和结果写成可以直接放进论文里的表述。读完你会发现算法选型并没有那么玄但确实需要一点对底层逻辑的理解。1. 选对算法之前先看懂GS预测精度波动的根源很多人有一个误解觉得“算法越复杂预测精度越高”。这句话在GS领域很大概率是错的。预测精度高不高核心取决于你的统计模型对标记效应的假设跟你做的是不是“Bayes”开头没有必然关系。基因组选择的基本逻辑是用覆盖全基因组的SNP标记去捕捉影响性状的QTL。因为标记足够密任何一个QTL都和某些标记存在连锁不平衡LD所以不需要预先知道哪些标记重要只要把全部标记放进模型就能间接估计出每个个体的基因组育种值GEBV。问题就在于“把全部标记放进模型”这件事不同算法对标记效应的先验分布有完全不同的假设。而这个假设跟真实遗传架构匹不匹配直接决定了预测效果。举个例子。产量、生育期这类性状通常由成千上万个微效基因控制效应分布比较均匀。这种时候RR-BLUP和GBLUP那种“所有标记效应方差相等”的假设就很稳。但如果你研究抗病性性状背后可能只有两三个主效基因加一堆微效修饰基因。要是还用“所有标记一视同仁”的模型主效信号就会被稀释在上万个标记里反而拉低预测精度。这时候BayesB、BayesCπ这些能做变量选择的算法就有优势——它们主动让一部分标记效应归零把信号留给真正相关的位点。另一个影响精度的重要变量是样本量。GS模型的标记数常常几万甚至几十万而样本量可能只有几百。高维小样本局面下复杂模型容易过拟合训练集精度很好验证集上一塌糊涂。所以当你手里只有150个个体的时候上来就跑BayesB大概率得不到什么惊喜反而GBLUP这种正则化十足的方法最稳。这不是算法不行是数据量喂不饱复杂模型的参数空间。理解了这层关系再去看算法之间的差异条理就清楚多了。2. 7种主流GS算法拆解先验假设决定行为2.1 RR-BLUP标记效应全相同方差的岭回归RR-BLUP全称Ridge Regression Best Linear Unbiased Prediction本质就是全标记岭回归。它假设每一个SNP的效应都服从一个均值为0、方差相同的正态分布换句话说所有标记的地位完全平等谁也别想搞特殊。这个假设对应的是“微效多基因”模型。优点是正则化强、计算非常快即使标记有几万几十万个也一样跑得动而且不需要做变量筛选把所有标记扔进去就行。在样本量和标记数比例悬殊的情况下RR-BLUP通常不会翻车。它估计出来的标记效应可以直接用来算个体GEBV也可以顺便看一下哪些标记效应比较大——虽然在这种情况下每个标记的效应都被压得很均匀定位精度一般。不过RR-BLUP对“少数基因控制大效应”的情况不友好。假如某个主效QTL的效应是100模型会硬生生把100拆成1万份平均分到1万个跟它存在LD的标记上。虽然预测GEBV时总和还是对的但单个标记的效应大小没有生物学上的直观意义。2.2 GBLUP从标记直接到亲缘矩阵GBLUPGenomic Best Linear Unbiased Prediction是动物育种上用得最广泛的模型。它的思路是绕开“标记效应”直接用标记构建个体之间的基因组亲缘关系矩阵G矩阵然后估计每个个体的随机育种值。模型写成 y Xb Zu e其中u的方差结构里放G矩阵而不是传统A矩阵。GBLUP和RR-BLUP在数学上是等价的。只要G矩阵的定义合理两者预测出的GEBV结果几乎一致。区别在于GBLUP不输出标记效应所以如果你想顺便找一找“哪些SNP贡献大”它帮不上忙。但它的计算效率非常高因为只需要处理一个个体的n×n矩阵而不用管几十万个标记。对于样本量几千、标记上百万的数据GBLUP是首选。日常操作中我用rrBLUP包的A.mat()函数从基因型矩阵构建G矩阵再用mixed.solve()做一次混合模型求解非常简单。2.3 BayesA给每个标记单独开一个方差账户BayesA开始进入“贝叶斯家族”。它的核心变化是每个SNP效应不再共享同一个方差而是各自拥有独立的方差每个方差的先验服从尺度化逆卡方分布。这意味着模型允许部分标记效应比较大部分标记效应比较小——相当于给每个标记“单独开账户”。这个设计比RR-BLUP灵活得多对混合遗传架构少数大效应多数微效应的适应能力更强。但BayesA没有真正的变量选择机制所有标记效应都不会归零只是大的大、小的小。它也是最基础的贝叶斯GS模型BayesB、BayesC都是对它做加法得到的。实现上BayesA在R的BGLR包里一行代码就能跑但MCMC的收敛速度比RR-BLUP慢得多。运行时间在样本量和标记数都大的时候会明显增加这是你要提前做好心理准备的地方。2.4 BayesB让一部分标记效应“强制归零”BayesB在BayesA的基础上加了一个关键的概率参数π表示每个SNP有π的概率效应为0有1-π的概率服从一个厚尾分布。这个设定天然带有变量选择功能模型会主动把一大批标记的效应压缩成0只保留少数“有关系”的标记。对有主效QTL的性状BayesB通常比RR-BLUP和GBLUP表现出更高的预测精度。但代价也很明显MCMC采样极慢在标记数量特别大时拟合一条链可能要跑几小时甚至一整天。早期实现中π是提前给定的一般设成0.95或0.99如果你的数据里QTL个数没那么多这个默认值还算合理但如果QTL比例很高指定一个过大的π反而会损失精度。2.5 BayesCπ共同方差加上动态πBayesCπ可以看成BayesB的“稳定性改良版”。BayesB允许每个非零效应标记有自己的方差信息会被高频标记分走BayesCπ则假设所有非零标记效应共享同一个共同方差只保留π决定哪些标记进入模型。更友好的一点是π不再是你拍脑袋给的常数而是作为未知参数由MCMC从数据中估计。实际使用中BayesCπ在我遇到的多数中等到高遗传力性状上都表现不错。它兼顾了变量选择和参数估计的稳定性而且对π先验不敏感是贝叶斯方法里比较省心的一款。在BGLR包里model参数直接填BayesC即可它实现的就是包含π估计的BayesCπ模型。2.6 Bayesian LASSOL1惩罚的贝叶斯版本Bayesian LASSOBL本质上是L1正则化的贝叶斯版本。它给每个SNP效应套了一个拉普拉斯先验效果类似把大部分效应往0方向压缩。和BayesB“硬归零”不同BL把效应压缩到接近0但不会精确成为0属于软收缩。这个特性让BL对异常值比较稳健在中等遗传力的复杂性状上往往表现不错收敛性也比BayesB顺畅。不过它没有真正的变量选择机制标记效应不会变成严格的0所以在“找出重要标记”这类需求上帮助有限。如果目的只是预测GEBVBL是个值得加入比较阵营的模型。2.7 BayesRSNP效应分成几个档次的混合模型BayesR是这几年比较热门的模型。它不是逐個标记估计效应的连续分布而是假设每个SNP效应来自几个不同方差的正态分布混合。以最常见的四分量版为例一部分SNP效应为0一部分效应小一部分中一部分大混合比例由数据自适应估计。BayesR的特点是灵活它能自动适配从“完全微效多基因”到“存在大效应QTL”的各种遗传架构。因为它把效应分档既不会像RR-BLUP那样把大效应均匀稀释也不会像BayesB那样对π过度敏感。像GCTB软件里实现的BayesR在人与动物基因组选择研究中应用越来越多。代价同样是计算量大不过GCTB用了一些加速策略实际体验比纯R跑BayesB还要快一点。2.8 七个算法的选择速查表算法标记效应假设是否做变量选择速度最适场景RR-BLUP所有标记同方差正态否快微效多基因、样本量小、快速基线GBLUP由G矩阵间接建模个体效应否很快大规模个体、不做标记定位BayesA每个标记独立方差否中混合遗传架构、标记数适中BayesB非零效应各带方差π归零是慢存在主效QTL的性状BayesCπ非零效应共享方差π数据估计是中慢主效QTL微效基因混合、稳健首选BL拉普拉斯先验软压缩软压缩中中等遗传力复杂性状BayesR效应分多个正态分布档次是分档归零中慢各种遗传架构、大样本这张表只是推荐起点具体数据里谁最优一定要跑交叉验证说话。跑之前记住一句话先快后慢。先用RR-BLUP和GBLUP摸清基线再上贝叶斯家族能省掉很多不必要的等待时间。3. 直接用你的数据跑通全流程R代码实战3.1 基因型数据和表型的准备进入代码前数据格式必须统一。我通常准备两个文件一个是基因型矩阵X行是样本列是SNP标记编码采用0/1/2方式0为AA纯合、1为AB杂合、2为BB纯合缺失值最好先做填补另一个是表型向量y建议用经过固定效应校正后的BLUE或BLUP值而不是原始观测值。如果直接把年份、地点这些环境因素丢进模型不是不行但BGLR里要为固定效应单独建ETA列表写起来会多一层。基因型质控在正式建模前完成最基础的三道门槛MAF小于0.05的标记删掉这些低频标记的效应极不稳定缺失率大于0.1的样和标记删掉或先填补严重偏离哈迪温伯格平衡的标记删掉通常是分型错误的信号。质控可以放在PLINK里做也可以用R完成。如果数据行数比较多建议用PLINK速度快很多。质控后剩下的标记数直接决定后续模型运行时间尤其是贝叶斯方法这一步省下来的时间非常可观。3.2 RR-BLUP和GBLUP的快速实现RR-BLUP和GBLUP使用rrBLUP包就能一次跑完。library(rrBLUP) # X: n行m列基因型矩阵0/1/2编码 # y: n长度表型向量 X - as.matrix(read.table(genotype.txt, header TRUE, row.names 1)) y - read.table(phenotype.txt, header TRUE)$blup # RR-BLUP: 直接估计标记效应 ans_rr - mixed.solve(y y, X X, K NULL) # ans_rr$u 为标记效应ans_rr$beta 为截距 # 训练群体个体GEBV X %*% ans_rr$u ans_rr$beta # GBLUP: 用A.mat构建亲缘关系矩阵再求解BLUP K - A.mat(X) # 默认按VanRaden方法计算 ans_gb - mixed.solve(y y, K K) # ans_gb$u 为每个个体的GEBV对这组代码解释几句。mixed.solve估计的是方差分量和BLUP解整个计算过程不需要做交叉验证时跑得飞快。A.mat()默认使用VanRaden方法构建G矩阵它的对角线和离对角都表示个体之间的基因组亲缘关系。GBLUP和RR-BLUP在这个框架下会给出几乎一致的GEBV排序如果你只想在论文里报告其中一个建议优先写GBLUP因为它更贴近常规BLUP育种模型语言审稿人更熟悉。3.3 用BGLR跑BayesA、BayesB、BayesCπ和BL贝叶斯家族我统一用BGLR包。所有模型共用同一个BGLR()函数只改ETA列表里的model参数。library(BGLR) nIter - 12000 burnIn - 2000 # 设置模型列表 models - c(BayesA, BayesB, BayesC, BL) results - list() for (m in models) { ETA - list(list(X X, model m)) fm - BGLR(y y, ETA ETA, nIter nIter, burnIn burnIn, thin 5) results[[m]] - fm # fm$yHat 是拟合值fm$ETA[[1]]$b 是标记效应后验均值 }跑完之后最常用的检查项是fm$fit的逐迭代数值看看MCMC链的RMSE是否在后期趋于稳定。如果burnIn之后还要很久才能稳定就把nIter加大到20000甚至30000。BGLR里BayesC对应的就是BayesCππ由数据估计不需要手动设。注意BayesB在这个代码下同样完成但它是最慢的一个。如果你的标记有10万以上样本有几百建议先在全部标记里跳过BayesB或先对标记做LD修剪再跑。否则一条链两三个小时很正常而五折交叉验证要跑五遍时间成本直接翻五倍。3.4 跑BayesR的GCTB命令行BayesR在R里没有特别统一易用的实现我通常直接用GCTB软件跑它在处理大数据时效率明显更高。gctb --bfile mydata --bayesR --pheno pheno.txt --chain-length 10000 --burn-in 2000 --out myrun # 输出文件: myrun.mcmc.snp标记效应, myrun.mcmc.samples方差分量采样输入基因型用PLINK二进制格式.bed/.bim/.fam表型文件是两列的文本文件第一列FID第二列IID第三列表型值缺测用-9表示。跑完后从.mcmc.snp文件里取每个标记的后验均值再和X矩阵相乘就能得到个体育种值。GCTB的BayesR还会输出每个SNP属于哪个方差分量的后验概率这非常有用可以帮你识别“关键标记”。3.5 一个完整的5折交叉验证代码框架交叉验证才是GS算法比较的重头戏。BGLR有个很省事的特性y向量里允许出现NA模型会基于后验分布自动预测缺失位置的样本值。这意味着可以把验证集的表型设成NA拟合完成后直接读预测值完全不需要手动实现BLUP预测公式。set.seed(123) n - length(y) nfolds - 5 folds - sample(rep(1:nfolds, length.out n)) models - c(BRR, BayesA, BayesB, BayesC, BL) # 说明BRR是RR-BLUP的贝叶斯等价版本用它做交叉验证等价于RR-BLUP cor_pred - matrix(NA, nrow nfolds, ncol length(models)) colnames(cor_pred) - models for (i in 1:nfolds) { test_idx - which(folds i) y_na - y y_na[test_idx] - NA for (m in models) { ETA - list(list(X X, model m)) fm - BGLR(y y_na, ETA ETA, nIter 6000, burnIn 1000, verbose FALSE) cor_pred[i, m] - cor(fm$yHat[test_idx], y[test_idx]) } } apply(cor_pred, 2, mean) # 每个算法的平均预测能力这个框架有几处需要注意。第一BRR模型和RR-BLUP预测结果高度一致写论文时可以直接说“使用RR-BLUP贝叶斯岭回归形式实现”。第二nIter设成6000只是为了快速看结果正式跑建议提到12000以上。第三cor()要求预测值和观测值都没有缺失如果y里有极端的确定性NA先处理好再跑。第四五折交叉验证的fold划分直接决定结果务必设定随机种子并在论文里写出种子号以便复现。GBLUP的交叉验证也可以放进同一个框架把model换成RKHSeta里提供线性核矩阵K即可。实际操作中我更倾向单独用rrBLUP做GBLUP验证因为A.mat()构建的G矩阵更标准而RKHS在BGLR里需要自己传核函数多一道手续。4. 预测精度计算和交叉验证设计里的坑4.1 表型用什么值BLUE还是BLUP这是很多新手翻车的第一站。如果你拿原始观测表型直接跑GS环境效应没扣除预测精度会被环境噪声拉低而且不同地点、年份的影响会让模型去拟合环境而不是基因型。标准做法是先用混合线性模型校正固定效应得到每个个体的BLUE值或者用单步法从动物模型里提取BLUP值作为响应变量。比如我的一个玉米数据集有3个地点2个年份我习惯先跑一个lme4混合模型library(lme4) fit - lmer(yield ~ (1|genotype) location year (1|location:year), data pheno) blup_vals - ranef(fit)$genotype fixef(fit)[1]然后把blup_vals当作y放进GS模型。注意如果后续要在论文里计算预测精度最好把遗传力估计也用同套数据报告出来因为预测精度的校正公式依赖遗传力。4.2 训练集和验证集的拆分逻辑随机5折交叉验证是最常见的做法但它默认假设验证集个体和训练集个体来自同一个随机交配群体亲缘关系比较接近。如果验证集里有训练集的子女或全同胞预测会变得异常容易得到的精度其实是一种“上限乐观值”。这就是为什么很多高水平文章会强调“预测场景”。要么做随机划分对应“群体内部未知个体的预测”要么做家系水平验证将整个家系同时归入验证集对应“新家系的早期选择”。这两种场景给出的精度含义完全不同。如果审稿人问你“验证集和训练集是否存在亲缘重叠”你必须有意识地回答。我通常建议论文正文用随机交叉验证展示算法比较的结果补充材料里放一个家系水平验证证明模型在新材料中仍然稳健。4.3 预测能力与预测精度怎么报告这里有一个术语陷阱predictive ability和predictive accuracy在不少文献中混用但它们其实不一样。预测能力predictive ability直接就是验证集预测GEBV和校正表型之间的Pearson相关系数通常写为r(y_pred, y_obs)。预测精度predictive accuracy则用r除以遗传力平方根得到acc r / sqrt(h²)。后者在理论上更接近于“真实育种值和估计值之间的相关”数值会偏高一些。如果文章里只能报告一个指标我建议两个都列出来并在表格里明确标注哪个是ability哪个是accuracy。不过要注意如果遗传力估计偏低校正后的精度可能大于1这时就尴尬了。现实中很多实证研究里这个情况并不少见妥帖的做法是在方法部分写清楚遗传力来源并在结果里注释“校正后精度在某些模型中接近或超过1可能是遗传力低估导致”这样就不会被当成错误。4.4 贝叶斯算法的MCMC收敛确认贝叶斯模型跑完了不等于结果可信。BGLR默认设置下nIter6000、burnIn1000这种参数组合对复杂模型是偏短的尤其BayesB这种变量选择模型标记效应的后验分布需要比较长的链才能稳定。判断收敛的土办法是看fm$fit每轮的值在控制台画出迭代轨迹如果前20%还在往下掉说明burnIn不够如果后段仍然大幅波动说明链太短需要加长。更严谨的做法是把同样的模型用不同随机种子跑两条独立链比较两条链在验证集上的预测相关性如果差异在0.01以内基本可以认为MCMC结果稳定了。论文里写MCMC设置时参考这种句式“使用12000次迭代前2000次作为预热保存间隔为5次并用两条独立链验证收敛性”就足够向审稿人交代。5. 发表导向论文里怎么描述GS算法才算专业5.1 材料与方法的可直接套用模板这部分是我觉得最实用的地方。拿到一个新的GS数据做完算法比较后论文的“基因组选择预测”小节照着这个框架改姓名即可采用基因组选择Genomic Selection, GS方法对目标性状的个体育种值进行预测。将全基因组SNP标记纳入统计模型并用以下7种模型进行比较RR-BLUPRidge Regression Best Linear Unbiased Prediction、GBLUPGenomic Best Linear Unbiased Prediction、BayesA、BayesB、BayesCπ、Bayesian LASSO和BayesR。RR-BLUP和GBLUP均假设所有标记效应服从同方差正态分布其中GBLUP通过基因组亲缘关系矩阵G矩阵直接估计个体效应。BayesA假设每个标记具有各自独立的效应方差BayesB在BayesA的基础上引入变量选择概率π以π的概率令标记效应为零BayesCπ将非零标记效应统一分配共同方差并将π作为未知参数由数据估计Bayesian LASSO基于拉普拉斯先验对标记效应进行压缩BayesR假设标记效应来自4个不同方差的正态分布混合比例由数据驱动估计。模型性能采用5折交叉验证重复5次进行评估。每次将全部个体随机划分为5个子集依次以4个子集作为训练集、1个子集作为验证集。预测能力定义为验证集个体预测GEBV与校正表型之间的Pearson相关系数预测精度由预测能力除以性状遗传力平方根获得。所有统计分析在R 4.3.1中完成RR-BLUP和GBLUP使用rrBLUP包实现贝叶斯模型使用BGLR包实现BayesR使用GCTB软件实现。MCMC参数统一设为12000次迭代、前2000次作为burn-in。这段文字把模型、版本、参数、评估指标全说清了审稿人无从挑刺。5.2 结果图和表的推荐呈现方式结果部分没有比“一个大表格一个箱线图”更高效的组合。表格列出每个算法的预测能力均值和标准差再加一列预测精度最后一列放平均运行时间。运行时间虽然看起来不核心但能在无形中告诉审稿人“我确实跑完了这些模型”。箱线图用ggplot2画横轴为算法名称纵轴为预测能力把各个fold的点画上去。这种图可以直观展现哪些算法稳定、哪些算法的fold间波动大。如果做了多个性状建议用分面图把不同性状放在同一行一眼就能看出不同遗传架构下算法排名的变化。5.3 审稿人常问的几个GS算法问题问题一“为什么用5折而不是10折交叉验证”回答要点5折在训练集比例80%和10折比例90%之间预测能力差异通常很小但5折计算量只有10折的一半。在贝叶斯方法耗时较多时5折是效率和可靠性之间的合理折中。问题二“训练集和验证集之间亲缘关系如何”回答要点明确写出随机划分下的个体可能具有不同程度亲缘关系报告训练集与验证集平均亲缘关系值有条件的再配上家系水平验证的结果。问题三“比较不同模型时是否确保标记数、质控标准一致”回答要点必须一致。所有模型使用同一套质控后的标记和同一套交叉验证划分而且随机种子完全相同这样可以保证算法间差异来自模型假设而不是数据处理流程。问题四“遗传力估计使用了什么方法是否影响预测精度校正”回答要点写明遗传力来自全基因组标记基于AI-REML或贝叶斯方差分量的估计并说明校正公式的局限性。这些问题提前在论文里主动交代清楚被审稿人追问的概率会小很多。这几年我自己跑GS数据的一个习惯是拿到数据先做一套完整的质控然后直接上RR-BLUP和GBLUP看基线如果时间充裕再挑2到3个遗传架构差异大的性状跑BayesCπ和BayesR。新增的算法如果只比GBLUP高0.01的精度我也会放进去但会在讨论里说明这个提升的实际意义要结合选择强度来看。真正的价值不在于用多复杂的模型而在于把数据质量、交叉验证设计和生物学解释串起来形成一个读者能信服的故事。算法比较本身只是工具箱里的几个选项别让它们互相打架先想清楚你的性状长什么样再决定让哪个模型上场。