基于PacBio HiFi数据的人类全基因组重测序变异分析流程

随着第三代测序技术，特别是PacBio HiFi（High Fidelity）测序技术的发展，我们能够获得兼具长读长和高准确度的测序数据。这为人类全基因组重测序（WGS）分析，尤其是复杂区域和结构性变异（Structural Variation, SV）的检测，带来了革命性的进步。本文旨在梳理一套完整的利用PacBio Sequel II或Revio平台产生的HiFi数据进行人类基因组变异分析的流程，详细介绍从原始数据处理、序列比对、变异检测、注释、过滤到可视化的各个环节，并涵盖所涉及的关键软件工具（如pbmm2, DeepVariant, pbsv, Annovar, SnpEff, AnnotSV, SnpSift, Slivar, IGV）的安装与使用细节。

一、引言

人类基因组蕴藏着生命的奥秘，个体间的遗传变异是疾病易感性、药物反应和表型多样性的基础。全基因组重测序（WGS）旨在全面检测个体相对于参考基因组的遗传变异，包括单核苷酸变异（SNV）、小片段插入缺失（Indel）和结构性变异（SV）。传统的短读长测序（如Illumina）在检测基因组重复区域和复杂SV方面存在局限性。PacBio SMRT（Single Molecule, Real-Time）测序技术产生的长读长能够跨越这些复杂区域，而其衍生的HiFi读长通过环化测序和一致性序列（Circular Consensus Sequencing, CCS）计算，准确度高达99.9%以上，完美结合了长读长和高精度的优势。这使得HiFi数据成为当前进行高质量WGS变异分析，特别是精确SV检测的理想选择。本流程将重点介绍如何利用PacBio HiFi数据和一系列生物信息学工具，完成从测序数据到最终变异解读的全过程。

二、分析流程概览

整个分析流程主要包括以下步骤：

1.数据生成与预处理: PacBio Sequel II/Revio 平台进行SMRT测序，通过CCS算法生成高精度HiFi reads。

2.序列比对: 使用 pbmm2 将HiFi reads比对到人类参考基因组。

3.变异检测:

• 使用 DeepVariant 检测SNVs和Indels。
• 使用 pbsv 检测结构性变异 (SVs)。

4.变异注释: 使用 Annovar, SnpEff (针对SNV/Indel) 和 AnnotSV (针对SV) 对检测到的变异进行功能和临床意义注释。

5.变异过滤: 使用 SnpSift 或 Slivar 根据质量、频率、注释信息等对变异进行过滤，筛选潜在的致病或重要变异。

6.变异可视化: 使用 IGV 对比对结果和变异位点进行可视化检查。

三、详细步骤与软件使用

1.数据生成与预处理 (HiFi Reads Generation)

构建文库后上机测序。原始数据（subreads）通过PacBio的 ccs 软件（通常集成在SMRT Link分析平台中，也可独立安装）进行处理，生成HiFi reads。这一步通常由测序服务提供商或测序平台完成，输出通常是 fastq.gz 或 bam 格式的HiFi reads文件。

CCS软件安装：推荐使用conda安装

conda create -n pacbio -c bioconda pbccs
conda activate pacbio

使用示例 (从Subreads BAM生成HiFi BAM):

ccs input.subreads.bam output.hifi.bam --min-passes 3 --min-rq 0.99 --chunk 1 24 # 示例参数，具体依需求调整
# --min-passes: 生成CCS所需的最少全长循环数
# --min-rq: 最低预测准确度 (e.g., 0.99 = QV20)
# --chunk: 并行处理设置 (可选)

# 输入: subreads.bam 文件 (包含ZMW信息)
# 输出: output.hifi.bam 或 output.hifi.fastq.gz 文件 (包含高精度HiFi reads)

2.序列比对 (Alignment)

使用pbmm2将高质量的HiFi reads精确映射到参考基因组上。pbmm2是PacBio基于minimap2优化的比对软件。

安装pbmm2: 推荐使用 conda安装。

conda install -c bioconda pbmm2 samtools # 同时安装samtools用于后续处理

使用示例:

#  索引参考基因组 (若未索引):
pbmm2 index ref.fasta ref.mmi

#  进行比对:
pbmm2 align ref.mmi input.hifi.bam output.aligned.bam --preset HiFi --sort -j $(nproc) --rg '@RG\tID:sample1\tSM:sample1\tPL:PACBIO'
# ref.mmi: 索引文件 (或直接使用 ref.fasta，pbmm2会自动处理)
# input.hifi.bam: HiFi reads (BAM 格式)
# output.aligned.bam: 输出的比对后BAM文件
# --preset HiFi: 使用针对HiFi reads优化的参数
# --sort: 对输出BAM按坐标排序
# -j: 使用的线程数
# --rg: 添加Read Group信息，对下游分析非常重要 (ID唯一标识, SM样本名, PL平台)

# 创建BAM索引:
samtools index output.aligned.bam

#输入: 参考基因组 FASTA 文件 (ref.fasta), HiFi reads (input.hifi.bam 或 fastq.gz)
#输出: 排序并索引的比对文件 (output.aligned.bam 和 output.aligned.bam.bai)

