PacBio HiFi测序分析¶

一句话概述¶

PacBio HiFi (High-Fidelity) 测序通过环形一致性序列(CCS)技术实现"又长又准"——白话说就是：让DNA分子在环形模板上被反复读取多次，多次结果取"投票共识"，最终得到10-25kb的长读长且准确率达99.9%的序列。

核心知识点表格¶

知识点	说明
CCS原理	同一DNA模板被反复测序（多pass），取共识序列，错误率降到<0.1%
HiFi reads	CCS产生的高保真reads，读长10-25kb，准确率≥99.9% (Q30)
CLR reads	单pass长读长（25-175kb），错误率8-15%，速度快
pbmm2	PacBio官方比对工具，本质是minimap2的wrapper
hifiasm	最主流的HiFi基因组组装工具，支持单倍型分离
DeepVariant	Google开发的AI变异检测工具，有HiFi专用模型
pbsv	PacBio官方结构变异检测工具
HiPhase	单倍型定相工具
Revio / Vega	PacBio最新测序仪，Revio高通量，Vega是短读长+长读长混合
适用场景	T2T基因组组装、结构变异、全长转录本、表观修饰

各步骤详解¶

第一步：理解HiFi测序原理¶

白话解释： PacBio测序时，DNA分子两端连上接头变成"甜甜圈"（SMRTbell）。测序酶沿着环走，每走一圈就读一次。一个15kb的DNA模板可能走5-10圈。虽然每一圈有错误（~15%），但多圈取"共识"后，错误率就降到了0.1%以下——就像让10个人各自抄一遍课文，虽然每个人都会抄错几个字，但把10份对比取多数结果，就几乎没错了。

PacBio两种测序模式对比：
CLR:  1次扫描 → 25-175kb → 错误率8-15% → 量大但不准
HiFi: 多次扫描 → 10-25kb  → 错误率<0.1% → 又长又准

第二步：CCS生成（subreads → HiFi reads）¶

白话解释： 把原始测序仪输出的多轮扫描数据（subreads）合并成高质量的HiFi reads。

# 安装PacBio工具集
conda install -c bioconda pbccs pbmm2 pbsv  # 安装核心工具

# 从subreads生成CCS/HiFi reads
# 输入：subreads.bam（测序仪原始输出）
# 输出：ccs.bam（HiFi reads）
ccs \                                    # CCS命令
  subreads.bam \                         # 输入subreads
  ccs.bam \                              # 输出HiFi reads
  --min-rq 0.99 \                        # 最低质量阈值Q20（99%准确率）
  --min-passes 3 \                       # 最少扫描3圈
  --num-threads 32                       # 32线程

# 注意：新版Revio测序仪可以直接输出HiFi reads
# 不需要再跑ccs命令

# 查看CCS统计信息
python -c "
import pysam                             # 导入pysam读取BAM
bam = pysam.AlignmentFile('ccs.bam')     # 打开BAM
reads = list(bam)                        # 读取所有reads
print(f'Total reads: {len(reads)}')      # 总读数
lengths = [r.query_length for r in reads]  # 读长列表
print(f'Mean length: {sum(lengths)/len(lengths):.0f}')  # 平均读长
print(f'N50: {sorted(lengths, reverse=True)[len(lengths)//2]}')  # N50
"

第三步：比对（pbmm2）¶

白话解释： pbmm2是PacBio对minimap2的封装，默认参数针对PacBio数据做了优化。也可以直接用minimap2。

# ===== 方法1：pbmm2比对（PacBio官方推荐）=====
# 建索引
pbmm2 index \                            # 建参考基因组索引
  /ref/GRCh38.fa \                       # 参考基因组FASTA
  ref.mmi                                # 输出索引文件

# 比对HiFi reads
pbmm2 align \                            # 比对命令
  ref.mmi \                              # 参考基因组索引
  ccs.bam \                              # HiFi reads（BAM格式）
  aligned.bam \                          # 输出比对后的BAM
  --sort \                               # 自动排序
  --preset HIFI \                        # 预设：HiFi模式（重要！）
  --strip \                              # 去除不必要的tag节省空间
  --unmapped \                           # 保留未比对的reads
  -j 16                                  # 16线程

