核心种质筛选-03操作流程

使用CoreHunter3提取core-set的分析流程。

分析准备

工具清单

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

上游数据处理与清洗

数据输入规范

已有 VCF/PLINK/矩阵：直接进入下文的“数据准入与过滤”阶段。

只有 FASTQ（需先做变异检测）：先完成以下上游步骤，再进入下文流程： 1) FASTQ 质控与过滤
2) 比对到参考（线性参考/图泛基因组按物种与材料差异选型）
3) 变异检测与联合分型
4) 变异过滤 → 得到 filtered.vcf 或 PLINK 文件

数据准入与基础过滤

这一阶段的目标是得到一个“能用于距离计算与优化”的稳定基因型矩阵。常见过滤项（按项目取舍）包括：

位点缺失率/样本缺失率过滤
MAF 过滤
（可选）仅保留双等位 SNP
（可选）HWE 或批次相关过滤（取决于数据类型与项目规范）

缺失处理（可选：高缺失/低深度时再做）

缺失过高会影响遗传距离与 EN/A-NE 等指标。若后续要用距离矩阵作为 Core Hunter 输入，建议在这里把缺失率控制到可接受范围。

步骤	工具	关键参数（示例）	输入/输出	说明
基因型填充	`Beagle 5.4` 等	`gt=filtered.vcf` `nthreads=10`	入：`filtered.vcf`；出：`imputed.vcf.gz`	分层很强时可先按结构/来源分组，再分别填充，提高稳定性与速度。

群体结构分析与矩阵构建

理解样本群体结构

这一步主要用于“看清楚数据”，顺便为后续可能的结构约束（例如每个亚群至少保留一定比例）提供依据。

LD pruning（可选）：用于 PCA/ADMIXTURE 等结构推断时，减少强 LD 对降维与聚类的影响。
PCA（建议）：
- 输入：VCF/GENO/PLINK（按实现）
- 输出：样本 PC 坐标表、解释率、作图用注释表（地区/亚群/批次等）
ADMIXTURE（可选）：
- 输入：PLINK bed/bim/fam
- 输出：Q 矩阵、交叉验证误差曲线（用于选 K）
树（可选）：
- 输入：距离矩阵或派生矩阵（按实现）
- 输出：Newick 树文件

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

如果用 Core Hunter 走“距离驱动”的优化路线，距离矩阵通常是核心输入之一。

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
输出：dist_matrix（数值型对称矩阵）

注：或通过plink来计算IBS距离矩阵。

# 1. 将 VCF 转换为 PLINK 二进制格式 (BED/BIM/FAM)
# --allow-extra-chr 强制 PLINK 接受非人类物种的染色体或 Contig 命名
# --double-id 确保样本 ID 被正确解析
plink --vcf imputed.vcf.gz --double-id --allow-extra-chr --make-bed --out plink_bin

# 2. 计算基于状态一致性 (IBS) 的距离矩阵
# --distance square 默认计算 1 - IBS，并输出一个完整的 N x N 对称矩阵
plink --bfile plink_bin --allow-extra-chr --distance square --out plink_dist

# 1. 读取样本 ID
# V1 通常包含样本的实际 ID (PLINK 的 Family ID 视 --double-id 参数而定，此处取 V1)
sample_ids <- read.table("plink_dist.mdist.id", header = FALSE, stringsAsFactors = FALSE)$V1

# 2. 读取纯数值矩阵并转换为 R 的 matrix 对象
dist_matrix <- as.matrix(read.table("plink_dist.mdist", header = FALSE))

# 3. 赋予矩阵精确的行名和列名
rownames(dist_matrix) <- sample_ids
colnames(dist_matrix) <- sample_ids

# 此时，dist_matrix 已经是一个数值型对称矩阵，可直接输入给 Core Hunter

核心集筛选与操作（Core Hunter 3）

核心种质筛选不仅仅是跑通代码，更核心的是分析策略的制定。有几个关键指标需要重点关注，我们可以通过组合目标函数并调节指标权重，从而得到符合业务需求的目标核心群体。

Core Hunter 3 数据结构与读入

Core Hunter 3 支持多种数据类型（基因型、表型、预计算距离矩阵）独立或组合输入，最终都可以打包成一个 chdata 对象供核心集筛选使用。

文件与通用约定（官网 Data 页）