3.变异检测 (Variant Calling)

由于SNV/Indel和SV的最佳检测算法不同，通常分开进行。

（1）SNV 和 Indel 检测

DeepVariant 是Google开发的基于深度学习的变异检测工具，对PacBio HiFi数据有专门优化模型。

DeepVariant安装: 推荐使用 Docker 或 Singularity 容器，或预编译的二进制文件。conda 安装有时可用但需检查官方文档。

# 拉取镜像
docker pull google/deepvariant
# 运行 (注意路径映射和参数)
# 需要设置 BIN_VERSION, INPUT_DIR, REF, READS, OUTPUT_DIR 等环境变量
sudo docker run \
  -v "${INPUT_DIR}":"/input" \
  -v "${OUTPUT_DIR}":"/output" \
  google/deepvariant:"${BIN_VERSION}" \
  /opt/deepvariant/bin/run_deepvariant \
  --model_type=PACBIO \
  --ref="/input/${REF}" \
  --reads="/input/${READS}" \
  --output_vcf="/output/output.dv.vcf.gz" \
  --output_gvcf="/output/output.dv.g.vcf.gz" \
  --num_shards=$(nproc)

使用:

# 假设使用预编译二进制文件或conda安装的版本
run_deepvariant \
  --model_type="PACBIO" \
  --ref=ref.fasta \
  --reads=output.aligned.bam \
  --output_vcf=output.dv.vcf.gz \
  --output_gvcf=output.dv.g.vcf.gz \
  --num_shards=$(nproc) # 使用的线程数
# --model_type="PACBIO": 指定使用PacBio模型
# --ref: 参考基因组
# --reads: 比对后的BAM文件
# --output_vcf: 输出VCF文件 (包含检测到的变异)
# --output_gvcf: 输出gVCF文件 (包含所有位点信息，用于联合分析)

#输入: 参考基因组 FASTA, 排序并索引的 BAM 文件 (output.aligned.bam)
#输出: VCF 文件 (output.dv.vcf.gz), gVCF 文件 (output.dv.g.vcf.gz)

（2）结构性变异 (SV) 检测

pbsv安装: 推荐使用 conda。

conda install -c bioconda pbsv

使用:pbsv 通常分两步：

#  发现 SV 信号 (Discover signatures): 
pbsv discover output.aligned.bam output.svsig.gz
# output.aligned.bam: 比对后的BAM文件
# output.svsig.gz: 存储检测到的SV信号的文件

#  调用 SV (Call variants):
pbsv call ref.fasta output.svsig.gz output.pbsv.vcf --tandem-repeats ref.trf.bed
# ref.fasta: 参考基因组
# output.svsig.gz: 上一步生成的信号文件
# output.pbsv.vcf: 输出的SV VCF文件
# --tandem-repeats: (可选但推荐) 提供串联重复注释文件 (如TRF生成的BED)，提高精度

#输入: 参考基因组 FASTA, 排序并索引的 BAM 文件, (可选) 串联重复注释文件
#输出: SV VCF 文件 (output.pbsv.vcf)

4.变异注释 (Variant Annotation)

为检测到的变异添加生物学信息，如所在基因、对外显子/内含子的影响、氨基酸变化、在人群数据库中的频率、临床数据库记录等。

软件:

• Annovar: 功能强大，数据库全面，但需要注册下载和手动更新数据库。
• AnnotSV: 专门针对SV进行注释，整合多种资源。

安装:

• Annovar: 访问官网注册下载，按说明安装和下载数据库。
• AnnotSV: conda install -c bioconda annotsv。首次运行可能需要下载注释数据库。

使用示例:

Annovar (注释 SNV/Indel VCF):

# 需要预先下载好数据库到 humandb/ 目录下
table_annovar.pl output.dv.vcf.gz humandb/ -buildver hg38 \
  -out annotated_annovar \
  -remove \
  -protocol refGene,cytoBand,exac03,avsnp150,dbnsfp41a,clinvar_20230325 \
  -operation g,r,f,f,f,f \
  -nastring . -vcfinput \
  --thread 4 # 使用线程数
# 输出: annotated_annovar.hg38_multianno.vcf (或 .txt)

AnnotSV (注释 SV VCF):

AnnotSV -SVinputFile output.pbsv.vcf \
  -outputFile annotated.annotsv.tsv \
  -genomeBuild GRCh38 \
  -annotationMode FULL \
  -svtBEDcol 3 # 根据pbsv VCF格式指定类型在INFO字段的位置
# 输出: annotated.annotsv.tsv (包含详细注释的TSV文件)
# 也可输出VCF格式，需调整参数

#输入: VCF 文件 (output.dv.vcf.gz, output.pbsv.vcf), 注释数据库
#输出: 包含注释信息的 VCF 或 TSV 文件

5.变异过滤 (Variant Filtering)

根据变异质量、测序深度、等位基因频率、注释结果（如是否位于编码区、预测的致病性、人群频率等）筛选出更可靠或更具潜在生物学意义的变异。

