||
有没有想过这个问题:做研究的时候该用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+
Archiver|手机版|科学网 ( 京ICP备07017567号-12 )
GMT+8, 2024-9-2 03:27
Powered by ScienceNet.cn
Copyright © 2007- 中国科学报社