面向超大规模数据集和复杂单步基因组模型的(协)方差组分估计新方法：蒙特卡洛单步基因组REML(MC-ss-GREML)的开发与应用

《Genetics Selection Evolution》：Estimation of (co)variance components for very large datasets and complex single-step genomic models

本研究针对大规模数据集或复杂模型下，由于计算限制而无法使用所有数据和信息进行无偏方差组分估计的问题，开发了蒙特卡洛单步基因组REML(MC-ss-GREML)方法。该方法通过模拟育种值并求解混合模型方程来近似计算涉及预测误差方差的迹，成功将蒙特卡洛REML扩展到包含基因组信息的情境。研究结果表明，MC-ss-GREML在保证估计准确性的同时，显著降低了计算时间和内存消耗（分别仅为传统精确方法的14%和1%），并能处理包含数百万个体谱系和数十万基因型数据的大规模数据集，为大规模基因组选择中的遗传参数估计提供了高效可行的解决方案。

　　

在动物育种和遗传改良领域，准确估计遗传参数的方差和协方差（即方差组分）是进行可靠遗传评估和做出正确选种决策的基石。传统的限制性最大似然法（REML）是估计这些参数的黄金标准。然而，随着基因组学时代的到来，数据量呈现爆炸式增长，现代育种项目往往涉及数百万个体的谱系和数十万个体的基因型数据。将如此海量的基因组信息整合进统计模型（例如单步基因组最佳线性无偏预测，即ssGBLUP模型）后，模型变得异常复杂，对计算资源提出了近乎苛刻的要求。传统的REML方法需要反复求逆一个巨大的混合模型方程系数矩阵，这导致计算时间漫长且内存消耗惊人，使得针对大规模数据的精确计算变得不切实际。

由于计算上的瓶颈，研究人员和育种公司在实际操作中常常被迫采取妥协方案：要么只使用数据的一个子集进行分析，要么简化模型（例如忽略基因组信息或剔除某些性状）。这种“将就”的做法虽然解了燃眉之急，但其后果是严重的——可能导致估计的方差组分产生偏差，进而使得估计的育种值、遗传进展评估以及综合选择指数的计算都不准确，最终可能导向错误的选种决策，影响育种效率。

为了突破这一计算困境，由Matias Bermann等人组成的研究团队在《Genetics Selection Evolution》杂志上发表了他们的最新研究成果。他们成功开发了一种名为“蒙特卡洛单步基因组REML”（Monte Carlo single-step genomic REML, MC-ss-GREML）的新方法，旨在以可承受的计算成本，为超大规模数据集和复杂的单步基因组模型提供准确的方差组分估计。

为了开展这项研究，研究人员主要应用了几个关键技术方法。核心是蒙特卡洛REML（MC-REML）框架，其基本思想是通过随机模拟来近似计算REML似然函数中难以直接求解的迹（trace）运算。本研究的关键创新在于将这一框架扩展至单步基因组模型（ssGBLUP）。他们改进了Pimentel等人提出的育种值模拟算法，使其能够高效地生成符合ssGBLUP模型假设的育种值向量，而无需显式构造或求逆庞大的基因组关系矩阵G或其扩展形式H矩阵。在优化策略上，研究结合了期望最大化（EM）算法和平均信息（AI）REML算法来迭代更新方差组分估计值。对于超大规模数据集，还采用了算法验证与年轻动物（APY）近似法来稀疏化基因组关系矩阵的逆，进一步提升计算效率。测试数据来自美国安格斯牛协会，包括一个用于验证方法准确性的小型数据集（约10万动物谱系，其中3.3万有记录，1万有基因型）和一个用于评估计算性能的大型数据集（约740万动物谱系，其中33万有基因型）。

研究结果

MC-ss-GREML估计精度验证

研究人员首先在一个三性状（出生重、断奶重、断奶后增重）的安格斯牛生长模型上测试了MC-ss-GREML的准确性。该模型包含直接遗传效应、母体遗传效应、母体永久环境效应等随机效应，共计需要估计14个方差协方差参数。他们将MC-ss-GREML的估计结果与使用精确迹计算的传统ss-GREML（exact ss-GREML）结果进行比较。结果表明，两种方法得到的参数估计值非常接近，统计检验显示没有显著差异。这证明了MC-ss-GREML在估计精度上可以达到与传统精确方法相当的水平。

计算性能显著提升

在计算效率方面，MC-ss-GREML展现出了巨大优势。对于上述三性状复杂模型，MC-ss-GREML仅消耗了传统ss-GREML所需计算时间的14%，同时内存占用仅为后者的1%。这种显著的性能提升主要归功于MC-ss-GREML避免了存储庞大的混合模型方程系数矩阵和进行耗时的矩阵求逆运算，转而采用迭代求解和数据模拟的策略。

大规模数据集应用实例

为了进一步测试MC-ss-GREML处理超大规模数据的潜力，研究团队将其应用于一个仅针对出生重性状的单性状模型，但数据量极大，包含约740万动物的谱系，其中近578万有记录，33.1万有基因型。在此规模下，MC-ss-GREML成功实现了收敛，整个过程经历了11轮迭代，总计算时间为121小时，峰值内存使用量为53 GB。这个结果证实了MC-ss-GREML有能力处理现实世界中育种项目所产生的大规模数据，并且在当前的高性能计算服务器上是可行的。

收敛性管理

研究还注意到，由于蒙特卡洛模拟引入的随机误差，基于参数估计值轮间变化（Δθ）的经典收敛准则可能难以达到。为此，他们提出结合评估近似的限制性对数似然值在连续几轮迭代中的变异系数（CV(log L)）作为辅助收敛标准。当Δθ小于一个严格阈值（如10^-12）或CV(log L)小于一个相对宽松的阈值（如10^-4）时，即可判定收敛。这种双重标准策略有效地解决了蒙特卡洛方法带来的收敛波动问题。

研究结论与意义

本研究成功开发并验证了MC-ss-GREML这一新型方差组分估计方法。该方法的核心贡献在于，它首次将蒙特卡洛REML技术有效地扩展至包含基因组信息的单步模型框架内，克服了传统方法在处理超大规模基因组数据时面临的计算瓶颈。

研究表明，MC-ss-GREML能够在保证估计准确性的前提下，极大地降低计算时间和内存需求。这使得利用全部可用数据（包括完整的谱系和基因组信息）进行复杂模型的方差组分估计成为可能，从而避免了因数据子集化或模型简化而可能导致的估计偏差。对于正处于高强度基因组选择下的种群而言，使用包含所有信息的完整模型进行参数估计尤为重要，因为选择压力可能会改变群体的遗传结构，忽略基因组信息可能导致有偏的估计。

此外，由于MC-ss-GREML采用迭代-on-data（IOD）技术，无需将整个模型矩阵载入内存，因此非常适合集成到常规的遗传评估流程中。该方法为动物育种行业应对日益增长的数据挑战提供了强有力的工具，有助于提高遗传评估的准确性和可靠性，最终推动更高效、更精准的育种实践。尽管确定最优的蒙特卡洛模拟次数等因素需要进一步研究以优化计算效率，但MC-ss-GREML无疑为大规模基因组选择时代的遗传参数估计开辟了一条新的可行路径。