文件格式：支持逗号分隔 csv 与制表符分隔 txt（GUI/命令行从文件读入；R 包也可从 data.frame/matrix 读入）。
必需列：几乎所有“表格型输入”都要求首列为 ID（唯一的 accession 标识符）。
可选列：可加第二列 NAME（不必唯一，也允许部分缺失）。
引号与空白（所有格式通用）：
- ID、NAME 以及任何文本值都可以用单引号 ' 或双引号 " 包裹；未被引号包裹的两侧空白会被去掉。
- 若文本内含单引号，则需用双引号包裹；若文本内含双引号，则需用单引号包裹。
- 同时包含单引号与双引号的值不被允许（无法用单一外层引号安全包裹）。

1. 基因型数据 (Genotypes) 使用 genotypes() 函数读入。支持三种格式：

default：默认格式（每个 marker 可多列）。一行一个样本、每个 marker 对应 1 列或多列“等位观测值”（例如二倍体每个 marker 两列）。值可以是等位基因名称/编号，或任何用于标识等位型的 token。不适用于 bulk（混样/池化）。
- 表头规则：首列必须 ID，可选 NAME；剩余列为 marker 列。
- marker 列命名：同一 marker 的多列可以在表头里重复 marker 名；也允许用 . / - / _ 加后缀区分多列（如 mk1.a/mk1.b、mk1-1/mk1-2、X_1/X_2）。Core Hunter 会将最后一次出现的 ./-/_ 之前的前缀识别为 marker 名。
- 示例（diploid：二倍体，每 marker 两列）：

ID	mk1.a	mk1.b	mk2.a	mk2.b	mk3.a	mk3.b	mk4.a	mk4.b
A	1	3	B	B	a1	a1	-	+
B	2	2	C	A	a1	a2	+	-
C	1	2	D	D	a2	a2	+	+
D	2	3	B	B	a2	a1	+	-
E	1	1	A	A	a1	a1	-	-

示例（含 NAME 且存在缺失）：

ID	NAME	mk1-1	mk1-2	mk2-1	mk2-2	mk3-1	mk3-2	mk4-1	mk4-2
A	Alice	1	3	B	B	a1	a1
B	Bob	2	2	C	A	a1	a2	+	-
C	Carol	1	2	D	D	a2	a2	+	+
D	Dave	2	3	B	B	a2	a1	+	-
E	Eve	1	1	a1	a1	-	-

示例（homozygous/haploid：每 marker 一列）：

ID	mk1	mk2	mk3	mk4
A	1	B	a1	-
B	2	C	a1	+
C	1	D	a2	+
D	2	B	a2	+
E	1	A	a1	-

frequency：等位基因频率格式（marker×allele 一列）。每行一个样本、每列对应“marker×allele”的组合，单元格为频率值；对同一个样本，在同一个 marker 下所有 allele 的频率应当加和为 1。适用于个体样本与 bulk 样本。
- 表头规则：首列必须 ID，可选 NAME；第一行表头写 marker 名（每个 marker 会跨多列重复出现）。
- 可选 ALLELE 头行：可在数据前追加一行（行首为 ALLELE）给每一列写 allele 名。
- 示例（不含 ALLELE 行）：

ID	mk1_1	mk1_2	mk1_3	mk2_1	mk2_2	mk3_1	mk3_2	mk3_3
A	0.50	0.50	0.00	0.50	0.50
B	1.00	0.00	0.00	0.50	0.50	0.00	0.50	0.50
C	0.60	0.00	0.40	0.50	0.50	0.00	0.50	0.50
D	1.00	0.00
E	0.33	0.33	0.33	0.50	0.50	0.00	0.50	0.50

示例（含 ALLELE 行）：

ID	mk1_1	mk1_2	mk1_3	mk2_1	mk2_2	mk3_1	mk3_2	mk3_3
ALLELE	a	b	c	+	-	1	2	3
A	0.50	0.50	0.00	0.50	0.50
B	1.00	0.00	0.00	0.50	0.50	0.00	0.50	0.50
C	0.60	0.00	0.40	0.50	0.50	0.00	0.50	0.50
D	1.00	0.00
E	0.33	0.33	0.33	0.50	0.50	0.00	0.50	0.50

示例（含 NAME + ALLELE，只展示前两组 marker）：

