做育种数据分析这几年我越来越觉得遗传相关是个“安安静静决定成败”的参数。表型相关很容易算肉眼也大概看得出来但真正决定我们能不能用一个易测性状去替难测性状做间接选择、能不能在多个目标性状之间找平衡的是那个看不见的遗传相关。前几天我们把一份多环境试验数据的遗传相关跑出来时发现产量和蛋白质含量存在明显的负向遗传相关——这解释了我们前几年为什么高产品系品质总是上不去。这篇文章就从这份工作出发完整演示一遍遗传相关的数据代码分析流程从数据整理、多性状混合模型构建到遗传相关矩阵的计算与解读再把收敛失败、异常估计、标准误这些我踩过的坑一并捋一遍。内容适合正在做植物或动物育种的同行也适合刚接触遗传参数估计、想快速跑通流程的朋友。1. 为什么育种人总绕不开遗传相关1.1 遗传相关是用一个性状替另一个性状做选择的底气遗传相关genetic correlation通常记作 (r_g)指的是两个性状之间由遗传原因引起的相关程度。说得再直白一点如果我们把一批材料的基因型换成另一批两个性状的表型均值会不会发生联动变化如果会这个联动有多强、方向是正是负这就是遗传相关要回答的问题。举个育种里几乎天天遇到的例子在鲜食玉米群体里产量和籽粒含水量往往存在正向遗传相关这意味着你在田间把产量选上去了含水量也可能跟着上去了收获时籽粒偏湿烘干成本上升。再比如小麦的产量和赤霉病抗性很多群体里存在负向遗传相关高产和抗病两件事常常“打架”。如果育种者对这些关系一无所知只按单个性状拼命选等发现目标性状被带偏时往往已经过去了三四个世代材料都进高代品比了非常被动。遗传相关更大的用处在于间接选择。有些目标性状是很难直接测的根系活力要挖根抗逆性要人工逆境鉴定品质指标要近红外扫描或者化学分析成本高、周期长、通量低。如果我们能找到另一个容易测、遗传力高、跟目标性状有较强遗传相关的辅助性状就可以在大群体里只测辅助性状间接对目标性状进行选择。间接选择响应的大致公式是[ CR i \times h_x \times r_g \times \sigma_{Ay} ]其中 (i) 是选择强度(h_x) 是辅助性状遗传力的平方根(r_g) 是辅助性状与目标性状的遗传相关(\sigma_{Ay}) 是目标性状加性遗传标准差。看得出遗传相关直接决定了间接选择的上限。(r_g) 如果是0.3间接选择基本没什么效率如果能到0.7以上就很值得认真考虑了。1.2 表型相关、环境相关与遗传相关的三角关系很多人拿到一份数据先算个皮尔逊相关系数就当作“相关”了这一步本身没错但那是表型相关(r_p)不是遗传相关。表型相关是两个性状在田间实际观测值的相关它是遗传和环境共同作用的结果三者之间有如下关系[ r_p r_g \times h_x \times h_y r_e \times e_x \times e_y ]这里 (h_x)、(h_y) 是两性状遗传力平方根(e_x)、(e_y) 是环境效应占比的平方根(r_e) 是环境相关。这个公式说明一件事表型相关里混着环境相关混着的比例取决于遗传力高低。遗传力越低表型相关越可能被环境相关“稀释”甚至“带偏方向”。我自己就踩过这样的坑。有一年看田间数据发现株高和倒伏率表型相关是正的直观感觉“越高越容易倒”。但用多性状模型拆完遗传相关之后发现两者在遗传层面其实接近零相关真正扯上关系的是水分管理和施肥差异造成的环境效应。如果按表型相关去设计抗倒伏选择方案很可能白忙一场。所以做育种决策尤其是间接选择和指数选择一定以遗传相关为准表型相关只能当作辅助参考。2. 遗传相关估计的方法论多性状混合模型到底在算什么2.1 从单性状到多性状混合模型的一次升级遗传相关没法像算普通相关系数那样把两个性状的观测值列出来直接求比值。原因在于田间试验里的每个观测值都由基因型效应、环境效应、区组效应、残差共同决定普通相关把这些混在一起了。要把遗传那一部分单独抽出来就得用混合模型。先看单性状的基本形式[ y Xb Zg e ](b) 是固定效应比如环境均值、区组效应(g) 是随机遗传效应(e) 是残差。到了两个性状模型扩成矩阵形式[ \begin{bmatrix} y_1 \ y_2 \end{bmatrix} \begin{bmatrix} X_1 0 \ 0 X_2 \end{bmatrix} \begin{bmatrix} b_1 \ b_2 \end{bmatrix} \begin{bmatrix} Z_1 0 \ 0 Z_2 \end{bmatrix} \begin{bmatrix} g_1 \ g_2 \end{bmatrix} \begin{bmatrix} e_1 \ e_2 \end{bmatrix} ]关键在随机效应部分的方差结构。遗传效应 (g_1)、(g_2) 的方差-协方差矩阵是一个 (2\times2) 矩阵[ G \begin{bmatrix} \sigma^2_{g1} \sigma_{g12} \ \sigma_{g12} \sigma^2_{g2} \end{bmatrix} ]这里的 (\sigma_{g12}) 就是两个性状的遗传协方差遗传相关则由这个公式得到[ r_g \frac{\sigma_{g12}}{\sqrt{\sigma^2_{g1} \times \sigma^2_{g2}}} ]所以做遗传相关的第一步不是去写一个“相关性计算函数”而是把一个多性状混合模型跑好估计出G矩阵里的方差和协方差组分再做一次标准化换算。很多人问“为什么遗传相关这么难算”根本原因就在这里它不是一个描述统计量而是一个模型参数。2.2 REML与贝叶斯Gibbs两条主流算法路线方差组分估计的主流路线有两条。第一条是REML即约束极大似然法。它在估计方差组分时先把固定效应对应的自由度“扣掉”再用极大似然去拟合剩余部分这样得到的方差组分受固定效应水平数的影响更小不容易被低估。REML配合BLUP最佳线性无偏预测使用先估方差组分再预测育种值是当前农林领域应用最广的流程。第二条是贝叶斯路线主要用Gibbs抽样。它的思路是为方差组分设定先验分布然后通过马尔可夫链蒙特卡洛抽样得到后验分布。贝叶斯的好处是能直接给每个参数一个完整的后验区间不用另外算标准误在样本量小、模型复杂的时候也比REML更容易通过先验信息把估计“稳住”。代价是计算量大而且先验怎么设定直接影响结果不太适合刚上手的人无脑跑。软件选型上我主力用的是R语言。开源免费这条路里sommer包对植物育种非常友好语法直观支持系谱A矩阵和分子标记G矩阵breedR也做得很扎实空间分析功能强MCMCglmm适合做贝叶斯多性状模型但参数设置门槛更高。商业软件如ASReml也是行业标杆速度和处理超高维数据的能力很强。对于团队联合分析、需要快速迭代的项目我目前最常用的组合是sommer跑REMLMCMCglmm做交叉验证。2.3 遗传力和遗传相关是同一模型产出的孪生参数跑完多性状模型之后我们通常先把遗传力也算出来。单性状遗传力 (h^2 \sigma^2_g / (\sigma^2_g \sigma^2_e))在多性状模型里对应每个性状分别取对角线元素。遗传力和遗传相关出自同一个G矩阵只是标准化方式不同遗传力用每个性状自身的方差做分母遗传相关用两个性状的标准差做分母。这个“孪生”关系提醒我们一件事如果某个性状遗传力低它的遗传方差估计就不稳定分母不稳定遗传相关的估计自然跟着波动。所以我在解读低遗传力性状的遗传相关时特别保守不会只看点估计就说“相关很强”一定要求看标准误。3. 数据代码演示从三份原始表到一张遗传相关矩阵3.1 数据准备表型、系谱、基因型矩阵缺一不可做遗传相关之前先确认手上有哪几类数据。最基础的是表型数据至少包含材料编号geno、环境编号env、区组或重复编号rep、以及至少两个性状的观测值。如果材料分散在多个环境、多个重复数据要按长表格式整理一行一个观测记录。第二类是亲缘关系信息。如果没有分子标记数据可以用系谱记录构建A矩阵如果有全基因组标记就用标记构建G矩阵。A矩阵和G矩阵的作用是描述材料之间的遗传关联亲缘关系越近遗传效应之间的协方差应该越大。这一步不是可有可无因为混合模型里随机遗传效应的方差结构必须乘以亲缘关系矩阵否则模型等同于假设材料间完全独立那就失去了育种数据分析的意义。第三类是可选的协变量信息比如播种密度、土壤肥力等级、采收日期等。这些如果已知对性状有系统性影响就放进固定效应处理。数据大概长这样genoenvrepYield(kg/ha)Protein(%)DH037E11864011.2DH037E12852011.5DH102E11912010.4DH102E21876010.9这份数据来自一个DH系群体试验材料数量不多不少120个DH系3个环境每个环境3次重复属于典型的育种早代材料多环境鉴定数据。数据量不大但足以跑出一个有参考意义的遗传相关。3.2 sommer包多性状模型完整可复现的代码我用R的sommer包来演示。如果你是第一次接触建议先把包装好然后按下面的顺序一步步走。install.packages(sommer) # 如果还没装过 library(sommer) # 1. 读入表型数据 pheno - read.csv(multi_pheno.csv, stringsAsFactors TRUE) # 2. 读入系谱并构建A矩阵 ped - read.csv(three_gen_pedigree.csv, stringsAsFactors TRUE) ped - pedigree_trim(ped, pheno$geno) # 把系谱裁到只保留试验材料及其祖先 A - A.mat(ped, min.common 2) # 基于系谱计算亲缘关系矩阵 # 3. 关键一步保证表型中的材料顺序和A矩阵行名一致 pheno$geno - factor(pheno$geno, levels rownames(A)) # 4. 多性状模型产量和蛋白质含量作为两个响应变量 mod - mmer(cbind(Yield, Protein) ~ Env, random ~ vsr(geno, G A), rcov ~ vsr(units), data pheno) summary(mod)这段代码里cbind(Yield, Protein)告诉模型我们要同时分析两个性状~ Env表示把环境作为固定效应random ~ vsr(geno, G A)指定材料的随机遗传效应并告诉模型用A矩阵定义材料间的协方差结构rcov ~ vsr(units)是残差项sommer会自动把残差展开成2x2的方差协方差结构。每次跑完模型我习惯看一眼输出末尾有没有警告信息。如果出现类似“cannot be fitted”或者“boundary”提示通常意味着数据或模型设置有地方不对绝对不能直接拿去算遗传相关。3.3 从方差组分中算出遗传相关矩阵模型跑通之后方差组分存在mod的对象里。用下面代码提取遗传效应方差协方差矩阵vc - summary(mod)$varcomp gen_comp - mod$sigma$geno print(gen_comp)假设输出长这样我这边跑完的典型结果Yield Protein Yield 456.32 -18.27 Protein -18.27 11.32这里对角线上是各性状的遗传方差副对角线是两性状的遗传协方差。遗传相关由以下代码计算g_var_yield - gen_comp[Yield, Yield] g_var_protein - gen_comp[Protein, Protein] g_cov - gen_comp[Yield, Protein] rg - g_cov / sqrt(g_var_yield * g_var_protein) rg把这组数字代进去会得到约-0.254。也就是说产量和蛋白质含量在这套DH群体里存在中等偏弱的负向遗传相关。这个数值意味着如果只按产量高强度选择蛋白质含量会被缓慢往下带但如果把两者放进指数选择也并非完全卡死因为相关不是强到-0.8那种程度。想得到完整的相关矩阵还可以批量处理# 以遗传标准差对角阵做标准化得到遗传相关矩阵 D - diag(1 / sqrt(diag(gen_comp)), nrow 2) rg_matrix - D %*% gen_comp %*% D rownames(rg_matrix) - colnames(rg_matrix) - c(Yield, Protein) print(rg_matrix)这一步其实就是在做公式里的标准化换算。很多新手搞不清楚为什么自己算出来的“遗传相关”和软件输出的不一样多半是在这里漏了标准化步骤。3.4 模型收敛诊断看到“成功”也别急着收工sommer跑完默认会报告是否收敛但我不会只看这个标志还会额外做两个检查。第一重新用不同的初值跑一遍看方差组分落到差不多的区间没有。如果两次运行的结果差异很大说明模型可能卡在局部最优解附近此时要考虑简化模型或者增修数据。第二检查对数似然值变化曲线。sommer里可以这样看迭代过程mod$monitor # 查看迭代历史如果logLik到后期还在明显爬升说明迭代没进入稳定区。如果已经变成一条平缓直线说明模型稳定了。不要相信“它有默认最大迭代次数所以肯定没问题”这种话混合模型包报“converged”只是一个数值判断不代表生物学上可靠。4. 遗传相关结果解读从数字到选择决策4.1 正相关、负相关与接近零相关的育种含义拿到遗传相关矩阵之后第一步不是写报告而是把每个关系放到育种背景下想清楚。正相关的情况好比“顺手牵羊”比如在这套DH群体里生育期和株高往往是正向遗传相关意味着选择早熟材料时株高也会倾向变矮——如果两个方向都是需要的那这就是一次选择能带来的双倍收益。负相关则意味着权衡。产量和蛋白质的例子就是典型。遇到负相关我的处理方式一般是先把相关强度确认清楚再看两性状的经济权重。如果都是核心目标通常走指数选择而不是单性状选择如果其中一个是次要性状那就定一个下限阈值在主性状上保持选择压力同时剔除掉低于阈值的个体。接近零相关反而常常是好消息。它说明两个性状基本可以独立选择不用互相妥协。比如某种材料的抗倒伏能力与产量遗传相关接近零那我们在选高产材料时可以放心淘汰倒伏敏感株不用担心把产量潜力也选没了。4.2 标准误别被一个点估计骗了遗传相关的标准误是很多人不重视、但恰恰是最要命的部分。点估计等于0.35如果标准误0.28那这个相关的可信区间无限逼近零基本等于没有如果标准误0.08才算得上一个可靠的弱相关。标准误的来源主要有三类。REML框架下可以用Delta法从方差组分的信息矩阵近似推导也可以用Bootstrap对材料重抽样重估贝叶斯路线则直接把后验样本的标准差当作标准误。sommer里想快速拿到标准误的一个粗暴办法是多次运行模型用估计结果的波动范围作为参考——虽然不严谨但很直观。影响标准误大小的因素里样本量是最主要的。120份材料的DH群体估出的遗传相关标准误通常就落在0.15到0.3之间只能说判断方向要稳定到0.1以下往往需要几百份材料加上合理的重复设置。另外遗传力越低遗传相关标准误越大两个性状之一遗传力接近零时标准误会大得离谱此时最好直接放弃报告该相关值。4.3 遗传相关如何落到田间选择跑遗传相关最终是为了用于选择方案。最简单的用法是间接选择用前面提到过的公式判断效率高不高。比如目标性状是根系抗逆性鉴定成本很高而另一个地上部形态性状容易测遗传相关达到0.6以上就可以在早代大群体里用地上部性状做初筛。更系统的用法是构造选择指数。选择指数本质上是给每个目标性状赋一个经济权重然后从遗传参数矩阵出发求解最优权重向量[ I b_1 x_1 b_2 x_2 \cdots b_n x_n ]这里的权重向量包含遗传方差、遗传协方差、经济权重信息。实际运行项目时我会先把多性状模型的G矩阵导出来作为指数计算的输入。这样一来从数据到田间决策就走通了完整链路。需要特别注意遗传相关有环境特异性在A环境算出的负相关到B环境可能变正所以指数建立在多年多点联合分析的基础上才靠谱单环境结果只用来做内部讨论。5. 实操中的四个大坑附完整复盘过程5.1 模型不收敛先检查量纲再谈迭代次数第一次跑多性状模型时我就栽在量纲差异上。产量按kg/ha记数值动辄几千蛋白质含量按百分比记数值只有十几。两者方差差了三个数量级方差协方差矩阵很容易变成病态矩阵模型反复警告不收敛。排查过程是这样的先看两个性状的描述统计发现方差差异巨大然后尝试把性状按单位标准差标准化重建数据后再跑收敛就顺畅多了。标准化不会改变遗传相关的值因为遗传相关本来就是标准化后的协方差也就是纯量纲比值。所以我的建议是多性状模型里尽量让各性状方差在一个数量级内用scale()处理一下既加速收敛也减少数值预警。输出遗传相关时标准化的影响会被约掉不用担心结果失真。5.2 遗传相关超出[-1,1]方差组分一定有问题遗传相关的定义决定它理论上必须在-1和1之间但我见过不少初学者跑出1.3、-1.7这类离谱值。出现这种情况一定不是公式错了而是方差组分估计出了问题。最常见的原因是模型没有真正收敛遗传协方差和遗传方差的比例在迭代末期还没稳定次常见的原因是样本量太小或亲缘关系矩阵奇异导致方差组分被估到边界。复盘我遇到的一次异常估计问题是系谱中有一个材料被错误登记为许多材料的亲本导致A矩阵里出现一批异常高的亲缘关系把遗传协方差撑大了。解决办法是重新核对系谱记录把可疑个体标记出来剔除后重跑。这类问题排查起来很耗时但结果验证值得做我把剔除前后两次的遗传相关画在同一张图上发现剔除异常系谱后估计值整体稳定了很多。5.3 固定效应设置不当会让遗传相关“虚胖”混合模型里固定效应怎么放直接影响方差组分在大类之间的分配。环境效应如果不放进固定效应田间环境差异会被吸收进随机效应和残差遗传方差的估计就会被抬高或扭曲遗传相关的点估计自然不可信。我现在的做法环境作为固定效应放进去重复、区组如果设计合理也作为固定效应或随机效应放进去有多年份重复的年份×环境交互作用通常按随机效应处理。一个实用的检查方法是把同一数据分别用“只有截距”和“包含环境固定效应”两个模型跑一遍看遗传相关变化有多大。如果变化剧烈说明固定效应set的好不好直接决定了故事的结论这种情况下务必谨慎。5.4 亲缘关系矩阵的坑A矩阵和G矩阵各有脾气A矩阵的坑主要在系谱质量。田间记录、人工授粉记录、编号差错率通常有5%左右单个错误系谱就足以让A矩阵某块数值偏离真实值。修剪系谱时min.common参数不要设太高否则祖先信息保留太少矩阵会变得稀疏。G矩阵的坑更隐蔽。基于分子标记计算亲缘关系时标记缺失率、等位基因频率、连锁不平衡都会影响G矩阵质量。用G矩阵做的遗传相关在小样本里常常和高通量标记的LD结构高度相关导致估计值偏乐观。我的交叉验证经验是同一批材料同时用A矩阵和G矩阵各跑一次模型如果两个来源的遗传相关方向一致结论就相对放心如果方向相反先别怀疑生物学先回头检查标记质控流程和系谱记录。6. 关于遗传相关我在实践里的几条体悟数据代码演示做到这一步核心流程已经完整了数据准备、多性状混合模型、遗传相关计算、结果解读、异常排查。最后分享几个我自己沉淀下来的习惯。第一永远不要只报一个遗传相关点估计。至少要报告标准误最好给出大致可信区间。没有不确定度量的遗传相关在审稿和内部决策中都缺乏说服力。第二遗传相关需要交叉验证。随机把数据分成两半分别估计遗传相关如果两次结果的符号和强度都接近这个相关才算初步站得住。尤其是那种“终于跑出一个理想相关”的时候更要做这个验证——数据挖掘里最容易骗人的就是好消息。第三把遗传相关和BLUP育种值结合起来看。遗传相关是一群材料的统计趋势BLUP是个体水平的预测值。两者结合才能回答“哪些单株既高产又高蛋白”这种具体问题。光有相关没有个体预测选择无处下手光有个体预测没有相关结构又难以理解为什么选了这个丢了那个。如果你正准备跑第一批多性状数据我的建议是不要在软件参数上花太多时间先把数据整理和固定效应设置想清楚。数据干净、模型合理遗传相关自然会出现一个可信的答案数据一塌糊涂再高级的算法也无非是把噪声算得更加精致。希望这篇演示能让你少走几段我走过的弯路。