1. 长读长单倍型分相工程导论

人类是二倍体生物，绝大多数常染色体都拥有两份来源不同的拷贝：一份来自父亲，一份来自母亲。这两份拷贝之间平均每 1000 个碱基就存在约 1 个差异位点，这些差异构成了个体的单倍型（haplotype）。然而，常规的二代测序（NGS）读长只有 100–300 bp，远短于基因组的重复区与结构变异跨度，导致测序得到的短读在比对后无法判断某个杂合位点究竟属于父源还是母源拷贝。这种"知道有差异、却不知道差异如何连锁"的状态，就是分相问题的核心。

本课程的目标，是让学习者掌握从长读长测序数据出发，构建高质量、染色体级别单倍型分型的完整工程能力。所谓"工程"，意味着我们不仅关心某个算法能否跑通，更关心在真实样本、真实覆盖度、真实错误率条件下，如何设计流程、选择工具、评估质量、排查异常，并最终交付可被下游分析（如等位基因特异性表达、疾病关联、群体遗传）直接使用的分相结果。

长读长测序技术概览

长读长测序是分相工程的物理基础。当前主流平台可分为两大类：以 PacBio 为代表的单分子实时测序（SMRT）和以 Oxford Nanopore Technologies（ONT）为代表的纳米孔测序。二者都能产生平均长度 10–20 kb、最长可超过 100 kb 的读段，足以跨越大多数重复序列和结构变异。

PacBio 的 HiFi 模式通过环形一致性测序（CCS）对同一分子多次读取后取共识，单次准确率可达 99.9% 以上，读长通常为 10–25 kb。HiFi 读段的高准确率使其在比对和变异检测阶段几乎不依赖复杂的纠错，是当前分相工程中最"省心"的输入。ONT 的读长优势更明显，超长读可超过 100 kb，但原始单次准确率约 90–97%，需要配合纠错或使用对错误容忍的分相算法。选择哪种平台，取决于目标区域的大小、样本质量以及预算。

无论哪种平台，分相工程都依赖两个关键信号：杂合位点和读段跨越关系。当一个读段同时覆盖两个或多个杂合位点时，这些位点就被"连锁"在一起，形成一条局部单倍型片段。读段越长、覆盖度越高，连锁的位点越多，单倍型块（haplotype block）就越长。

单倍型分相的生物学意义

分相不是纯粹的计算游戏，它直接决定了许多生物学结论的正确性。在医学遗传学中，复合杂合突变（compound heterozygous）的两个致病位点分别位于父源和母源拷贝上时才会致病；如果错误地分到同一条单倍型上，就会把致病判断变成误判。在等位基因特异性表达分析中，只有知道每个读段来自哪条单倍型，才能准确计算两个等位基因的表达比例。在群体遗传与进化研究中，单倍型是重组、选择、基因流分析的基本单位。

此外，分相还能帮助组装端粒到端粒（T2T）级别的二倍体基因组，识别大规模结构变异和拷贝数变异的等位基因归属，这些都是短读技术难以完成的任务。

工程化分相流程的整体架构

一个可复现的分相工程，通常包含五个相互衔接的阶段：数据质控与比对、杂合位点识别、读段连锁与分相、单倍型块构建与定相、质量评估与结果交付。下图展示了这一整体架构以及各阶段之间的数据流关系。

流程的起点是长读 FASTQ 数据。以 ONT 为例，典型样本要求覆盖度 20×–30×，读长 N50 不低于 15 kb；若使用 HiFi，覆盖度 15×–20× 即可获得较稳定的分相效果。比对阶段推荐使用 minimap2 的 map-ont 或 map-hifi 预设，对高度重复区域可改用 winnowmap 以降低错配。比对完成后，用 samtools 排序并建立索引，得到 BAM 文件。

杂合位点识别阶段，HiFi 数据可直接用 DeepVariant 或 bcftools call 生成 VCF，并过滤掉低深度、低质量及偏离哈迪-温伯格平衡的位点。ONT 数据则建议先用 Clair3 等对长读友好的变异检测工具。最终保留的杂合位点数量通常在 150 万–250 万个之间（以人类基因组为例），这些位点就是后续分相的"锚点"。

分相阶段是工程的核心。以 whatshap 为例，其基本命令为：

whatshap phase \
  --reference GRCh38.fa \
  --output phased.vcf \
  --indels \
  sample.vcf.gz \
  sample.bam

该命令会遍历每个读段，根据读段上覆盖的杂合位点等位基因，将读段分配到单倍型 H1 或 H2，再通过图优化算法求解全局一致的分相方案。对于读长较短或覆盖度不足的区域，可启用 --ped 参数引入家系信息，或使用 longshot 进行基于读段重叠的局部定相。

最后是质量评估。分相质量的核心指标有两个：switch error rate（切换错误率）和单倍型块 N50。切换错误率衡量相邻杂合位点被错误分配到不同单倍型的比例，理想值应低于 1%；单倍型块 N50 则反映分相结果的连续性，人类样本中优秀结果可达数 Mb 甚至染色体级别。评估时可用 whatshap stats 生成统计报告，并与已知真值（如 GIAB 样本）比对验证。

需要强调的是，分相工程不是一次性的线性任务。当质量评估发现某些区域切换错误率偏高时，往往需要回到比对或变异检测阶段，检查是否存在比对错误、低覆盖或假杂合位点，然后重新分相。这种"评估—反馈—优化"的迭代，正是工程化与简单脚本调用的本质区别。