基于长读长测序与 cuteHap 的单倍型分型结构变异检测

一、背景

结构变异（Structural Variation, SV）指基因组上 ≥50 bp 的序列改变，包括缺失（DEL）、插入（INS）、倒位（INV）、重复（DUP）和易位（TRA）五种主要类型。SV 影响细胞过程、表型与疾病，是基因组研究绕不开的一层。长读长测序让 SV 检测和单倍型（haplotype）分析达到了前所未有的分辨率——而单倍型层面的信息对解释等位基因特异性效应、诊断隐性遗传病至关重要。

然而现有工具各有短板：

• 比对法检测器（cuteSV、Sniffles2 等）：在多等位基因区域识别不准，单倍型分析能力有限；
• 集成（ensemble）方法：能提升整体准确率，但天花板取决于底层检测器，单倍型分配同样存在缺陷；
• 从头组装法：分辨率最高，但计算资源消耗巨大，难以日常使用；
• 已有的相位整合（phasing-integrated）工具：跨平台兼容性差、计算负担重、单倍型信息利用不充分。

cuteHap 是为解决上述问题设计的单倍型感知 SV 检测框架。本 Protocol 给出一条端到端流程：SNV 检测 → 读段分型（read phasing）→ 单倍型分型 SV 检测与基因型判定，兼容 PacBio HiFi 与 ONT 等各类长读长平台的数据。文中以 HG002 样本的 PacBio HiFi（Revio）数据 + GRCh38 参考基因组为演示，其他个体的长读长数据同样适用。

创新点

cuteHap 重新设计了相位信息的利用方式，两大算法核心：

1. 单倍型特异性自适应聚类（haplotype-specific self-adaptive clustering）；
2. 簇可信度优先的束搜索（cluster credibility-prioritized beam search）。

两者配合，在不牺牲计算效率的前提下实现更优的 SV 发现与精确的分相基因型判定。整条流程不需要 GPU、不需要大规模集群，一台 32–64 GB 内存的多核服务器甚至高配台式机即可。

二、关键资源表

材料与设备：样本的人类参考基因组 FASTA 文件 + 长读长比对 BAM 文件；本流程所需软件（至少包括 Clair3、LongPhase、Samtools、cuteHap）。

三、开始之前（Before you begin）

1. 创建工作目录（耗时 1 min）

建立目录结构，用于存放临时文件与结果：

mkdir -p svwork
mkdir -p svwork/{ref,alns,snvs,phased,svs,benchmarks,models}

2. 下载数据集（耗时 1–12 h，视网速）

① 下载长读长测序比对文件（PacBio 或 ONT 平台）。本例使用已比对到 GRCh38 的公开 PacBio HiFi 数据：

wget -c -P alns https://downloads.pacbcloud.com/public/revio/2022Q4/HG002-rep1/analysis/HG002.m84011_220902_175841_s1.GRCh38.bam
wget -c -P alns https://downloads.pacbcloud.com/public/revio/2022Q4/HG002-rep1/analysis/HG002.m84011_220902_175841_s1.GRCh38.bam.bai

⚠️ CRITICAL：下一步下载的参考基因组必须与比对文件使用的参考基因组完全一致（本例两者均为 GRCh38）。

② 下载人类参考基因组及其索引：

wget -c -P ref http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/GRCh38_full_analysis_set_plus_decoy_hla.fa
wget -c -P ref http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/GRCh38_full_analysis_set_plus_decoy_hla.fa.fai

③ 下载 HG002 基准真值（ground truth）与高置信区间（用于后续评估）。本例使用基于 GRCh38 发布的 T2T Q100 版本：

wget -c -P benchmarks https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/data/AshkenazimTrio/analysis/NIST_HG002_DraftBenchmark_defrabbV0.019-20241113/GRCh38_HG2-T2TQ100-V1.1_stvar.vcf.gz
wget -c -P benchmarks https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/data/AshkenazimTrio/analysis/NIST_HG002_DraftBenchmark_defrabbV0.019-20241113/GRCh38_HG2-T2TQ100-V1.1_stvar.vcf.gz.tbi
wget -c -P benchmarks https://ftp-trace.ncbi.nlm.nih.gov/ReferenceSamples/giab/data/AshkenazimTrio/analysis/NIST_HG002_DraftBenchmark_defrabbV0.019-20241113/GRCh38_HG2-T2TQ100-V1.1_stvar.benchmark.bed

④ 下载 Clair3 对应的预训练模型。比对数据来自 HiFi Revio，故选 hifi_revio 模型：