# 建索引
samtools index aligned.bam               # 为BAM建索引

# pbmm2与minimap2的默认差异：
# --sort       → 自动排序
# -Y           → 使用soft clipping
# --eqx        → 使用X/=代替M（更精确的CIGAR）
# --secondary=no → 不输出次优比对

# ===== 方法2：直接用minimap2 =====
minimap2 \                               # 比对命令
  -ax map-hifi \                         # HiFi预设
  -t 16 \                                # 16线程
  /ref/GRCh38.fa \                       # 参考基因组
  ccs.fastq.gz \                         # HiFi reads（FASTQ）
  | samtools sort -@ 8 -o aligned.bam    # 排序输出
samtools index aligned.bam               # 建索引

# 比对质量评估
samtools flagstat aligned.bam            # 查看比对率统计
# HiFi reads比对率通常 > 98.9%

第四步：变异检测¶

白话解释： 找出样本和参考基因组的差异。HiFi数据的高准确率使得变异检测非常可靠。

# ===== SNP/Indel检测（DeepVariant）=====
# DeepVariant有PacBio HiFi专用模型

# Docker运行DeepVariant
docker run \                             # Docker运行
  -v /data:/data \                       # 挂载数据目录
  google/deepvariant:latest \            # DeepVariant镜像
  /opt/deepvariant/bin/run_deepvariant \ # 运行脚本
  --model_type=PACBIO \                  # 模型类型：PacBio（重要！）
  --ref=/data/GRCh38.fa \                # 参考基因组
  --reads=/data/aligned.bam \            # 输入BAM
  --output_vcf=/data/output.vcf.gz \     # 输出VCF
  --output_gvcf=/data/output.g.vcf.gz \  # 输出gVCF
  --num_shards=16                        # 并行线程数

# ===== 结构变异检测（pbsv）=====
# pbsv是PacBio官方SV检测工具

# 步骤1：发现SV信号
pbsv discover \                          # SV信号发现
  aligned.bam \                          # 输入BAM
  sv_signatures.svsig.gz                 # 输出SV信号文件

# 步骤2：调用SV
pbsv call \                              # SV调用
  /ref/GRCh38.fa \                       # 参考基因组
  sv_signatures.svsig.gz \               # SV信号文件
  sv_calls.vcf                           # 输出VCF

# ===== 单倍型定相（HiPhase）=====
hiphase \                                # 单倍型定相
  --bam aligned.bam \                    # BAM文件
  --vcf output.vcf.gz \                  # SNP VCF
  --reference /ref/GRCh38.fa \           # 参考基因组
  --output-bam phased.bam \              # 定相后的BAM
  --output-vcf phased.vcf.gz             # 定相后的VCF

第五步：基因组从头组装（hifiasm）¶

白话解释： 不依赖参考基因组，直接从HiFi reads拼出完整基因组。hifiasm可以同时组装父源和母源两个单倍型。

# 安装hifiasm
conda install -c bioconda hifiasm        # conda安装

# ===== 基本组装（仅HiFi）=====
hifiasm \                                # 组装命令
  -o sample_asm \                        # 输出前缀
  -t 32 \                                # 32线程
  ccs.fastq.gz                           # HiFi reads

# 输出文件（GFA格式）：
# sample_asm.bp.p_ctg.gfa  → 主要contigs（primary）
# sample_asm.bp.hap1.p_ctg.gfa → 单倍型1
# sample_asm.bp.hap2.p_ctg.gfa → 单倍型2

# GFA转FASTA
awk '/^S/{print ">"$2"\n"$3}' \          # 从GFA提取序列
  sample_asm.bp.p_ctg.gfa \              # 输入GFA
  > assembly.fasta                       # 输出FASTA

# ===== 三代组装（HiFi + Hi-C 联合）=====
# 加入Hi-C数据可以改善单倍型分离
hifiasm \                                # 组装命令
  -o sample_hic_asm \                    # 输出前缀
  -t 32 \                                # 32线程
  --h1 hic_R1.fastq.gz \                 # Hi-C Read1
  --h2 hic_R2.fastq.gz \                 # Hi-C Read2
  ccs.fastq.gz                           # HiFi reads

