核心种质筛选-03操作

Posted on | Post modified | In

1. 解决问题

  • 目标：给出从输入数据到核心集交付物的一套“可执行流程骨架”，让项目可以按步骤复算与复用。
  • 约束：本篇只写“怎么做”，指标含义与算法原理集中在第 05 篇，交付物验收集中在第 04 篇。

2. 主要节点流程图（线性步骤）

常用软件（参考）

流程环节 软件/工具 主要用途
变异过滤/基础质控 vcftools / PLINK 缺失率/MAF/LD 处理、样本与位点过滤
核心集筛选 Core Hunter 3 / GenoCore 多目标优化或基于指标的核心集构建
群体结构 ADMIXTURE 祖先成分推断、结构分层
PCA/矩阵计算 SNPRelate（R） 高性能 PCA 与相关计算
系统发育/树 FastTree 快速构建树并辅助解释
标识生成 qrencode 将指纹信息生成二维码等输出

2.1 输入形态 A：已有 VCF/PLINK/矩阵

1) 数据准入检查（缺失率/MAF/样本异常）
2) 变异与样本过滤（形成“分析用矩阵”）
3) （可选）LD pruning（用于结构推断/降维）
4) 群体结构理解（PCA/聚类/树）
5) 核心集优化（目标：代表性/多样性；可加结构约束）
6) 梯度规模模拟（推荐）与拐点选择
7) 核心集验证（指标 + 可视化）
8) 输出交付物（清单/表/图/配置）

2.2 输入形态 B：FASTQ（需从头做变异检测）

1) FASTQ 质控与过滤
2) 比对到参考（线性/图泛基因组按物种选择）
3) 变异检测与联合分型
4) 变异过滤 → 进入 2.1 的矩阵流程

3. 贝叶斯批判性视角下的四阶段流程

思路：把每一步都视为“证据生成”，从先验（样本/数据质量）到后验（核心集是否真代表全集），逐步排除“这套核心集不靠谱”的备择假设。

阶段一：数据预处理与填充（上游准备）

目的：得到一个“对 Core Hunter 友好”的、尽量完整且已控噪的基因型矩阵。缺失值过多会直接干扰遗传距离计算与 EN/A-NE 指标。

步骤 软件/函数 关键参数（示例） 输入/输出 目的
1. 质控 vcftools / bcftools --maf 0.05 --max-missing 0.8 等 入：raw.vcf；出：filtered.vcf 剔除低质量/极低频位点，减少后续填充噪声。
2. 基因型填充 Beagle 5.4 等 gt=filtered.vcf nthreads=10 入：filtered.vcf；出：imputed.vcf.gz 在给定先验 LD 结构下，用 HMM 填补缺失，获得完整矩阵。

若群体强分层，可结合第 3.4.3 节采用“先聚类、后填充”的策略；所有阈值应在项目配置中固化。

阶段二：核心种质筛选（Core Hunter 3 核心）

目的：在可接受规模下，最大化“代表性 + 内部多样性”，并把这一目标显式写进 Core Hunter 的目标函数与约束里。

2.1 距离矩阵构建（R：SNPRelate 示例）

1
2
3
4
5
6
7
8
library(SNPRelate)

snpgdsVCF2GDS("imputed.vcf.gz", "data.gds")
genofile <- snpgdsOpen("data.gds")

## 示例：构建 GRM 后转成距离矩阵
grm <- snpgdsGRM(genofile, method = "GCTA")
dist_matrix <- 1 - grm$grm
  • 输入：imputed.vcf.gz（VCF 4.2）
  • 输出：dist_matrix（数值型对称矩阵，可直接供 Core Hunter 使用）

2.2 在 R 中调用 Core Hunter 3

1
2
3
4
5
6
7
8
9
10
11
12
13
library(corehunter)

d <- distanceData(dist_matrix)
ch_data <- coreHunterData(d)

## 示例：选取 10% 样本，兼顾 EN（内部多样性）与 GD（遗传距离）
result <- sampleCore(
  ch_data,
  size = 0.10,
  obj  = objective("EN", "GD")
)

core_ids <- result$sel  # 被选中的样本索引或 ID
  • 关键参数：
    • size：核心集比例（通常 0.05–0.20，具体由梯度模拟确定，见第 4.2 节）。
    • obj：优化目标组合；可根据业务重点（代表性 vs 多样性）调整为 EN/A-NE/NA/RA 等不同权重组合。
  • 输出：核心集样本列表（作为下游所有分析的“训练/代表集合”）。

阶段三：下游应用一——精准育种与 GS

目的：把核心集上学到的“基因型→表型”关系，外推到全群体（或更大候选集合），实现不测表型的预测与筛选。

应用步骤 软件/函数 简化模型 输入/输出 目的
1. 模型训练 R 包 rrBLUP::mixed.solve 等 (y = X\beta + Zu + e) 入：核心集基因型 + 表型；出：SNP 效应值 (\beta) 学习每个位点对目标性状的贡献。
2. 基因组预测 矩阵乘法 (GEBV = G_{\text{all}} \times \beta) 入：全群体基因型矩阵 (G_{\text{all}})；出：每个个体的 GEBV 预测值 在不测表型的前提下，筛选高产/抗病/目标性状优良个体。

这里核心集的“代表性”会直接决定 GS 模型在全集上的外推能力，这是前两阶段所有证据的业务落点。

阶段四：下游应用二——资源库管理与溯源工具

目的：把核心集与指纹信息沉淀成“可交付、可操作”的管理与合规工具。

