跳转至

图基因组比对分析

一句话概述

利用vg、GraphAligner和minigraph等图基因组工具，将测序reads比对到包含群体变异信息的图参考基因组上，克服传统线性参考基因组的参考偏好性(reference bias)问题，实现更准确的变异检测和基因型分析，是下一代基因组学分析的前沿方法。

核心知识点表格

知识点关键内容常用工具
图基因组概念多路径DAG表示群体遗传多样性-
图构建从VCF/多序列比对构建变异图vg construct, minigraph
短读比对短reads比对到图参考vg giraffe, GraphAligner
长读比对长reads比对到图GraphAligner, minigraph
变异检测从图比对中genotypingvg call, vg pack
参考偏好性线性参考导致非参考等位基因mapping偏差-
泛基因组图人类/作物泛基因组参考HPRC, pggb
结构变异图基因组的SV检测优势vg, pangenie

各步骤详解

第一步：理解图基因组原理

白话解释： 传统参考基因组是一条线性序列，代表一个"平均"个体。但人群中存在大量变异——如果你的基因组在某个位置和参考不同，reads在那里就很难正确比对（参考偏好性）。图基因组把所有已知变异都编码进去，形成一个有分支和汇聚的"路线图"。这样reads可以选择最匹配自己的路径，无论它是否与"主参考"一致。

技术细节： - 图基因组是一个有向无环图(DAG)：节点是序列片段，边表示可能的连接 - 每个单倍型(haplotype)对应图中的一条路径 - 参考偏好性：线性参考对非参考等位基因的比对率降低~1-3% - 这个偏差在SV(结构变异)附近最严重 - 人类泛基因组参考(HPRC) Release 1(2023)包含47个个体的94个单倍型基因组；Release 2(2025年5月)已扩展至230+个样本，覆盖更广泛的人类遗传多样性

第二步：构建图基因组

白话解释： 可以从两种来源构建图基因组：1）线性参考+VCF变异文件；2）多个完整基因组的全基因组比对。前者简单但只包含已知变异；后者全面但计算量大。

代码示例：

# 安装vg
conda install -c bioconda vg

# === 方法1: 从VCF构建图 ===
# 需要: 参考基因组FASTA + 变异VCF

# 构建图（单条染色体）
vg construct \
    -r reference.fa \
    -v variants.vcf.gz \
    -R chr1 \                    # 指定染色体
    -C \                         # 使用alt路径
    -a \                         # 保留alt allele信息
    > chr1.vg

# 查看图统计信息
vg stats -z chr1.vg
# 输出: 节点数、边数、路径数、序列总长

# === 方法2: 使用minigraph从多基因组构建 ===
conda install -c bioconda minigraph

# minigraph: 增量式图构建
# 第一个FASTA作为骨架，后续FASTA的变异增量添加
minigraph -cxggs \
    reference.fa \
    sample1.fa \
    sample2.fa \
    sample3.fa \
    > pangenome.gfa

# 转换GFA为vg格式
vg convert -g pangenome.gfa -v > pangenome.vg

# === 方法3: 使用pggb从头构建泛基因组图 ===
# pggb: 全基因组比对→图构建
conda install -c bioconda pggb

pggb -i all_genomes.fa.gz \
    -o pggb_output/ \
    -n 47 \                      # 单倍型数
    -t 32 \
    -p 98 \                      # 序列相似性阈值
    -s 100000 \                  # segment大小
    -k 79 \                      # 最小匹配长度（min-match-length）
    --multiqc

# === 方法4: 使用HPRC预建图 ===
# 人类泛基因组参考联盟已发布预建图
# 下载地址: https://github.com/human-pangenomics/hpp_pangenome_resources
wget https://s3-us-west-2.amazonaws.com/human-pangenomics/pangenomes/freeze/freeze1/minigraph-cactus/hprc-v1.1-mc-grch38/hprc-v1.1-mc-grch38.gbz

第三步：图基因组索引

白话解释： 就像BWA需要对线性参考建索引一样，图比对工具也需要对图基因组建索引。vg giraffe使用GBWT (Graph BWT) 索引来高效搜索图中的路径，这是实现快速比对的关键。

代码示例：

# vg giraffe需要的索引（推荐方式）

# 步骤1: 从VG创建GBZ (Graph BZ) 格式
# 先创建GBWT（记录已知单倍型路径）
vg gbwt -x chr1.vg -o chr1.gbwt \
    --vcf-input variants.vcf.gz \
    --preset 1000gp                # 使用1000基因组preset

# 步骤2: 创建GBZ文件（vg giraffe的主索引）
vg gbz -x chr1.vg -g chr1.gbwt \
    --gbz-format -o chr1.gbz