ID	NAME	mk1_1	mk1_2	mk1_3	mk2_1	mk2_2
ALLELE		a	b	c	+	-
A	Alice	0.50	0.50
B	Bob	1.00	0.00	0.00	0.50	0.50
C	Carol	0.60	0.00	0.40	0.50	0.50
D	Dave	1.00	0.00
E	Eve	0.33	0.33	0.33	0.50	0.50

biparental：二亲/双等位计数格式（0/1/2 计数）。每行一个样本、每个 marker 1 列，取值 0/1/2 表示“某个参考等位基因出现次数”。0 与 2 表示两种纯合，1 表示杂合。仅适用于每 marker 至多两种等位的情况，不适用于 bulk。
- 表头规则：首列必须 ID，可选 NAME；表头可写 marker 名（也可以省略，由软件按列序号处理，官网示例是提供 marker 名）。
- 示例：

ID	mk1	mk3	mk4	mk5	mk6	mk7
A	1	2	1	1	0	0
B	2	2	0	1	2	1
C	1	0	1	1	0
D	1	1	1	1	2
E	1	0	2	0

示例（含 NAME）：

ID	NAME	mk1	mk3	mk4	mk5	mk6	mk7
A	Alice	1	2	1	1	0	0
B	Bob	2	2	0	1	2	1
C	Carol	1	0	1	1	0
D	Dave	1	1	1	1	2
E	Eve	1	0	2	0

# 示例：从文件读入双等位 SNP 数据
geno <- genotypes(file = "genotypes-biparental.csv", format = "biparental")

2. 表型数据 (Phenotypes) 使用 phenotypes() 函数读入，一行一个样本、一列一个性状，首列必须 ID（可选 NAME）。

强烈建议显式写 TYPE 行：可在表头下增加一行（行首为 TYPE），为每个性状指定变量类型与（可选）数据类型。若不提供，变量类型会尝试推断；但在 GUI 中，常会把未注明类型的列默认当作“nominal 字符串”处理，因此建议总是写 TYPE。
- 变量类型（type）编码：

变量类型	Code	默认数据类型
Nominal（名义）	N	String
Ordinal（序数）	O	Integer
Interval（区间）	I	Integer
Ratio（比率）	R	Double

数据类型（data type）编码（可附加在类型后，例如 NB、NI）：

数据类型	Code
Boolean	B
Short	T
Integer	I
Long	L
Float	F
Double	D
Big Integer	R
Big Decimal	M
Date	A
String	S
None	X

可选 MIN/MAX 行：对 O/I/R（序数/区间/比率）性状，可再增加 MIN 与 MAX 行提供“指示性范围”。如果不写，会由数据推断；若数据超出所写范围，会自动把范围扩展到实际的 min/max。
示例（含 TYPE 行；包含 N/O/I/R/NB）：

ID	trait 1	trait 2	trait 3	trait 4	trait 5
TYPE	N	O	I	R	NB
A	A	3	4	1.4	false
B	B	1	5	0.5	true
C	A	0	6	0.5	true
D	C	2	9	0.5	false
E	B	2	1	1.3	true

示例（含 NAME）：

ID	NAME	trait 1	trait 2	trait 3	trait 4	trait 5
TYPE	N	O	I	R	NB
A	Alice	A	3	4	1.4	false
B	Bob	B	1	5	0.5	true
C	Carol	A	0	6	0.5	true
D	Dave	C	2	9	0.5	false
E	Eve	B	2	1	1.3	true

示例（显式 MIN/MAX）：

ID	NAME	trait 1	trait 2	trait 3	trait 4	trait 5
TYPE	N	O	I	R	NB
MIN			0	0.0
MAX			10	2.0
A	Alice	A	3	4	1.4	false
B	Bob	B	1	5	0.5	true
C	Carol	A	0	6	0.5	true
D	Dave	C	2	9	0.5	false
E	Eve	B	2	1	1.3	true

# 示例：从文件读入表型数据
pheno <- phenotypes(file = "phenotypes.csv")

3. 距离矩阵 (Distances) 如果通过其他软件（如 SNPRelate、PLINK）预先计算好距离矩阵，可直接作为输入。

矩阵形状：一行一列一个样本（同一顺序），本质上应为对称矩阵（距离矩阵一般是对称的）。
表头与行名：首列必须为 ID（可选 NAME）；列名可以选择性地重复写 ID。
可截断：每一行允许在对角线处或之后截断（对角线元素若出现应为 0）；若提供上三角的一部分，会校验其与下三角一致性；缺失的上三角会用下三角拷贝补全。
示例（完整矩阵）：

