微生信分享 http://blog.sciencenet.cn/u/chinapubmed 0代码在线绘制160+科研图

博文

做研究,该用hg19还是hg38基因组?一行代码将hg19的bw文件转成hg38

已有 3615 次阅读 2022-8-28 09:49 |个人分类:生物信息|系统分类:科研笔记

有没有想过这个问题:做研究的时候该用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号染色体100200区域中的信号是5.2

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

fig1.png

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

 

其中bigWigToBedGraphfetchChromSizesbedGraphToBigWig都可以在UCSC下载

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

 

最后,让我们来检验下hg19hg38的转换效果吧:

fig2.jpg

2. Hg19可视化

fig3.jpg

3. Hg38可视化

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

注意:有些位点会转化失败,因此最佳解决方案还是重新使用hg38基因组从原始数据进行处理。

 

回答开头的问题:

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

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

 

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



https://blog.sciencenet.cn/blog-707141-1352956.html

上一篇:在线绘制圆形miRNA-靶基因,转录因子-gene等网络图
下一篇:分久必合 -- 合并相同基因/探针的表达量
收藏 IP: 117.143.114.*| 热度|

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-4-28 01:21

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部