mkdir models/hifi_revio
wget -c -P models/hifi_revio/ -r -np -nH --cut-dirs=3 -A .pt \
  https://www.bio8.cs.hku.hk/clair3/clair3_models_pytorch/hifi_revio/

注：Clair3 为不同平台和仪器提供多个预训练模型，必须与输入数据匹配。例如 ONT PromethION R9.4.1 流动槽 Super Accuracy（SUP）模式应选 r941_prom_sup_g5014，并相应修改目录与下载地址。

3. 安装软件、配置环境（耗时 1–2 h）

需要安装：Clair3（SNV 检测）、LongPhase（读段分型）、Samtools（比对索引）、cuteHap（单倍型分型 SV 检测）、Truvari（基准评估）。推荐用 conda/mamba 管理环境：

mamba create -n clair3_env   -c conda-forge -c bioconda -y clair3
mamba create -n longphase_env -c conda-forge -c bioconda -y longphase samtools
mamba create -n cutehap_env  -c conda-forge -c bioconda -y cutehap
mamba create -n truvari_env  -c conda-forge -c bioconda -y truvari samtools

四、分步方案

步骤 1：对样本进行 SNV 检测（耗时 3–6 h）

输入：长读长比对 BAM + 参考 FASTA；输出：全部检测到的 SNV（VCF）。

SNV（如 A→G）是基因组中最丰富的变异形式。正因为其常见且分布广泛，绝大多数读段分型工具依赖准确检出的 SNV 作为标记，来区分读段究竟来自父源还是母源单倍型。推荐用 Clair3（成熟的长读长 SNV 检测器）；也可换成 miniSNV、DeepVariant、Longshot 等。

conda activate clair3_env
run_clair3.sh \
  --bam_fn=alns/HG002.m84011_220902_175841_s1.GRCh38.bam \
  --ref_fn=ref/GRCh38_full_analysis_set_plus_decoy_hla.fa \
  --threads=16 --platform=hifi \
  --model_path=models/hifi_revio --output=snvs

注：--platform 与 --model_path 须与数据平台匹配。Clair3 支持 hifi 与 ont 两大平台，模型包括 hifi_revio、hifi_sequel2、r941_prom_hac_g360+g422、r941_prom_sup_g5014 等，随 Clair3 一起装在 ${CONDA_PREFIX}/bin/models/ 下，找不到可另行下载。

替代方案（以 miniSNV 为例）：

mamba install -c conda-forge -c bioconda -y whatshap
git clone https://github.com/CuiMiao-HIT/miniSNV.git
cd miniSNV/Release && make -j 12 && cd ../..
mkdir minisnv_dir
python miniSNV/chr_chunk_task.py \
  --fin_ref ref/GRCh38_full_analysis_set_plus_decoy_hla.fa \
  --fin_bam alns/HG002.m84011_220902_175841_s1.GRCh38.bam \
  --workDir minisnv_dir --human \
  --dup_bed miniSNV/GRCh38_segdups.bed \
  --homo_bed miniSNV/GRCh38_AllHomopolymers_gt6bp_imperfectgt10bp_slop5.bed

步骤 2：用检出的 SNV 对比对读段分型（耗时 10–30 min）

输入：同一份 BAM + 参考 FASTA + SNV VCF；输出：每条读段带 haplotag 标签的 BAM。

分型过程给每条读段打上 haplotag——标明它来自父源、母源，还是无法判定（unknown）。推荐 LongPhase（专为长读长数据设计，可高效利用 SNV 标记）；WhatsHap 亦可。

① phase 模式，对 Clair3 生成的 SNV 集合分相：

conda activate longphase_env
longphase phase -s snvs/merge_output.vcf.gz \
  -b alns/HG002.m84011_220902_175841_s1.GRCh38.bam \
  -r ref/GRCh38_full_analysis_set_plus_decoy_hla.fa -t 16 \
  -o phased/longphase --pb

注：--pb 对应 PacBio 数据；若为 ONT 数据，改成 --ont。

② haplotag 模式，给输入比对打标签：

longphase haplotag \
  -r ref/GRCh38_full_analysis_set_plus_decoy_hla.fa \
  -s phased/longphase.vcf \
  -b alns/HG002.m84011_220902_175841_s1.GRCh38.bam -t 16 \
  -o phased/longphase.tag

③ 为带标签 BAM 建索引：

samtools index phased/longphase.tag.bam

替代方案（以 WhatsHap 为例）：

