机器学习做基因组选配文献笔记

机器学习做基因组选配文献笔记

读了一篇用集成5种机器学习算法做基因组选配的文献，后代基因型通过亲本基因型模拟得到，每队亲本模拟2000个后代基因型，通过预测这2000个后代基因型的表型预测值的分布来判断是否选配。

文献1

摘要

核心背景与痛点

• 传统局限： 基因组选择（GS）以往大多只关注“基因型到表型”的预测模型，往往忽略了决定后代选择潜力的最佳杂交设计（Optimal Cross Design）。

EMLGP 框架的核心机制

• 方法结合： EMLGP（集成机器学习基因组预测）首先利用机器学习算法在种质群体中训练精准的表型预测模型，随后结合基因组模拟技术，预测育种群体中的最佳杂交组合。

突出性能与优势

• 准确率大幅提升： 在 14 个大豆性状的预测训练中表现优异，最高预测准确率（相关系数）达 0.92，比传统 GBLUP 方法最高提升了 35.85%。
• 高鲁棒性（稳定性）： 模拟研究证实，在中小样本量（300–5000）、中低遗传力（0.4–0.6）以及复杂遗传结构（100 个因果位点）下，该框架依然表现稳健。
• 多物种验证： 在水稻、玉米、棉花、高粱和柳枝稷的真实数据验证中，其表现持续优于 GBLUP 和深度学习方法。

育种潜力与结论

• 超亲分离潜力： 在分析的 14 个大豆性状中，有 13 个在后代中表现出超亲分离潜力。其中，华北地区大豆种子的亚麻酸含量重组潜力最高，超过亲本最大值 16.89%。
• 总结： EMLGP 优化了亲本选择和表型预测，为高效、智能驱动的作物育种提供了一个强有力的框架。

引言

GBLUP 和 Bayes 算法

略

机器学习

机器学习方法可以捕捉复杂的基因模式，目前常用的机器学习方法有 SVR 和 RF 等。但是单个机器学习有2个问题：

1. 预测效果显著受到数据集的影响
2. 超参数调参让人烦扰

这两个问题限制了机器学习的广泛应用。

因此，集成学习通过使用多个模型来提高模型鲁棒性和预测效果，已经变成了现在的主流策略。

深度学习

深度学习可以捕捉非线性关系和基因组上长距离的互作效应。一开始使用深度学习是 DeepGS 包，之后是 DNNGP，然后是 SoyDNGP 采用类似图片的基因型输入格式。GPformer 方法利用 Transformer 框架捕捉 DNA 序列中的长距离相互作用和依赖关系。

方法

EMLGP 框架

EMLGP 框架有2个程序：

1. 创建一个基因组选择模型
2. 预测后代潜力

表型拆分为 均值 + 基因型值（total genotypic value） + 残差。

基因型值 ($g_{i}$) 可以表示为关于所有基因位点集合的函数。
$$
\mathbf{g}i = f(\mathbf{x}{i1}, \mathbf{x}{i2}, \cdots, \mathbf{x}{im})
$$
其中 $x_{i1}$ 等是每个位点的基因型，取值为 0 1 2。

阶段1：在种质资源群体中训练GS模型

我感觉他这里前两步都是为了降维。

1. 降维：基于连锁不平衡（LD）筛选标签SNP（PLINK软件 --indep 和 --show-tags all）。
2. 特征选择：使用 Boruta 算法（基于随机森林）进一步筛选与目标性状相关的关键标记。
3. 多模型超参数调优（训练过程）：包括 Lasso、Ridge（ridge regression）、SVR、RF、XGBoost （使用 scikit-learn 软件 ）。研究人员使用 scikit‑learn 库中的 GridSearchCV 工具，通过 网格搜索 和 交叉验证，分别对 Lasso、Ridge、SVR、RF、XGBoost 这 5 种模型进行超参数调优，并以 RMSE（均方根误差） 作为评估指标（就是 MSE ，只是开方值），最终选出每个模型中表现最好的那组超参数。
4. 堆叠泛化（预测过程）：使用一个 线性回归模型 作为“元学习器”，把这些基础模型的预测值作为输入，最终输出一个综合的、更准确的预测结果。这个组合过程由 scikit‑learn 的 StackingRegressor 完成，所有参数保持默认。

堆叠泛化的具体解释

