历史数据的整合分析往往涉及不同时期使用不同参考基因组版本产生的数据。例如，某些较早的研究可能基于 hg19/GRCh37 进行分析，而新的研究多采用 hg38/GRCh38。为了进行比较分析，需要将这些数据统一到同一个版本。

坐标转换是通过 chain file 来实现的，该文件记录了两个基因组版本之间的对应关系，包含了基因组重排、插入和删除等信息。转换工具会根据这些信息将原始坐标映射到目标版本的相应位置。

正文

本文将介绍如何使用两款主流工具 CrossMap 和 GATK LiftoverVcf 以及最近的一个新工具 BCFtools liftover 进行 VCF 文件的坐标转换。在说明之前，首先简要介绍三个软件的差别：

  • CrossMap 除了 VCF 的坐标映射外，也支持很多其他文件格式的坐标转换。对于 VCF 文件而言，CrossMap 并不要求其具有完备的 Meta-information 和 Header line 部分，只需要突变数据行符合格式即可。
  • GATK LiftoverVcf 则主要支持对 VCF 文件的坐标进行转换。其要求 VCF 文件具有完整的 Meta-information 和 Header line 部分。相对 CrossMap 而言对 VCF 文件格式的要求更加严格。
  • BCFtools liftover 在前两者的基础上，对 InDel 的转换逻辑进行了完善，从而支持更准确的 InDel 坐标转换，并且相对而言功能更加完善齐全（例如多等位基因的支持、INFO 字段的更新等）& 效率更高。其对 Meta-information 的要求更加严格，需要包含原始坐标的 contig 长度信息。
各种坐标转换工具的内存使用和时间开销对比

VCF 文件格式详解：https://samtools.github.io/hts-specs/VCFv4.2.pdf

Meta-information & Header line 示例可见以下文件开头以 # 起始的行：

https://github.com/vcflib/vcflib/blob/master/samples/sample.vcf

CrossMap

github 链接：https://github.com/liguowang/CrossMap

官方文档链接：https://crossmap.readthedocs.io/en/latest/

使用 pip 安装 CrossMap：

1
2
3
pip3 install git+https://github.com/liguowang/CrossMap.git # 从 github 上安装，相对而言可能更新一些
或者
pip3 install CrossMap # 从 PyPI 安装，按照版本发布情况安装

需准备的输入文件：

  • 指定了转换方向的 chain 文件，可以从 UCSC 网站下载（所有物种的下载信息详细可见 UCSC download 中的 LiftOver files）。以人类基因组为例，hg19 到 hg38 的 chain file 下载链接: hg19ToHg38.over.chain.gz。
  • 目标参考基因组文件，例如转坐标是 hg19 到 hg38 的话，则准备 hg38 的参考基因组文件。

运行命令示例（链文件和参考基因组文件按需求更改为自己的路径）：

1
CrossMap vcf [chain file] [input vcf] [target genome] [output vcf]

各框注部分需填写的文件路径：

  • [chain file]：填入使用的链文件路径，例如 hg19ToHg38.over.chain.gz。
  • [target genome]：填入转换的目标基因组文件，例如 hg38.fa。
  • [input vcf]：填入需要进行转换的 VCF 文件路径
  • [output vcf]：填入输出的转换后的 VCF 文件路径。

可选参数如下：

  • --chromid：指定输出的染色体 id 格式，详细可见 CrossMap vcf -h 说明。
  • --no-comp-alleles：指定该参数后，CrossMap 会取消对 ref & alt 的检查。默认不指定该参数时，如果转换后 ref 和 alt 一致则会被认为转换失败。
  • --compress：指定该参数后，CrossMap 会调用 gzip 对输出的 VCF 文件进行压缩。

未转换成功的变异信息可见 [output vcf].unmap。

GATK LiftoverVcf

官方文档链接：https://gatk.broadinstitute.org/hc/en-us/articles/360036831351-LiftoverVcf-Picard

使用 conda 安装 GATK：

1
conda install -c bioconda gatk4

需准备的输入文件与 CrossMap 相同（可见上文），以下是一个运行示例：

1
2
3
4
5
6
gatk LiftoverVcf \
  -I [input vcf] \
  -O [output vcf] \
  -CHAIN [chain file] \
  -R [target genome]  \
  --REJECT [unmap vcf]

各框注部分需填写的文件路径：

  • [input vcf]：填入需要进行转换的 VCF 文件路径
  • [output vcf]：填入输出的转换后的 VCF 文件路径。
  • [unmap vcf]：指定转换失败的变异输出文件。
  • [chain file]：填入使用的链文件路径，例如 hg19ToHg38.over.chain.gz。
  • [target genome]：填入转换的目标基因组文件，例如 hg38.fa。