mamba install -c conda-forge -c bioconda -y whatshap
whatshap phase snvs/merge_output.vcf.gz \
  alns/HG002.m84011_220902_175841_s1.GRCh38.bam \
  --reference=ref/GRCh38_full_analysis_set_plus_decoy_hla.fa \
  --distrust-genotypes --ignore-read-groups \
  -o phased/whatshap.vcf.gz
bcftools index phased/whatshap.vcf.gz
whatshap haplotag phased/whatshap.vcf.gz \
  alns/HG002.m84011_220902_175841_s1.GRCh38.bam \
  --reference=ref/GRCh38_full_analysis_set_plus_decoy_hla.fa \
  --ignore-read-groups -o phased/whatshap.tag.bam \
  --output-haplotag-list=phased/phasedinfo.tsv
samtools index phased/whatshap.tag.bam

步骤 3：用带标签比对进行单倍型分型 SV 检测（耗时 3–10 min）

输入：带 haplotag 的 BAM + 参考 FASTA；输出：含全部单倍型分型 SV 的 VCF。

SV 检测在分型后的 BAM 上进行。cuteHap 识别全部五种常见 SV 类型（DEL、INS、INV、DUP、TRA），并为每个 SV 判定分相基因型——因为输入是带 haplotag 的比对，输出天然是单倍型分型结果。

conda activate cutehap_env
cuteHap phased/longphase.tag.bam \
  ref/GRCh38_full_analysis_set_plus_decoy_hla.fa \
  svs/cutehap.vcf svs -s 5

注：-s 指定一个 SV 需要的最少支持读段数，低于阈值即过滤（排除测序错误与比对伪影导致的假阳性）。建议设为平均测序深度的约 1/6。cuteHap 还有其他常用参数，见表 1。

表 1：cuteHap 主要参数及默认值

五、预期结果

整个流程的设计目标是解决一个普遍痛点：常规 SV 检测方法无法区分两条同源染色体上的拷贝。借助长读长分相技术，本流程在准确率和单倍型分辨率上都优于传统 SV 检测——相位信息减少假阳性、并把变异正确分配到对应单倍型，最终得到父源、母源染色体各自独立的 SV 集合。

• 快：绕开了全基因组组装的高昂时间与算力开销，直接从比对文件出发；三个主步骤（SNV 检测 → 读段分型 → SV 检测）在典型数据集上一个自然日内可完成（不含下载与环境搭建的一次性开销）。
• 轻：兼容主流 CPU 计算环境，无需 GPU 或集群；32–64 GB 内存的普通多核服务器即可。
• 产出：
  • 步骤 1 输出 Clair3 结果目录，核心为压缩 VCF merge_output.vcf.gz 及索引 .tbi；
  • 步骤 2 输出带标签比对 longphase.tag.bam（每条读段记录含 HP 字段）及索引 .bai；
  • 步骤 3 输出最终结果——包含样本全部单倍型分型 SV 及分相基因型的 VCF 文件。

六、定量与统计分析（可选，耗时 30–60 min）

拿到单倍型分型 SV 集合（VCF）后，可用基准评估衡量检出质量。评估工具将检出集合与"基准真值"比对——本例为 HG002 T2T 组装经人工验证的高置信 SV 集合，并配合对应的高置信区间（ground truth 可信的基因组区间）。推荐使用广泛使用的 Truvari，计算 precision、recall、F1 score 与基因型一致率（genotype concordance）。

两个命令：truvari bench（必需，直接比对）与 truvari refine（可选，按 SV 表示方式的差异对 bench 结果做精细化调整，在做不同 caller 间对比时推荐使用）：

conda activate truvari_env
bgzip -c svs/cutehap.vcf > svs/cutehap.vcf.gz
tabix -p vcf svs/cutehap.vcf.gz
truvari bench -b benchmarks/GRCh38_HG2-T2TQ100-V1.1_stvar.vcf.gz \
  -c svs/cutehap.vcf.gz --includebed \
  benchmarks/GRCh38_HG2-T2TQ100-V1.1_stvar.benchmark.bed \
  -o truvari_cmp --pick ac --passonly -r 2000 -C 5000 \
  -f ref/GRCh38_full_analysis_set_plus_decoy_hla.fa
truvari refine --use-original-vcfs --align mafft \
  -f ref/GRCh38_full_analysis_set_plus_decoy_hla.fa \
  -r truvari_cmp/candidate.refine.bed truvari_cmp -t 16 \
  --mafft-params "--auto --thread 16"

注：-b 与 --includebed 指定真值。真值 VCF 与高置信 BED 必须与被评估样本、以及全流程使用的参考基因组一致——否则指标毫无意义。