简单地说，就是对所有单个模型的预测值做一个加权回归，得到一个最终的值。

加权平均相比于简单平均的优势将下表。

简单平均	加权平均（线性回归）
所有模型权重相等，即使某个模型很差也会拖累结果	线性回归自动学习最优权重，表现差的模型权重低甚至为 0
无法利用模型之间的互补性	可以捕捉不同模型在不同样本上的相对优势
对异常预测不鲁棒	线性回归带截距，可以自动校准偏差

堆叠泛化的具体做法（结合 EMLGP 流程）

第一步：准备数据

• 训练集：基因型矩阵 $X_{\text{train}}$ 和对应表型向量 $y_{\text{train}}$。
• 测试集：基因型矩阵 $X_{\text{test}}$（无表型）。

第二步：训练基础模型（已经调优好的最优版本）

分别训练：

• Lasso（最佳 $\alpha$）
• Ridge（最佳 $\alpha$）
• SVR（最佳 $C, \epsilon, \gamma$）
• RF（最佳树数量、深度等）
• XGBoost（最佳学习率、树深度等）

每个模型都在 相同的训练集 上拟合。

第三步：生成第一层预测（用于训练元模型）

• 对训练集，需要得到每个基础模型的 交叉验证预测（否则会导致元模型过拟合）。
  常见做法：5 折交叉验证，每折中，模型在 4 折上训练，对剩余 1 折预测。最后把 5 折的预测拼接起来，得到 训练集的次级特征 $Z_{\text{train}}$（形状为 $n_{\text{samples}} \times 5$）。
• 对测试集，每个基础模型都在 全部训练集 上重新训练，然后预测测试集，得到 $Z_{\text{test}}$（形状为 $n_{\text{test}} \times 5$）。

第四步：训练元学习器（线性回归）

• $$
  y_{\text{final}} = \beta_0 + \beta_1 \cdot \hat{y}{\text{Lasso}} + \beta_2 \cdot \hat{y}{\text{Ridge}} + \dots + \beta_5 \cdot \hat{y}_{\text{XGBoost}}
  $$

• 优化目标是最小化均方误差（MSE）。

第五步：在测试集上预测

• 使用训练好的线性回归模型，将 $Z_{\text{test}}$ 作为输入，得到最终的集成预测值。

第六步：scikit-learn 的 StackingRegressor 自动完成这些

from sklearn.ensemble import StackingRegressor
from sklearn.linear_model import LinearRegression
from sklearn.svm import SVR
from sklearn.ensemble import RandomForestRegressor
from xgboost import XGBRegressor
from sklearn.linear_model import Lasso, Ridge

# 定义基础模型（已经调好参的）
base_models = [
    ('lasso', Lasso(alpha=0.01)),
    ('ridge', Ridge(alpha=1.0)),
    ('svr', SVR(C=1.0)),
    ('rf', RandomForestRegressor(n_estimators=100)),
    ('xgb', XGBRegressor(learning_rate=0.1))
]

# 定义元模型（线性回归）
meta_model = LinearRegression()

# 创建 StackingRegressor，使用默认参数
stack = StackingRegressor(
    estimators=base_models,
    final_estimator=meta_model,
    cv=5  # 默认就是 5 折交叉验证
)

# 训练
stack.fit(X_train, y_train)

# 预测
y_pred = stack.predict(X_test)

参数 cv=5 就是内部自动做交叉验证生成第一层训练集预测，避免过拟合。final_estimator 默认是 LinearRegression()，与文中一致。

阶段2：在育种群体中预测最优杂交组合

