做研究用hg19orhg38基因组？一行代码将bw转hg38

有没有想过这个问题：做研究的时候该用hg19还是hg38基因组？

人类基因组版本现状

对于同一个物种，数据库中存在不同的基因组版本，以人类（Homo Sapiens）为例，UCSC基因组浏览器中有多个版本：Dec. 2013 (GRCh38/hg38)、Feb. 2009 (GRCh37/hg19)、Mar. 2006 (NCBI36/hg38)等，括号前面的是组装（Assembly）日期。是不是有点惊讶！都2022年了，默认的基因组还是9年前的，更令人惊讶的是，好多网站现在默认使用的还是hg19，也就是13年前的基因组版本，更离谱的是，GEO数据库中存在大量hg19，甚至hg18的数据集。那么在GEO数据挖掘的在日常工作中，经常会遇见以下两种场景：

1）hg19 -> hg38，例如文献中使用的是hg19，你自己的测序数据是hg38

2）hg38 -> hg19，例如你师兄的师兄留给你的数据是hg19，而文献是hg38

解决方案：UCSC提供的一个工具liftover

Liftover简介

http://genome.ucsc.edu/cgi-bin/hgLiftOver

将不同版本的染色体上的位置一一对应，官方定义为：This tool converts genome coordinates and genome annotation files between assemblies.

不难想象，该工具至少需要3个参数：hg19的坐标文件，hg38的坐标文件，以及两者之间关系的数据库文件（chain文件）

bw文件简介

bw文件是bigwig的缩写，存储了坐标及对应信号信息。然而，bw是一种二进制文件，不能用liftover直接处理。因此，要将hg19基因组的bw文件转成hg38，需要找个代理，这个代理就是bedgraph文件，包含4列，例如

chr1 100 200 5.2

表示：1号染色体100到200区域中的信号是5.2

bedgraph格式可以被liftover直接处理。

图1. 转换示意图

转换代码

前人栽树后人乘凉，python 包CrossMap可以直接处理bw文件的转化问题。

因此一行代码的转化过程如下：

1，安装CrossMap

pip install CrossMap

2，下载hg19-hg39转化对应的数据库文件

http://hgdownload.cse.ucsc.edu/goldenpath/hg19/liftOver/hg19ToHg38.over.chain.gz

3，一行代码转化

CrossMap.py bigwig hg19ToHg38.over.chain input.bw output.bw

然后就可以导入到IGV进行查看和比较了。

当然，也可以逐步进行

bigWigToBedGraph input.bw input.bedGraph

liftOver input.bedGraph hg19ToHg38.over.chain input_hg38.bedgraph

fetchChromSizes hg38 > hg38.chrom.sizes

sort -k1,1 -k2,2n input_hg38.bedgraph > input_hg38.sorted.bedgraph

bedGraphToBigWig input_hg38.sorted.bedgraph hg38.chrom.sizes output.bw

其中bigWigToBedGraph、fetchChromSizes、bedGraphToBigWig都可以在UCSC下载

http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/

最后，让我们来检验下hg19和hg38的转换效果吧：

图2. Hg19可视化

图3. Hg38可视化

一模一样，肉眼看不出差别，说明结果正确。

注意：有些位点没有对应关系的话，会转化失败，因此最佳解决方案还是使用hg38基因组从原始数据重新处理。

回答开头的问题：

现在包括UCSC、TCGA、Ensembl等大型数据库均以hg38作为默认基因组，因此，用hg38就对了，还在用hg19的研究者，请赶紧升级！

数据分析的时候，一定要看清楚，网上的数据到底是hg38还是hg19！因为除了大型数据库外，其他的小型数据库、网站经常是发了文章就不再更新，甚至是发了文章6个月以后就找不到那种，一定要“三思而后行”，避免浪费时间！

微生信助力发文章，谷歌引用590+，知网引用450+