做育种数据分析这几年，我越来越觉得遗传相关是个“安安静静决定成败”的参数。表型相关很容易算，肉眼也大概看得出来，但真正决定我们能不能用一个易测性状去替难测性状做间接选择、能不能在多个目标性状之间找平衡的，是那个看不见的遗传相关。前几天我们把一份多环境试验数据的遗传相关跑出来时，发现产量和蛋白质含量存在明显的负向遗传相关——这解释了我们前几年为什么高产品系品质总是上不去。这篇文章就从这份工作出发，完整演示一遍遗传相关的数据代码分析流程：从数据整理、多性状混合模型构建，到遗传相关矩阵的计算与解读，再把收敛失败、异常估计、标准误这些我踩过的坑一并捋一遍。内容适合正在做植物或动物育种的同行，也适合刚接触遗传参数估计、想快速跑通流程的朋友。

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跑REML，MCMCglmm做交叉验证。

2.3 遗传力和遗传相关是同一模型产出的孪生参数

跑完多性状模型之后，我们通常先把遗传力也算出来。单性状遗传力 (h^2 = \sigma^2_g / (\sigma^2_g + \sigma^2_e))，在多性状模型里对应每个性状分别取对角线元素。遗传力和遗传相关出自同一个G矩阵，只是标准化方式不同：遗传力用每个性状自身的方差做分母，遗传相关用两个性状的标准差做分母。

这个“孪生”关系提醒我们一件事：如果某个性状遗传力低，它的遗传方差估计就不稳定，分母不稳定，遗传相关的估计自然跟着波动。所以我在解读低遗传力性状的遗传相关时特别保守，不会只看点估计就说“相关很强”，一定要求看标准误。

3. 数据代码演示：从三份原始表到一张遗传相关矩阵

3.1 数据准备：表型、系谱、基因型矩阵缺一不可

做遗传相关之前，先确认手上有哪几类数据。最基础的是表型数据，至少包含：材料编号（geno）、环境编号（env）、区组或重复编号（rep）、以及至少两个性状的观测值。如果材料分散在多个环境、多个重复，数据要按长表格式整理，一行一个观测记录。

第二类是亲缘关系信息。如果没有分子标记数据，可以用系谱记录构建A矩阵；如果有全基因组标记，就用标记构建G矩阵。A矩阵和G矩阵的作用是描述材料之间的遗传关联：亲缘关系越近，遗传效应之间的协方差应该越大。这一步不是可有可无，因为混合模型里随机遗传效应的方差结构必须乘以亲缘关系矩阵，否则模型等同于假设材料间完全独立，那就失去了育种数据分析的意义。

第三类是可选的协变量信息，比如播种密度、土壤肥力等级、采收日期等。这些如果已知对性状有系统性影响，就放进固定效应处理。

数据大概长这样：

这份数据来自一个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是个体水平的预测值。两者结合，才能回答“哪些单株既高产又高蛋白”这种具体问题。光有相关没有个体预测，选择无处下手；光有个体预测没有相关结构，又难以理解为什么选了这个丢了那个。

如果你正准备跑第一批多性状数据，我的建议是不要在软件参数上花太多时间，先把数据整理和固定效应设置想清楚。数据干净、模型合理，遗传相关自然会出现一个可信的答案；数据一塌糊涂，再高级的算法也无非是把噪声算得更加精致。希望这篇演示能让你少走几段我走过的弯路。