1. 子代基因型模拟：对每一对亲本组合，模拟2000个纯合后代（从 F1 开始，自交多代直到所有位点纯合的过程）。

   • 在连锁模型中，他们用一个泊松过程来模拟减数分裂时染色体上发生交叉互换的位置和次数。染色体越长（单位 Morgan），平均交换次数越多。并且假设交换事件相互独立（没有干涉）。（如果不是纯合子，其实是需要先定相的）
   • 作为对比，他们还采用了一个简化模型——假设所有 SNP 位点之间没有连锁关系，完全独立地分离。这相当于把所有基因座看作分布在不同的染色体上（或者认为它们之间自由组合）。这个模型计算快，但生物学上不够真实。

   估计大豆亲本基因型应该都是纯合子，这样 F1 才是完全确定的基因型。从 gitee 上的示例数据来看确实全是纯合子。

   这个连锁模型没太看懂，AI 解释如下

   它把一个亲本的染色体传给后代的遗传过程，用计算模型拆解为三个主要步骤：定次数、定位置、定基因。

   1. 第一步：定次数——交叉互换的次数
      • 原理：一条染色体上发生几次交叉互换，是随机的，但总体趋势符合 泊松分布 (Poisson Distribution)。遗传长度（单位：Morgan）就是平均交叉互换的次数（参数 $\lambda$）。
      • 一个例子：假设一条染色体遗传长度为 1 Morgan，则它平均发生1次交换。但实际模拟时，有 36.8% 的概率发生0次，36.8% 的概率发生1次，18.4% 的概率发生2次……。
      • 一句话总结：遗传图距越长，染色体越长，平均发生的交叉互换次数就越多。
   2. 第二步：定位置——交叉互换发生的位置
      • 模拟方法：一旦用泊松分布确定了某条染色体的具体交叉次数（比如 k 次），这 k 个交叉点的位置就会独立且均匀地分布在整个染色体上。
      • 关键词：整个过程是随机的，并且假设 无交叉干涉 (No Crossover Interference)，意味着一个交叉的发生不会影响另一个交叉。
   3. 第三步：定基因——生成后代的配子和基因型
      • 按点切分：在染色体上，根据第二步确定的交叉点位置，将染色体切成一段一段的区间。
      • 交替拼接：将这些区间按照“来自父本A”、“来自父本B”、“来自父本A”……的顺序交替排列，就组装成了一条全新的、包含两个亲本遗传信息的重组染色体（配子）。（反过来，另一个配子就是父本B，父本A，父本B的顺序）
      • 获得纯合后代：将两个亲本各自的模拟配子结合，得到F1代。但育种需要的是性状稳定的纯合系（RIL）。所以，模型会在这个F1基础上，模拟连续多代的自交过程，直到所有位点纯合。你看到的“2000个纯合后代”，就是对这个过程的程序化重复。

   具体做法：

   1. 染色体遗传长度（单位：Morgan，1 Morgan = 100 cM 。人类基因组：平均约 1 cM ≈ 1 Mb）得到泊松分布的唯一参数 $\lambda$ ，然后基于泊松分布抽一次样得到具体交叉次数 K 。
   2. 计算物理位置：根据第一步确定的次数 k，在 [0, 染色体总长度] 的区间内，生成 k 个独立的、服从均匀分布的随机数。

   python 示例代码如下

   import numpy as np
   
   # 设置模拟参数
   chr_length_morgan = 1.0      # 染色体遗传长度: 1 Morgan
   chr_physical_length = 100.0  # 染色体物理长度（可设为任意单位）
   
   # 1. 使用泊松分布随机确定交叉互换次数
   n_crossovers = np.random.poisson(chr_length_morgan)
   
   # 2. 若发生交叉，确定具体位置
   if n_crossovers > 0:
    # 在染色体物理长度范围内随机均匀分布生成交叉点
    positions = np.random.uniform(0, chr_physical_length, n_crossovers)
    positions.sort()  # 按顺序排列
   
    print(f"发生了 {n_crossovers} 次交叉")
    print(f"交叉点位: {positions}")
   else:
    print("未发生任何交叉")

   这个连锁模型动物上怎么做？

   第一步先定相，后面一样

   为什么要模拟 完全纯合 的后代？（这是核心问题）

   育种中，从杂交后代选育新品种时，最终获得的稳定品系都是 纯合系（自交多代后所有位点纯合）。大豆是自花授粉作物，育种流程通常是：
   亲本A × 亲本B → F1（杂合） → 自交多代（F2，F3，...，F6） → 获得纯合的重组自交系（RIL）
   这些 RIL 一旦纯合，就可以稳定遗传，成为新品种或育种材料

2. 子代表型预测：使用训练好的GS模型预测每个后代的表型。

3. 超越潜力评估：计算后代表型的95%分位数等指标，判断是否超越亲本最大值。具体解释，取 95% 分位数，这个组合中，意味着表现最好的 5% 后代的表型下限，也就是后代极端优秀部分的潜力（如果是负向选择，那就是 5% 分位数）。AI 的解释如下

   育种实践中，育种家不会把所有后代都保留，而是只选择 前 5% 或前 10% 的优良单株。因此，95% 分位数正是对应“保留最优秀 5% 后代”这一选择强度下的预期表型值。数值越高，说明该组合越有可能通过高强度选择获得突破性品种。

