基因组版本升级实战：bed与vcf文件坐标转换全攻略

1. 为什么需要基因组版本升级？

做生物信息分析的朋友们应该都遇到过这样的场景：好不容易在hg19参考基因组上做完变异检测，突然发现最新的数据库和工具都转向hg38了。这时候就面临一个头疼的问题——如何把已有的bed和vcf文件从hg19转换到hg38？我当年第一次遇到这种情况时，整整折腾了两天才搞定。今天就把这些年积累的实战经验分享给大家，手把手教你用UCSC的liftover和picard工具完成坐标转换。

基因组版本升级的核心难点在于不同版本间染色体坐标系统的变化。举个生活中的例子，就像城市道路重新规划后，原来的"中山路123号"可能变成了"解放路56号"。基因组版本更新时，由于序列修正、填补gap等原因，基因和变异的位置坐标也会发生变化。根据我的经验，hg19到hg38的转换中，约85%的坐标变化在±100bp范围内，但也有约5%的变异会因序列重构而无法直接映射。

2. 准备工作：获取必要工具和文件

2.1 下载chain文件

chain文件是坐标转换的"翻译字典"，记录了不同基因组版本间的对应关系。UCSC提供了各种版本间的chain文件，下载方法很简单：

# hg19转hg38的chain文件 wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/liftOver/hg19ToHg38.over.chain.gz # hg38转hg19的chain文件（反向转换时使用） wget http://hgdownload.soe.ucsc.edu/goldenPath/hg38/liftOver/hg38ToHg19.over.chain.gz

这里有个实用小技巧：我习惯把常用的chain文件都下载到本地建立一个资源库，比如~/genome_liftover/目录下。这样下次需要时就不用重复下载了。记得用gunzip解压后再使用，可以提升转换速度。

2.2 安装UCSC工具集

liftover工具包含在UCSC的工具集中，可以通过以下命令安装：

# 下载UCSC工具集 wget http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/liftOver chmod +x liftOver

实测发现，不同操作系统的二进制文件可能不兼容。我在CentOS 7和Ubuntu 20.04上都测试过，建议使用Linux系统进行操作。如果遇到权限问题，记得用chmod +x给执行权限。

3. bed文件转换实战

3.1 基本转换命令

准备好chain文件和liftover工具后，bed文件转换就很简单了：

./liftOver input_hg19.bed hg19ToHg38.over.chain.gz output_hg38.bed unmapped.bed

这个命令会产生两个输出文件：

• output_hg38.bed：转换成功的坐标
• unmapped.bed：无法映射的坐标

我处理过的一个RNA-seq项目中，约12%的peak区域无法映射到新版本。这时候需要手动检查unmapped.bed文件，有些可能是因为基因组结构发生了重大改变。

3.2 处理转换失败的区域

对于无法映射的区域，我有几个实用建议：

1. 尝试放宽bed文件的范围，比如上下游各扩展50bp
2. 检查是否是新版本中填补的gap区域
3. 使用IGV基因组浏览器可视化查看具体变化

这里分享一个实际案例：我在处理ChIP-seq数据时，有个重要的转录因子结合位点无法映射。后来发现是因为hg38在这个区域填补了一个原先的gap，导致坐标系统完全改变。最后通过BLAST比对找到了对应的新位置。

4. vcf文件转换进阶操作

4.1 Picard工具安装

vcf文件转换推荐使用Picard工具，它是Broad Institute开发的Java工具包：

# 下载Picard wget https://github.com/broadinstitute/picard/releases/download/2.27.5/picard.jar

建议使用最新版本，我测试过2.27.5版本在hg19-hg38转换中表现稳定。记得提前安装Java运行环境（Java 8+）。

4.2 准备参考基因组

转换vcf时需要新版本的参考基因组fasta文件：

# 下载hg38参考基因组 wget http://hgdownload.soe.ucsc.edu/goldenPath/hg38/bigZips/hg38.fa.gz gunzip hg38.fa.gz # 建立索引 samtools faidx hg38.fa

这里有个坑我踩过：一定要确保fasta文件和chain文件的版本完全匹配。有一次我用UCSC的chain文件搭配GRCh38的fasta文件，结果转换完全失败。

4.3 执行vcf转换

完整的Picard转换命令如下：

java -jar picard.jar LiftoverVcf \ I=input_hg19.vcf.gz \ O=output_hg38.vcf.gz \ CHAIN=hg19ToHg38.over.chain.gz \ REJECT=unmapped_variants.vcf \ R=hg38.fa \ CREATE_INDEX=true

关键参数说明：

• CREATE_INDEX=true：自动创建索引文件
• REJECT：保存无法映射的变异
• R：指定新版本的参考基因组

5. 转换后的质量检查

完成转换后，强烈建议做以下检查：

1. 统计转换成功率：比对input和output的记录数
2. 随机抽查若干位点，用IGV确认坐标正确性
3. 检查unmapped文件，分析失败原因

我在处理千人基因组项目的vcf时，发现转换后约7%的SNP丢失。进一步分析发现，其中大部分是hg38新版本中修正的错误位点，这反而是好事——说明数据质量得到了提升。

6. 常见问题解决方案

问题1：转换后部分区域坐标偏移较大解决：这通常是由于基因组结构变化导致。建议使用BLAST比对确认新位置，或者考虑使用保守性更高的区域。

问题2：Picard报错"Invalid variant position"解决：先检查vcf文件是否有效（用bcftools view），有时是因为vcf头信息不完整。

问题3：链特异性数据转换解决：对于RNA-seq等链特异性数据，转换后要特别注意链的方向。可以使用bedtools的flank命令辅助处理。

记得转换完成后，一定要更新分析流程中的所有注释文件（如基因注释、功能预测等），确保使用同一版本的基因组数据。我见过不少项目因为混用不同版本的注释数据而导致分析结果不可靠。