利用蒙特卡洛抽样近似多性状基因组预测模型中的预测误差方差与可靠性

《Theoretical and Applied Genetics》：Approximating prediction error variances and reliabilities in a multiple-trait genomic prediction model using Monte Carlo sampling

　　

核心信息：研究人员评估了四种现有基于蒙特卡洛（Monte Carlo, MC）的方法，用于近似多性状基因组预测模型中的预测误差方差（prediction error variance, PEV）和可靠性。所有四种方法均获得了与精确值一致的近似结果。

摘要：基因组预测模型，如基因组最佳线性无偏预测（genomic best linear unbiased prediction, GBLUP），利用基因组数据来提高估计遗传值的准确性。随着基因型和性状数量的增加，由于需要反转混合模型方程组（mixed model equations, MMEs）的系数矩阵，而该矩阵的维度随个体数和性状数直接增加，预测误差方差（PEVs）和可靠性的精确计算在计算上变得不可行。本研究的目的是评估蒙特卡洛（MC）抽样方法在与杂种育种相关的多性状GBLUP框架中近似PEVs和可靠性的适用性。MC方法通过从假定分布中重复抽样遗传值来近似PEVs，从而避免直接的矩阵求逆。研究人员采用四种先前发表的公式应用MC方法近似PEVs和可靠性。所有公式均产生了一致的PEVs和可靠性估计，其收敛速率取决于公式、可靠性水平以及MC样本量。

### 研究背景与问题

在基因组预测中，预测误差方差（prediction error variance, PEV）和可靠性是评估估计遗传值（estimated genetic value, EGV）精度的关键指标。PEV定义为真实遗传值与其估计值之差的方差，可靠性则定义为真实遗传值与估计值之间相关性的平方，取值范围为0到1。在混合模型框架下，这些量的精确计算需要反转混合模型方程组（mixed model equations, MMEs）的系数矩阵，而该矩阵的维度随个体数和性状数的增加而增大，导致计算量急剧上升。对于多性状基因组最佳线性无偏预测（multiple-trait genomic best linear unbiased prediction, MT-GBLUP）模型，当性状数和基因型个体数较多时，直接矩阵求逆在计算上变得不可行。

已有多种基于蒙特卡洛（Monte Carlo, MC）抽样的方法被提出用于近似PEV和可靠性，这些方法通过重复模拟遗传值并求解MMEs来避免直接矩阵求逆。此类方法最初主要用于动物育种模型，其在植物育种特别是杂种育种框架下的适用性尚需系统评估。植物育种与动物育种在群体结构、交配设计以及非加性遗传效应的重要性等方面存在差异，因此有必要在植物育种背景下验证这些方法的有效性。

该研究发表于《Theoretical and Applied Genetics》，其目的并非提出新理论，而是评估已有MC近似方法在MT-GBLUP模型中的表现。研究人员采用一个具有三个随机效应组的MT-GBLUP模型，即两个杂种亲本群体的一般配合力（general combining ability, GCA）效应与杂交种的特殊配合力（specific combining ability, SCA）效应，将四种已有的PEV近似公式应用于杂交育种场景，并与精确值进行比较。

### 关键技术与方法

研究人员基于García-Cortés等人和Hickey等人提出的MC算法，通过三个步骤近似PEV：模拟遗传值、残差和表型；从模拟数据求解MMEs；利用模拟真实值与预测值计算PEV近似值。遗传值的模拟可采用基于标记矩阵的策略（Marker-based MC）或基于基因组关系矩阵Cholesky分解的策略（Cholesky MC）。MMEs使用预条件共轭梯度（preconditioned conjugate gradient, PCG）迭代算法求解。研究比较了四种近似公式：PEV1、PEV2、PEV3和NF2，其中PEV3是PEV1与PEV2按经验方差加权的组合。模拟数据使用AlphaSimR软件生成，模拟玉米基因组，包括10条染色体、共50,000个SNP标记。两个杂种亲本群体各含1000个双单倍体（double haploid, DH）系，随机交配产生10,000个F1杂交种。表型记录仅针对杂交种，包含两个性状在两个环境中的表达（四个性状状态），每个杂交种每个性状有10次重复观察，10%的观察随机缺失。遗传协方差矩阵设定为：性状A窄感遗传力约0.28、环境间相关0.80，性状B窄感遗传力约0.17、环境间相关0.70，两性状间轻微负相关，SCA方差设为0.10。通过反转MME系数矩阵获得精确PEV和可靠性，并将MC近似值与精确值进行比较。