EMLGP 框架

EMLGP 软件 gitee 网址

https://gitee.com/njau-sri/EMLGP

模型比较

和下面的模型进行比较预测准确性。采用5折交叉验证获取准确性。

• GBLUP
• DNNGP
• scikit-learn 软件中的机器学习算法，包括 Ridge, Lasso, EN, SVR, RF, XGBoost, and LightGBM 。
  • 这里 Ridge, Lasso, and EN 都是正则化的线性回归方法，只是正则化的方式不同。
  • RF, XGBoost, and LightGBM 这三个是决策树算法。

结果

没看，让 AI 汇总的。我看计算效率上相比于 GBLUP 还是慢很多啊。

一、EMLGP 在 14 个大豆性状上的预测精度显著优于传统方法

• 在 CSGP1（821份栽培大豆，9个性状）和 CSGP2（750份野生+栽培大豆，5个性状）两个群体中，EMLGP 的预测精度（Pearson相关系数 PCC）在多数性状上达到 0.75–0.92。
• 开花期（DTF） 和 成熟期（DTM） 的 PCC 最高，分别达到 0.92 和 0.90。
• 相比经典方法 GBLUP，EMLGP 的预测精度最大提升 35.85%（例如大豆亚油酸含量 LA：PCC 从 0.53 提升至 0.72）。
• 同时，EMLGP 在所有测试的机器学习方法（Lasso、Ridge、SVR、RF、XGBoost）及深度学习方法（DNNGP）中均表现最佳。

二、在不同作物（水稻、玉米、棉花、高粱、柳枝稷）中稳健领先

• 棉花纤维长度/强度：EMLGP（PCC=0.83）优于 GBLUP（0.78–0.79）和 DNNGP（0.69–0.70）。
• 玉米株高（HT）：EMLGP 预测精度（0.85）相比 GBLUP（0.45）提升 88.9%。
• 柳枝稷：株高稳定性（ST）PCC=0.88，开花期（FT）PCC=0.92。
• 高粱：籽粒水分（GM）和株高（HT）的 PCC 分别为 0.78 和 0.81，显著高于 GBLUP/DNNGP（0.55–0.57）。
• 结论：EMLGP 具有跨物种的泛化能力和稳定性。

三、在小样本、低遗传力、复杂遗传结构下依然稳健

通过模拟不同样本量（300–5000）、遗传力（0.4–0.8）和 QTL 数量（10、100）：

• 当 QTL 数量从 10 增加到 100（更复杂遗传结构）时，GBLUP 和 EMLGP 的预测精度均上升，但 EMLGP 始终更高。
• 在 样本量仅 300、遗传力 0.4、100 个 QTL 的严苛条件下，EMLGP 仍保持较好的预测能力。

四、成功预测大豆 13 个性状的超越潜力，筛选出优良亲本

• 在 14 个大豆性状中，有 13 个性状 的后代表现出超越双亲的潜力（即存在超亲个体）。
• 华北地区（NC） 的 亚麻酸含量（LNA） 重组潜力最大，比双亲最大值高出 16.89%。
• 基于三个生态区（华北 NC、华中 CC、华南 SC），给出了高蛋白、高油及双高（蛋白+油）的推荐亲本材料：
  • 高蛋白潜力：N25486（NC）、N23633（CC）、N23630（SC）
  • 高油潜力：N25318（NC）、N24192（CC）、N25151（SC）
  • 双高材料：NC 区有 9 份（如 N25324、N24461 等），SC 区有 4 份（如 N25403、N27998 等）

五、计算效率（尽管较高，但可接受）

在 750 个样本、274 万个标记的数据集上：

• EMLGP 总运行时间 8.09 小时，峰值内存 2.14 GB。
• 作为对比，GBLUP 耗时 0.63 小时（内存 3.65 GB），DNNGP 耗时 3.25 小时（内存 0.26 GB）。
• 虽然 EMLGP 耗时最长，但考虑到其显著的精度提升和育种实用性，作者认为可以接受。

参考文献

1. Feng W, Tai S, Liu F, et al. An efficient machine-learning framework for genomic selection of optimal crosses in soybean germplasm population[J]. The Crop Journal, 2026.