4.1 指纹图谱与热图可视化

  • 工具：R 包 ComplexHeatmap 等。
  • 做法：提取核心集或位点面板的基因型矩阵，按群体结构/地理来源注释行/列，绘制聚类热图与树。
  • 诊断性证据：核心集应在热图与 PCA 上“均匀撑开”原始群体的结构，而非集中在某一两个亚群。

4.2 二维码/条码化交付

1
2
# 示例：将某样本的核心 SNP 串联成字符串后生成二维码
echo "ID001_VarietyName_ATGCCG..." | qrencode -o ID001_QR.png
  • 输入：编码了样本 ID / 品种名 / 关键位点型的字符串。
  • 输出：.png 二维码，可贴在资源袋或出库单上，作为“物理世界的指纹索引”。

阶段性诊断：贝叶斯式“结果复盘”

在项目收尾前，建议在第 04 篇的验收指标基础上，显式做两类“备择假设检验”：

  • 代表性测试：计算核心集的等位基因覆盖度（NA/RA 等），例如“是否覆盖了 ≥95% 的稀有等位基因”。
  • 结构一致性测试：在同一 PCA 空间中同时绘制全群体与核心集，观察核心集是否像“伞骨”一样均匀撑开整体结构、是否存在明显空白亚群。

4. 核心集筛选（Core Hunter 3 / 其他）

3.1 VCF 基础过滤（示例：vcftools / bcftools）

（示例，按项目替换）

  • 位点缺失率/样本缺失率过滤
  • MAF 过滤
  • 仅保留双等位 SNP（若需要）

3.2 PLINK：LD pruning（用于 PCA/结构推断的常见做法）

（示例，按项目替换窗口大小与 (r^2) 阈值）

  • 生成 pruned in 位点列表
  • 导出用于 PCA 的子矩阵

3.3 群体结构分析（至少一种）

PCA（R：SNPRelate，建议用于大矩阵）

  • 输入：VCF/GENO
  • 输出：PC1/PC2/…、样本坐标表、解释率

群体结构（ADMIXTURE，可选）

  • 输入：PLINK bed/bim/fam
  • 输出：Q 矩阵、交叉验证误差曲线（用于选 K）

树（FastTree，可选）

  • 输入：距离矩阵/比对/变异派生矩阵（按实现）
  • 输出：Newick 树文件

3.4 可选技术路线（按物种与数据形态选型）

3.4.1 比对：线性参考 vs 图泛基因组（Graph Pangenome）

  • 问题背景：结构变异（SV）丰富的物种中，强行比对到单一线性参考会引入“参考偏差”，远缘材料比对率与变异检测会系统性变差。
  • 选型建议：
    • 常规场景：短读长大规模项目，优先选择成熟线性参考方案（成本/复杂度更可控）。
    • 强 SV/远缘材料占比高：考虑图泛基因组路线（如 vg / minigraph 生态），以降低参考偏差（需同步评估工程复杂度与交付周期）。

3.4.2 变异检测：传统统计模型 vs 深度学习检测器

  • 问题背景：多倍体/高杂合群体中，重复序列与等位基因剂量不确定性会推高假阳性/假阴性风险。
  • 选型建议：
    • 常规 WGRS：以项目稳定性为先，选择主流统计模型工具链并加强过滤/验证。
    • 复杂基因组/高精度诉求：评估深度学习变异检测器（如 DeepVariant/DeepTrio）在目标物种数据上的可迁移性与成本；若采用，需固化模型版本与训练/微调口径。

3.4.3 缺失处理：基因型填充（Imputation）

  • 适用：低深度或高缺失矩阵；核心集构建与结构分析前常需提高矩阵完整性。
  • 常用工具：Beagle 等（若群体强分层，可采用“先聚类、后填充”的策略以提升速度与稳定性）。

4. 核心集筛选（Core Hunter 3 / 其他）

4.1 目标与约束（必须在项目里明确）

  • 目标 1（代表性）：优先优化 A-NE（见第 05 篇定义）
  • 目标 2（内部多样性）：优化 E-NE/内部距离类指标
  • 约束（可选）：PC/MC（按亚群/地理来源最低占比、至少 1 个代表等）

4.2 梯度筛选（推荐做法）

  • 做法：对多个核心集规模（如 5%、10%、15%、20%、…）分别筛选核心集并计算指标。
  • 选型：以 A-NE/等位基因覆盖/结构覆盖的“收益拐点”确定最终规模，同时满足业务约束（如必须覆盖特定地区/类型）。

4.3 输出（进入交付目录）

  • 核心集样本清单（每个规模一份）
  • 每个规模的指标表（与全体对比）
  • 核心集与全体的 PCA/聚类/树对比图

5. 关键参数说明（对结果影响最大）

  • 过滤阈值：缺失率、MAF、（可选）HWE、（可选）LD pruning 参数
  • 距离度量与输入矩阵：用遗传距离/表型距离/混合距离会改变优化方向
  • 优化目标权重：代表性 vs 多样性权重变化会影响核心集组成
  • 结构约束：PC/MC 约束会影响是否能覆盖小亚群/稀有类型
  • 随机种子/迭代次数：启发式优化需保证可重复

6. 与其他笔记的引用（避免重复）

  • 模拟退火/并行回火：见 模拟退火与并行回火.md
  • 指标体系（PIC/MAF/(H_e)/Shannon 等）与“代表性 vs 多样性”取舍：见第 05 篇