软件:

• SnpSift: SnpEff 套件的一部分，使用类似SQL的表达式进行过滤。
• Slivar: 灵活的VCF过滤工具，特别适用于家系分析，但也可单样本过滤，支持复杂的Javascript表达式。

安装:

• SnpSift: 通常随 SnpEff 一起安装。
• Slivar: conda install -c bioconda slivar 或下载二进制文件。

使用示例:

SnpSift (过滤注释后的 SNV/Indel VCF):

java -jar SnpSift.jar filter \
  "( QUAL > 30 ) & ( DP > 10 ) & ( AF[0] >= 0.01 ) & ( exists ANN[*] ) & ( ANN[?('EFF[*].IMPACT == \'HIGH\' || EFF[*].IMPACT == \'MODERATE\')].exists() )" \
  annotated.snpeff.vcf > filtered.snpsift.vcf
# 过滤条件示例:
# QUAL > 30: 质量值大于30
# DP > 10: 深度大于10
# AF[0] >= 0.01: 第一个等位基因频率大于等于0.01
# exists ANN[*]: 必须有SnpEff注释
# ANN[?('EFF[*].IMPACT == \'HIGH\' || EFF[*].IMPACT == \'MODERATE\')].exists(): 变异影响为HIGH或MODERATE

Slivar (过滤注释后的 SNV/Indel VCF):

# 1. 创建一个简单的 pedigree 文件 (即使是单样本)
# family_id sample_id paternal_id maternal_id sex phenotype
# sample1   sample1   0           0           0   0
echo -e "family_id\tsample_id\tpaternal_id\tmaternal_id\tsex\tphenotype\nsample1\tsample1\t0\t0\t0\t0" > sample.ped

# 2. 使用表达式过滤 (假设VCF已包含如gnomAD频率注释)
slivar expr --vcf annotated_annovar.hg38_multianno.vcf \
  --ped sample.ped \
  --pass-only \
  --info 'INFO.gnomad_genome_AF < 0.01 || INFO.gnomad_genome_AF == "."' \
  --filter "QUAL > 30 && INFO.DP > 10" \
  -o filtered.slivar.vcf
# --pass-only: 只保留通过所有过滤器(包括VCF本身的FILTER列)的记录
# --info: 基于INFO字段的过滤条件
# --filter: 基于QUAL, FILTER, INFO字段的通用过滤条件

#输入: 注释后的 VCF 文件
#输出: 过滤后的 VCF 文件 (filtered.*.vcf)

6.变异可视化 (Variant Visualization)

在基因组浏览器中直观地检查比对情况和变异位点，确认可疑变异，理解其基因组上下文。

软件:IGV (Integrative Genomics Viewer, Broad Institute)

安装: 从 IGV 官网下载适用于操作系统的版本 。

使用:

（1）启动 IGV。

（2）加载参考基因组 (File > Load Genome... 或选择预置的，如 hg38)。

（3）加载比对文件 (File > Load from File... 选择 output.aligned.bam 文件，确保 .bai 索引文件在同一目录)。

（4）加载VCF文件 (File > Load from File... 选择 filtered.*.vcf 或原始的 *.vcf.gz 文件，确保 .tbi 索引文件在同一目录)。

（5）在顶部输入基因名或基因组坐标进行导航。

（6）检查特定变异位点：查看reads的比对情况（颜色表示链方向、错配等），覆盖深度，VCF记录的详细信息（通过点击VCF轨道上的变异查看），以及注释信息。对于SV，可以观察到支持SV的reads特征，如跨越断裂点的reads，异常插入大小的read pairs (虽然HiFi主要是单分子长reads，但比对模式能反映SV)。

输入: 参考基因组, BAM+BAI 文件, VCF+TBI 文件。

输出: 可视化检查结果，截图用于报告或进一步分析。

四、结论与讨论

利用PacBio HiFi数据进行人类WGS变异分析，结合pbmm2, DeepVariant, pbsv, Annovar/SnpEff/AnnotSV, SnpSift/Slivar 和 IGV 等工具，可以构建一个强大而全面的分析流程。该流程能够高效、准确地检测SNVs、Indels，并在检测复杂基因组区域的结构性变异方面展现出显著优势。

关键优势:

• HiFi reads的长读长和高精度特性显著提高了比对准确性，尤其是在重复和低复杂度区域。
• 极大提升了对各类结构性变异（缺失、插入、重复、易位、倒位）的检测灵敏度和断点解析精度。
• 有助于单倍型分型（phasing），更好地理解等位基因特异性表达或复合杂合变异。

注意事项:

• 计算资源: WGS分析，特别是长读长数据分析，需要大量的计算资源（CPU、内存）和存储空间。
• 参数调优: 各个软件的参数设置可能需要根据具体数据质量和研究目标进行调整。
• 数据库更新: 注释数据库需要定期更新以获取最新的功能、频率和临床信息。
• 结果验证: 对于关键发现，可能需要通过其他方法（如PCR测序、其他平台数据）或在更大队列中进行验证。