输出目录 truvari_cmp 中包含：检出集合的真阳性/假阳性、真值集合的真阳性/假阴性（均为压缩 VCF，可逐个 SV 查看），以及两个 JSON 指标文件——summary.json（bench 直接比对结果：precision、recall、F1、基因型一致率、TP/FP/FN 计数）与 refine.variant_summary.json（refine 精修后的 precision、recall、F1）。

七、局限性

1. 纯命令行：本流程全部通过 Linux 命令行运行，用户友好但无图形界面。
2. 依赖工具选择：整体性能取决于所选生信工具。文中选用 Clair3、LongPhase、cuteHap 三件套，但同类先进工具可直接替换，预期结果无显著差异——流程不绑定特定软件组合。
3. 仅适用二倍体：流程以人类样本为例，但同样适用于任何二倍体样本（如其他动物或植物）；多倍体不适用——现有分相工具按"两条单倍型"设计，三条及以上单倍型无法准确分相。
4. 仅适用长读长数据：短读长分相工具虽存在，但短读长无法可靠跨越长而复杂的结构变异，故本流程不支持短读长比对。

八、Troubleshooting

问题 1：下载的输入文件损坏或不完整（对应"下载数据集"）

解决：下载时始终加 -c 选项——wget 会检查本地已有的部分文件，从断点续传而非重新开始。下载完成后，与数据源提供的 MD5 校验值（如有）比对，确认文件完整。

问题 2：Clair3 环境检查后立即退出，报 WARNING: "No contig intersection found, output header only in snvs/merge_output.vcf.gz"（对应步骤 1）

解决：输入参考基因组与生成比对文件时的参考基因组不一致（contig 名称或序列不匹配）。确保下载的参考版本与 reads 最初比对所用版本完全相同。不确定时，用 samtools view -H 查看 BAM 头部——其中通常记录了 contig 名称与基因组信息。参考文件匹配是 Clair3 找到重叠 contig、正常产出 SNV 的前提。

问题 3：编译 miniSNV 时报错 "form.hpp:35:5: error: 'uint32_t' does not name a type"（对应步骤 1 替代方案）

解决：miniSNV 原基于 GCC 4.8.5 开发，新版 GCC 编译更严格导致该类型未识别。用文本编辑器打开 miniSNV/src/cmtools/form.hpp，在文件顶部（其他 #include 旁）加入一行 #include <cstdint>，保存后重新 make -j 12 即可。

问题 4：运行 cuteHap 报 KeyError: 'GL000008.2'（对应步骤 3）

解决：提供给 cuteHap 的参考基因组不含 contig GL000008.2，但 BAM 中的 reads 有比对到该 contig——即参考与生成 BAM 时的参考不完全一致。确保 cuteHap 所用 FASTA 与 BAM 头部出现的全部 contig 一致；若个别 contig 仍缺失，要么下载完整参考，要么用一致的完整参考重新比对。

问题 5：运行 cuteHap 报 FileNotFoundError: [Errno 2] No such file or directory: 'xxx.pickle'（对应步骤 3）

解决：多线程并行时的竞态条件所致。减少线程数重跑即可——文件访问冲突消除后错误通常不再出现，且不影响最终结果。

问题 6：运行 truvari refine 报 "[E::fai_build_core] File truncated at line 1 / pysam.utils.SamtoolsError: 'samtools returned with error 1: ... Could not build fai index /tmp/xxxx.fa.fai'"（对应定量分析）

解决：truvari refine 调用 pysam.samtools.faidx 提取参考子序列；pysam ≥0.20 版本把 FASTA 输出导向 stdout 而非文件，大段 FASTA 直接打印到终端导致失败。降级 pysam 可绕过，但更推荐安装 Truvari 的 develop 版（已修复且不改其他依赖）：

mamba create -n truvari_develop -c conda-forge -c bioconda -y pyabpoa mafft
git clone https://github.com/ACEnglish/truvari.git
cd truvari && python3 -m pip install .

九、资源与数据可用性

• 分析用源数据：PacBio 云数据库（downloads.pacbcloud.com/public/revio/2022Q4/HG002-rep1/analysis/）与 1000 Genomes FTP；所有工具均开源。
• 本文产出的 SV 集合：Zenodo（doi.org/10.5281/zenodo.19728430）。
• 基金：国家自然科学基金（62472120）、国家重点研发计划（2024YFC3406303）、中央高校基本科研业务费（XNJKKGYDJ2026010）等。