ID	A	B	C	D	E
A	0.0	0.2	0.4	0.6	0.8
B	0.2	0.0	0.2	0.4	0.6
C	0.4	0.2	0.0	0.1	0.4
D	0.6	0.4	0.1	0.0	0.2
E	0.8	0.6	0.4	0.2	0.0

示例（含 NAME）：

ID	NAME	A	B	C	D	E
A	Alice	0.0	0.2	0.4	0.6	0.8
B	Bob	0.2	0.0	0.2	0.4	0.6
C	Carol	0.4	0.2	0.0	0.1	0.4
D	Dave	0.6	0.4	0.1	0.0	0.2
E	Eve	0.8	0.6	0.4	0.2	0.0

示例（截断到对角线；上三角可省略）：

ID	NAME	A	B	C	D	E
A	Alice	0.0
B	Bob	0.2	0.0
C	Carol	0.4	0.2	0.0
D	Dave	0.6	0.4	0.1	0.0
E	Eve	0.8	0.6	0.4	0.2	0.0

# 示例：直接读入距离矩阵
dist <- distances(file = "distances.csv")

4. 组合多维度数据 实际业务中，我们常需要兼顾基因型和表型。使用 coreHunterData() 可以将多种数据打包对齐（要求各数据的样本 ID 必须一致）：

# 整合基因型与表型数据为一个综合对象
my_data <- coreHunterData(genotypes = geno, phenotypes = pheno)

算法

Core Hunter 3 的底层搜索算法默认使用高级的平行回火算法（Parallel Tempering），也可切换为快速的随机爬山算法（Stochastic Hill-climbing）。相较于早期版本，其在计算效率上有显著提升，同时能兼顾运算速度与遗传变异覆盖的完整度。

(注：在一些文献中常简写的 “EN-MR” 实际上是指以 Entry-to-Nearest (EN) 为优化目标，并采用 Modified Roger’s distance (MR) 作为遗传距离度量的组合策略，而非某种特定的搜索算法名称。)

目标函数简称与组合策略

在构建核心集时，Core Hunter 允许我们组合不同维度的目标函数。我们通常采用组合目标（多目标优化）以兼顾全面：

EN (Entry-to-Nearest)：核心集内部个体到最近邻个体的距离。最大化此距离可提高核心集内部的多样性（Diversity），避免冗余。
AN (Accession-to-Nearest)：全集样本到核心集最近个体的距离。最小化此距离可提高核心集的代表性（Representativeness），确保所有原始样本在核心集中都能找到相似的代表。
MR (Modified Roger’s distance) / GD (Gower’s distance) / PD (Precomputed Distance)：
- MR：改进的罗杰斯距离，是 Core Hunter 处理基因型 (Genotypes) 数据的默认距离度量。
- GD：Gower 距离，是处理表型 (Phenotypes) 数据的默认距离度量。
- PD：预计算距离矩阵 (Precomputed Distance)，当传入外部算好的距离矩阵时使用。

经典组合策略示例： 若希望兼顾表型多样性（EN）与遗传代表性（AN）：

EN + GD：在表型层面（GD），最大化内部多样性（EN）。
AN + MR：在遗传层面（MR），最大化全集代表性（AN）。

(关于这些目标函数的详细算法原理与权衡逻辑，请参考 [05篇：底层算法])

library(corehunter)

# 假设我们已经按照上文组合好了 my_data 对象（包含基因型与表型数据）
# 示例：采用组合策略，兼顾表型多样性与遗传代表性
core <- sampleCore(my_data, obj = list(
  objective("EN", "GD", weight = 0.3), # 表型 (Gower's distance) 的多样性
  objective("AN", "MR", weight = 0.7)  # 遗传 (Modified Roger's distance) 的代表性
))

core_ids <- core$sel  # 被选中的样本索引或 ID

(注：虽然在一些育种文献中 MR 指 Marker Retention（标记保留率），但在 Core Hunter 的参数与目标函数语境下，MR 指代的是 Modified Roger’s distance，这是一种遗传距离算法。)

核心集规模选择（梯度规模模拟，推荐）

建议对多个核心集规模分别筛选并计算指标，再结合业务约束确定最终规模。

做法：5%、10%、15%、20%…分别跑核心集筛选
对比：每个规模输出一份核心清单 + 一张指标表 + 一组对比图
决策：看代表性/覆盖类指标的收益拐点，同时满足必须覆盖的亚群/地区/类型约束