需注意 GATK 会识别文件尾缀以判断输出的文件类型，因此相关的 VCF 文件名应当以 .vcf 结尾。

一些可能比较常用的参数（详细信息请使用 gatk LiftoverVcf -h 查看）：

  • --ALLOW_MISSING_FIELDS_IN_HEADER：允许 header 中缺失 INFO & FORMAT 字段的说明，但仍需其他的部分（例如 header line 等）。
  • --WRITE_ORIGINAL_ALLELES & --WRITE_ORIGINAL_POSITION：输出原始的等位基因信息 & 输出原始的坐标信息。

BCFtools liftover

该工具是 BCFtools 的插件，因此需要先安装 BCFtools 后再使用。

github 链接：https://github.com/freeseek/score?tab=readme-ov-file#liftover-vcfs

安装方法：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
# 首先需要安装 bcftools，可以自行编译或者使用 conda 安装
conda install -c bioconda bcftools

# 自行选择一个合适的路径，开始下载插件文件 (Linux)
# 插件也可以自行编译，有需求请在上方 github 链接中自行查阅
wget https://software.broadinstitute.org/software/score/score_1.20-20240927.zip
mkdir -p bcftools_plugin
unzip -d bcftools_plugin score_1.20-20240927.zip

# 将插件的路径设置为 BCFtools 插件的环境路径
export BCFTOOLS_PLUGINS=$(pwd .)/bcftools_plugin
# 或者保存到 .bashrc 文件中保证以后能直接调用
echo "export BCFTOOLS_PLUGINS=$(pwd .)/bcftools_plugin" >> ~/.bashrc
source ~/.bashrc

查看是否安装成功：

1
bcftools +liftover -h

需准备的输入文件与 CrossMap 相同，在此基础上还需要准备源基因组版本的基因组序列文件，以下是一个运行示例：

1
2
3
4
5
6
bcftools +liftover -Ou [input vcf] -- \
    -s [source genome] \
    -f [target genome] \
    -c [chain file] \
    --reject [unmap vcf] --write-reject | \
    bcftools sort -Ou -o [output vcf]

各框注部分需填写的文件路径：

  • [input vcf]：填入需要进行转换的 VCF 文件路径
  • [output vcf]：填入输出的转换后的 VCF 文件路径。
  • [unmap vcf]：指定转换失败的变异输出文件。
  • [chain file]：填入使用的链文件路径，例如 hg19ToHg38.over.chain.gz。
  • [source genome]：填入转换的源基因组文件，例如 hg19.fa。
  • [target genome]：填入转换的目标基因组文件，例如 hg38.fa。

更多可选项的说明详细可见 bcftools +liftover -h。以上命令行在输出后对 VCF 进行了排序并以未压缩格式输出（-Ou），可根据自身需求调整以上命令的工作流。

后记

关于三个工具，就我个人的使用经验而言（主要是一些 SNV），三者的结果相差并不大，所以目前使用较多的是更方便的 CrossMap（因为它不对 header 做要求）。如果在你的应用场景中涉及到较多的 InDel，可以考虑使用专门为 InDel 进行过优化的 BCFtools liftvcf 进行。以下是 BCFtools liftvcf 文章中的一个表格，展示了各个工具的功能和限制，可按照自己的需求进行挑选：

Tool BCFtools liftvcf GATK LiftoverVcf CrossMap
Version tested 2023-12-06 3.1.1 0.6.6
GitHub username freeseek broadinstitute liguowang
GitHub repository score picard CrossMap
License MIT MIT GPLv3
Main features
Can reverse-complements alleles Yes Yes Yes
Handles multi-allelic records Yes No No
Can swap SNV alleles Yes Bi-allelic only No
Can swap indel alleles Yes No No
Can add new reference allele Yes No No
Can recover SNVs at chain gaps Yes No No
File input/output options
Sort records after liftover No Yes No
Left-aligns indels after liftover Yes Yes No
Can record the original position Yes Yes No
Flexible with contig names Yes No Yes
Loads full reference in memory No Yes No
Can input VCF as a file stream Yes No Yes
Can output VCF as a file stream Yes No No
Can input/output binary VCFs Yes No No
Variants annotations updates
Updates INFO/END field Yes No No
Updates Number=G/R fields Yes PL and AD No
Updates AC-like fields Yes No No
Updates AF-like fields Yes Yes No
Updates GWAS-VCF fields Yes No No
Source: https://academic.oup.com/view-large/438467641

但也正如该文章中结尾的部分所说，不管哪个工具能做得更好，liftover 终究是一个有损过程，无论什么情况下，重新比对进行 call variant 总会产生比 liftover 更好的结果。当然有些情况下我们也不可避免地只能选择 liftover 这一途径（比如说原始测序数据无法访问或需申请权限等）。

相关文章链接