### 研究结果

#### 近似的一致性

研究发现，四种MC近似方法均产生了与精确PEV和可靠性一致的近似结果。在500个MC样本下，除PEV2外，各方法近似可靠性值与精确值之间的相关均超过0.99，均方根误差（root mean squared error, RMSE）最大为0.047，最大绝对差（MAX）介于0.055和0.251之间，回归分析的截距接近0、斜率接近1。这表明MC方法能够较为准确地近似MT-GBLUP模型中的PEV和可靠性。

#### 收敛速率与可靠性水平的关系

不同方法的收敛速率存在差异，这与各方法的渐近抽样方差有关。PEV1在可靠性较低时近似效果较好，PEV2在可靠性较高时表现更好，而PEV3和NF2在所有可靠性水平下均具有相对较低的抽样方差，因此稳定性更高。随着MC样本量增加，各方法的RMSE均下降，但下降趋势因方法和组别而异：对于两个亲本群体，PEV2、PEV3和NF2的RMSE趋势相似，PEV1的RMSE始终较高；对于杂交种群体，PEV2的RMSE则高于其他方法。研究还发现，当按群体-性状组合分别分析时，由于各组合可靠性范围较窄，相关系数、截距和斜率等指标可能低估或高估方法的真实表现，因此应同时报告RMSE和MAX。

#### 两种模拟策略的比较

为验证SCA互作标记矩阵的近似构造不影响结果，研究人员比较了Marker-based MC与Cholesky MC两种策略。结果显示，两种策略在不同群体和方法下产生的RMSE几乎相同，说明遗传值模拟策略的选择对近似精度没有实质影响。计算复杂度方面，Cholesky MC在标记数m满足n_G<2m时计算量较小，而Marker-based MC在标记数较少时更具优势。

#### 计算时间与可扩展性

计算时间测试表明，精确方法的时间随数据规模增大呈立方级增长，而PCG迭代算法近似呈二次方增长。对于最大的数据集（MME系数矩阵维度160,000），精确求逆需要172.70分钟，而单次PCG求解仅需1.92分钟。由于MC样本相互独立，可将抽样与求解任务分配到多个计算节点并行执行。以50个节点（每节点10核）和500个MC样本为例，预计总耗时约35分钟，比精确方法快四倍以上。此外，PCG算法支持“基于数据的迭代”（iteration on data），无需在内存中形成完整的MME系数矩阵，从而显著降低对计算机内存的需求。

### 讨论与结论

讨论部分指出，MC近似方法的主要优势在于通用性强，可较容易地适配多种模型，并且能够避免直接矩阵求逆和降低内存需求；其主要缺点是需要多次求解MMEs，但可通过并行计算有效缓解。研究也指出了一些局限性：在模拟中遗传和残差（协）方差被假定为已知且正确设定，未考虑方差组分估计不确定性对PEV和可靠性近似的影响。实际应用中，方差组分的估计质量会直接影响近似结果的准确性。另外，本研究的最大数据规模为160,000，外推到更大数据集时需谨慎，但鉴于矩阵求逆的时间复杂度为立方级，而PCG近似为平方级，迭代求解方法的相对优势预计会随数据规模增大而更加明显。研究还讨论了其他不依赖MC抽样的可靠性近似方法，如单步GBLUP中的近似方法。

研究结论可概括为：研究人员调查了MT-GBLUP模型中基于MC抽样的PEV和可靠性近似方法，结果表明源于动物育种的方法同样适用于植物育种。这些方法能够产生与精确PEV和可靠性一致的近似值，该结论与先前研究一致。然而，需多次求解MMEs所带来的计算负担可能限制这些方法在超大规模数据集中的应用。