# 步骤3: 创建minimizer索引
vg minimizer -d chr1.dist \      # 距离索引
    -o chr1.min \                 # minimizer索引
    chr1.gbz

# 或使用自动化索引构建
vg autoindex \
    --workflow giraffe \
    -r reference.fa \
    -v variants.vcf.gz \
    -p index_prefix \
    -t 16

# 输出文件:
# index_prefix.giraffe.gbz
# index_prefix.min
# index_prefix.dist

# 对于minigraph-cactus图（如HPRC）
# 已经提供预建索引，直接使用即可

第四步：短读长比对（vg giraffe）

白话解释： vg giraffe是当前最快的图比对工具，速度接近BWA-MEM。它先用minimizer快速定位reads可能来自图中的哪些区域，然后在局部图中精确比对。输出GAM/GAF格式（图比对格式），也可转为标准BAM。

代码示例：

# 使用vg giraffe比对短reads
vg giraffe \
    -Z index_prefix.giraffe.gbz \
    -m index_prefix.min \
    -d index_prefix.dist \
    -f reads_R1.fastq.gz \
    -f reads_R2.fastq.gz \
    -o gaf \                      # 输出GAF格式
    -t 16 \
    -p \                          # 进度显示
    > aligned.gaf

# 转换GAF到BAM（用于传统工具兼容）
vg surject \
    -x index_prefix.giraffe.gbz \
    -t 16 \
    -b \                          # 输出BAM
    -N sample_name \
    -R "ID:sample LB:lib PL:illumina SM:sample" \
    aligned.gaf | \
    samtools sort -@ 4 -o aligned_surjected.bam

samtools index aligned_surjected.bam

# 比较图比对和线性比对的统计
# 图比对通常有:
# - 更高的mapping率(+0.5-2%)
# - 更少的mapping质量为0的reads
# - 在SV区域更准确的比对

# 查看比对统计
vg stats -a aligned.gaf

第五步：长读长比对（GraphAligner）

白话解释： 长reads（PacBio/Nanopore）因为长度长、错误率高，比对策略与短reads不同。GraphAligner专门处理长reads到图基因组的比对，使用seed-chain-extend策略。

代码示例：

# 安装GraphAligner
conda install -c bioconda graphaligner

# 长reads比对到GFA图
GraphAligner \
    -g pangenome.gfa \
    -f long_reads.fastq.gz \
    -a aligned_long.gaf \
    -x vg \                       # 使用vg compatible模式
    -t 16 \
    --seeds-minimizer-count 5 \   # minimizer seed参数
    --seeds-minimizer-length 19

# 使用minigraph比对长reads
minigraph -cxlr \                 # lr = long reads模式
    pangenome.gfa \
    long_reads.fastq.gz \
    -t 16 \
    > aligned_minigraph.gaf

# 解析GAF格式
# GAF格式类似PAF，但路径信息包含图中的节点遍历
# 列: read_name, read_len, query_start, query_end, strand,
#     path, path_len, path_start, path_end, matches, block_len, mapq
head aligned_long.gaf

第六步：变异检测与基因型分析

白话解释： 从图比对结果中检测变异和确定基因型。vg可以直接从图比对（GAM/GAF）中进行genotyping——因为图中已经包含了已知变异，所以只需要确定reads更支持哪个等位基因路径。这比从头发现变异更准确。

代码示例：

# 方法1: vg pack + vg call（从图比对genotyping）

# 计算每个节点/边的reads支持度
vg pack \
    -x index_prefix.giraffe.gbz \
    -a aligned.gaf \
    -o aligned.pack \
    -t 16

# 基因型判定
vg call \
    index_prefix.giraffe.gbz \
    -k aligned.pack \
    -t 16 \
    -s sample_name \
    > genotyped_variants.vcf

# 方法2: DeepVariant + 图比对BAM
# 先surject到线性参考
vg surject -x index_prefix.giraffe.gbz -b aligned.gaf | \
    samtools sort -o surjected.bam
samtools index surjected.bam

# 用DeepVariant调用变异（surjected BAM）
docker run google/deepvariant:latest /opt/deepvariant/bin/run_deepvariant \
    --model_type=WGS \
    --ref=reference.fa \
    --reads=surjected.bam \
    --output_vcf=deepvariant_output.vcf.gz \
    --num_shards=16

# 方法3: pangenie（从图基因组genotyping短reads）
conda install -c bioconda pangenie

pangenie \
    -i reads_R1.fastq.gz \
    -r reference.fa \
    -v panel_variants.vcf \       # 包含SV的VCF
    -o pangenie_output \
    -t 16 \
    -j 16