# ===== 超长模式（HiFi + ONT ultra-long）=====
hifiasm \                                # 组装命令
  -o sample_ul_asm \                     # 输出前缀
  -t 32 \                                # 32线程
  --ul ont_ultralong.fastq.gz \          # ONT超长reads
  ccs.fastq.gz                           # HiFi reads

# 组装质量评估
quast assembly.fasta \                   # QUAST评估
  -r /ref/GRCh38.fa \                    # 参考基因组
  -o quast_results/ \                    # 输出目录
  --large                                # 大基因组模式

# BUSCO完整性评估
busco \                                  # BUSCO评估
  -i assembly.fasta \                    # 输入组装
  -l mammalia_odb10 \                    # 哺乳动物数据集
  -o busco_results \                     # 输出目录
  -m genome                              # 基因组模式

常见报错与解决¶

报错	原因	解决方案
`ccs: too few passes`	DNA片段太长导致无法多圈	调低--min-passes到2，或改用CLR模式
`pbmm2: preset not recognized`	版本太旧不支持HIFI预设	更新pbmm2到v1.13.0+
`hifiasm: killed`	内存不足	人类基因组组装需要≥64GB内存
`DeepVariant: low quality calls`	覆盖度不足	HiFi建议≥30x覆盖度做变异检测
`pbsv: no SVs found`	SV过滤太严格	降低最小SV长度阈值
`CCS yield too low`	文库片段太大，pass数太少	优化文库构建，目标片段15-20kb

速查表¶

# ========== PacBio HiFi分析速查 ==========

# 1. CCS生成（新仪器可跳过）
ccs subreads.bam ccs.bam --min-rq 0.99 --min-passes 3

# 2. 比对
pbmm2 align ref.mmi ccs.bam aligned.bam --sort --preset HIFI
# 或：minimap2 -ax map-hifi ref.fa ccs.fq | samtools sort -o aligned.bam

# 3. 变异检测
# SNP: DeepVariant --model_type=PACBIO
# SV:  pbsv discover + pbsv call
# 定相: hiphase

# 4. 组装
hifiasm -o asm -t 32 ccs.fq.gz          # 基本组装
hifiasm -o asm -t 32 --h1 hi_c_R1.fq --h2 hi_c_R2.fq ccs.fq.gz  # +Hi-C

# 5. 质量评估
quast assembly.fasta -r ref.fa           # 组装质量
busco -i assembly.fasta -l lineage -m genome  # 完整性

面试高频问题¶

Q1: PacBio HiFi和CLR有什么区别？

CLR（Continuous Long Reads）是单次通过，读长可达25-175kb但错误率8-15%。HiFi是多次环形读取后取共识，读长10-25kb但准确率≥99.9%。选择取决于需求：需要超长读长（>50kb）选CLR，需要高准确率选HiFi。目前HiFi是主流选择。

Q2: pbmm2和minimap2有什么区别？

pbmm2本质上是minimap2的包装器(wrapper)，为PacBio数据设置了优化的默认参数：启用soft clipping(-Y)、使用X/=代替M(--eqx)、禁用次优比对(--secondary=no)。直接用minimap2加-ax map-hifi预设也可以，效果类似。

Q3: hifiasm组装的输出是什么？

hifiasm输出GFA格式文件，包含三组contigs：(1) primary contigs（bp.p_ctg）——主要组装，最常用；(2) hap1.p_ctg和hap2.p_ctg——两个单倍型。加入Hi-C或trio数据后，单倍型分离质量更好。输出的GFA需要转换为FASTA才能用于后续分析。

Q4: HiFi数据做基因组组装需要多少深度？

一般推荐30-40x HiFi覆盖度。人类基因组（3Gb）需要约90-120Gb HiFi数据。如果加入ONT超长reads（10-20x）或Hi-C数据，可以大幅改善组装连续性，接近T2T（端粒到端粒）水平。

Q5: PacBio和Nanopore长读长测序怎么选？

PacBio HiFi优势：准确率极高（>99.9%），适合精确变异检测和基因组组装。Nanopore优势：读长更长（可达Mb级），设备便携（MinION），可实时测序，能直接检测碱基修饰。趋势：两者结合使用——HiFi提供高准确度骨架，ONT ultra-long提供跨越重复区域的读长。