关键参数总结（对结果影响最大）

过滤阈值：缺失率、MAF、（可选）HWE、（可选）LD pruning 参数
核心集比例（size）：常见 0.05–0.20，建议用梯度规模模拟寻找指标的拐点。
距离度量与输入矩阵：用遗传距离/表型距离/混合距离会改变优化方向
优化目标与权重：代表性 vs 多样性权重变化会影响核心集组成
结构约束：按亚群/地区的最低占比、至少 1 个代表等约束会改变结果
随机种子/迭代次数：启发式优化需要保证可重复

验证评估与下游应用

验证与评估

跑完工具拿到核心集列表后，通常还需要进行质量验证与可视化，证明“选得对”。

指标计算：需要计算一系列指标（如等位基因覆盖度、多样性保留率等）来量化评估结果。(详细的验收标准与指标解读，请见 [04篇：结果表征与业务解读]；底层的距离定义见 [05篇：底层算法])
可视化验证：
1. 结构一致性检查 (PCA/ADMIXTURE)：在同一 PCA 空间中绘制全体与核心集，检查核心集是否均匀覆盖所有独立亚群。
2. 指纹图谱与热图：提取核心集的基因型矩阵绘制聚类热图，直观展示多样性保留情况。

(关于交付物标准目录结构，请参见 [04篇：结果表征与业务解读])

下游应用（可选）

指纹图谱标记开发与筛选

DNA 是遗传物质，携带有生物材料的完整的遗传信息和物种个体间的差异信息。生物材料的 DNA 差异可以稳定遗传，受外界条件的干扰较少，随着生物测序技术的快速发展，各种分子标记被开发出来用于种质筛选、品种鉴定、产权保护等。这些标记包括 RFLP、AFLP、SSR、InDel 和 SNP 等，与形态标记相比，基于 DNA 序列差异的分子标记具有以下优点：

不受环境影响；
数量多、分布广；
多态性高、自然材料存在丰富的变异；
多数标记类型表现为共显性，能够鉴定出纯合基因型和杂合基因型，为育种材料的选择和鉴定提供了极大的便利。

DNA 指纹图谱标记由少数具有较强代表性的标记构成，这些标记组合可以对同种内的不同材料进行区分，结果稳定，鉴定方便，因此在作物种质资源多样性的研究、品种产权确定和保护等方面广泛应用。

DNA 指纹图谱标记开发主要目的在于对推广品种进行产权确认和保护、对生产种子进行纯度鉴定，所以用于标记开发的材料应当包含当前生产使用的骨干育种亲本、重要的育种资源材料及育种应用潜力较大的野生材料，以保证筛选的位点的代表性。

标记筛选

指纹图谱标记应尽量避免冗余，基于非遗传群体材料构成的群体开发的 SNP，应当具有标记质量高、代表性强、标记（组合）的材料区分度高、在基因组上分布均匀、特异性强等特点。基于这些原则进行筛选（Wang et al., 2021）：

标记在基因组上分布均匀；
标记不存在缺失位点，即位点完整度 100%；
丢弃最小等位基因频率（Minor Allele Frequency, MAF）低于 20% 的位点；
丢弃多态信息含量（PIC）小于 0.35 的位点；
哈迪-温伯格检验，保留 p-value 大于 0.01 的位点；
相邻标记距离大于 250bp。

精准育种与 GS

把核心集上学到的“基因型→表型”关系外推到全群体，实现不测表型的预测筛选。

应用步骤	工具/函数	简化模型	输入/输出	说明
模型训练	R 包 `rrBLUP::mixed.solve` 等	$y = X\beta + Zu + e$	入：核心集基因型 + 表型；出：SNP 效应值 $\beta$	核心集代表性越好，外推越稳。
基因组预测	矩阵乘法	$GEBV = G_{\text{all}} \times \beta$	入：全群体基因型矩阵 $G_{\text{all}}$；出：每个个体的 GEBV	用于候选筛选与排序。

资源库管理与溯源（二维码/条码）

# 示例：将某样本的核心 SNP 串联成字符串后生成二维码
echo "ID001_VarietyName_ATGCCG..." | qrencode -o ID001_QR.png

输入：编码了样本 ID / 品种名 / 关键位点型的字符串
输出：.png 二维码，可贴在资源袋或出库单上