# pangenie特别适合结构变异genotyping

# 比较图比对vs线性比对的变异检测改进
# 通常在以下方面有改进:
# - SV检测灵敏度提高10-30%
# - 杂合位点genotyping准确率提高
# - 在高度多态区域（如MHC/HLA）改进最大

# 评估改进
bcftools stats genotyped_graph.vcf genotyped_linear.vcf > comparison_stats.txt

第七步：图基因组优势分析

白话解释： 图基因组相比线性参考的核心优势是消除参考偏好性。在某些区域（如MHC、着丝粒附近的高度多态区、结构变异断点），改进尤其明显。

代码示例：

# 量化参考偏好性改善
import pandas as pd
import matplotlib.pyplot as plt

# 比较线性比对和图比对在已知杂合位点的等位基因平衡
linear_balance = pd.read_csv("linear_allele_balance.tsv", sep="\t")
graph_balance = pd.read_csv("graph_allele_balance.tsv", sep="\t")

# 理想的杂合位点: REF reads ≈ ALT reads (比值0.5)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

ax1.hist(linear_balance["alt_fraction"], bins=50, alpha=0.7)
ax1.axvline(x=0.5, color="red", linestyle="--")
ax1.set_title("Linear reference alignment")
ax1.set_xlabel("ALT allele fraction at het sites")

ax2.hist(graph_balance["alt_fraction"], bins=50, alpha=0.7)
ax2.axvline(x=0.5, color="red", linestyle="--")
ax2.set_title("Graph genome alignment")
ax2.set_xlabel("ALT allele fraction at het sites")

plt.tight_layout()
plt.savefig("reference_bias_comparison.pdf")

# 图比对的ALT fraction应更接近0.5（更对称）
# 线性比对通常向REF偏移（ALT fraction < 0.5）


# 比较mapping率
echo "=== Linear alignment ==="
samtools flagstat linear_aligned.bam | head -5

echo "=== Graph alignment ==="
samtools flagstat graph_surjected.bam | head -5

# 在困难区域(MHC)的比较
echo "=== MHC region mapping rate ==="
samtools view -c -F 4 linear_aligned.bam chr6:28000000-34000000
samtools view -c -F 4 graph_surjected.bam chr6:28000000-34000000

实战命令

#!/bin/bash
# === 图基因组分析完整流程 ===
set -e
THREADS=32
SAMPLE="HG002"

# 1. 下载HPRC图基因组
wget -P ref/ https://s3-us-west-2.amazonaws.com/human-pangenomics/pangenomes/freeze/freeze1/minigraph-cactus/hprc-v1.1-mc-grch38/hprc-v1.1-mc-grch38.gbz
wget -P ref/ https://s3-us-west-2.amazonaws.com/human-pangenomics/pangenomes/freeze/freeze1/minigraph-cactus/hprc-v1.1-mc-grch38/hprc-v1.1-mc-grch38.min
wget -P ref/ https://s3-us-west-2.amazonaws.com/human-pangenomics/pangenomes/freeze/freeze1/minigraph-cactus/hprc-v1.1-mc-grch38/hprc-v1.1-mc-grch38.dist

# 2. 图比对
vg giraffe -Z ref/hprc-v1.1-mc-grch38.gbz \
    -m ref/hprc-v1.1-mc-grch38.min \
    -d ref/hprc-v1.1-mc-grch38.dist \
    -f ${SAMPLE}_R1.fq.gz -f ${SAMPLE}_R2.fq.gz \
    -o gaf -t $THREADS > ${SAMPLE}.gaf

# 3. Surject到线性参考
vg surject -x ref/hprc-v1.1-mc-grch38.gbz -b -t $THREADS \
    -N $SAMPLE -R "ID:${SAMPLE} SM:${SAMPLE} PL:illumina" \
    ${SAMPLE}.gaf | samtools sort -@ 8 -o ${SAMPLE}_surjected.bam
samtools index ${SAMPLE}_surjected.bam

# 4. 变异检测
vg pack -x ref/hprc-v1.1-mc-grch38.gbz -a ${SAMPLE}.gaf -o ${SAMPLE}.pack -t $THREADS
vg call ref/hprc-v1.1-mc-grch38.gbz -k ${SAMPLE}.pack -s $SAMPLE -t $THREADS > ${SAMPLE}_variants.vcf

echo "Pipeline complete!"

面试常问点

Q1: 图基因组相比线性参考的核心优势是什么？ A: 1）消除参考偏好性：非参考等位基因的reads不再因为与参考不匹配而被错误比对或丢弃；2）结构变异检测改进：图中显式编码SV，reads可以直接比对到SV路径；3）难区域改进：MHC、着丝粒等高度多态区域显著改善；4）更公平的群体代表性：不再以单一个体为参考。

Q2: vg giraffe和BWA-MEM的速度差异？ A: vg giraffe的速度已接近BWA-MEM——在映射到包含数千个单倍型的泛基因组图时，速度与BWA-MEM映射到单一线性参考大致相当（处理30x WGS约需2-6小时）。这是相对于早期vg map（比BWA-MEM慢一个数量级以上）的巨大进步。NVIDIA Parabricks还提供了GPU加速版vg giraffe，可进一步提升速度。

Q3: 图基因组在哪些区域改进最大？ A: 1）MHC/HLA区域（人类最多态的区域）：genotyping准确率提升10-30%；2）已知SV断点附近：mapping率和genotyping显著改善；3）高重复区域：图中冗余路径帮助区分重复拷贝；4）群体特异的common变异位点。对高度保守区域改进很小。

Q4: HPRC泛基因组参考是什么？ A: HPRC (Human Pangenome Reference Consortium) 发布的人类泛基因组参考。Release 1（2023年）包含47个个体的94个单倍型解析基因组；Release 2（2025年5月）扩展至200+个个体（400+个单倍型），覆盖更广泛的人类遗传多样性。用minigraph-cactus方法构建，预计将逐步替代GRCh38作为标准参考。

Q5: 图比对的BAM文件能直接用GATK等传统工具分析吗？ A: 可以。通过vg surject将图比对投影(surject)到线性参考坐标系，生成标准BAM文件。这个BAM可以用GATK/DeepVariant等传统工具分析。但这样会丢失一部分图比对的优势（尤其是SV区域）。理想方案是使用vg原生的变异检测流程。

易错点

  1. 索引版本不匹配：GBZ、minimizer和距离索引必须从同一个图构建，不能混用不同版本。

  2. 内存不足：人类全基因组图索引加载需要约30-50GB内存。确保机器有足够内存。

  3. GAF格式不兼容：GAF（图比对格式）不能直接用samtools等传统工具处理，需要先surject到BAM。

  4. VCF坐标系统差异：从图比对得到的VCF可能使用图内部坐标而非传统线性坐标，需要注意转换。

  5. 图中缺失变异：如果个体携带的变异未被包含在图中，图比对的优势就体现不出来。图的全面性很重要。

  6. 过度复杂的图：包含过多罕见变异会使图过于复杂，反而降低比对速度和准确性。HPRC使用了合理的变异频率过滤。

补充知识

图基因组工具生态

工具功能适用
vg全功能图基因组工具套件短reads/长reads/变异检测
vg giraffe快速短reads图比对WGS/WES短读长
minigraph增量式图构建+长reads比对长reads/基因组比较
GraphAligner长reads图比对PacBio/Nanopore
pggb全基因组图构建从头构建泛基因组
pangenie图genotypingSV genotyping
odgi图操作和可视化图编辑/统计/可视化

图基因组数据格式

GFA (Graphical Fragment Assembly):
  - 文本格式，定义图的节点(S)和边(L)
  - GFA1: 基本格式，包含S(segment)、L(link)、P(path)行；v1.1新增W(walk)行，v1.2新增J(jump)行
  - GFA2: GFA1的超集，用E(edge)行替代L/C行，用O/U行替代P行，新增F(fragment)和G(gap)行

GBZ (Graph BWT ZIP):
  - vg giraffe的压缩索引格式
  - 包含图+GBWT单倍型索引

GAF (Graph Alignment Format):
  - 类似PAF但包含图路径信息
  - 路径用节点遍历表示: >1>2>4>5

GAM (Graph Alignment/Map):
  - vg的二进制比对格式
  - protobuf编码，紧凑但不可读

从线性到图基因组的过渡路线图

当前(2024-2026):
  GRCh38/T2T仍是主要参考
  图基因组用于特殊区域(MHC/SV)
  HPRC图参考作为补充

近期(2026-2028):
  HPRC图参考逐步成为标准
  临床基因组学开始采用图参考
  工具链成熟化

远期(2028+):
  图参考成为默认选择
  每个人的基因组直接映射到图
  个性化图参考(加入家系信息)

📚 相关文章推荐

46. 泛基因组分析（Pan-genome Analysis）
38% 相关
全基因组测序WGS变异检测
24% 相关
全外显子组测序（WES）分析全流程
22% 相关
BWA-MEM2 比对加速
19% 相关
泛基因组分析
18% 相关