(14)不同基因坐标转换-生信菜鸟团博客2周年精选文章集

主流有3个,我只介绍了两个:

用crossmap代替liftover做基因组坐标转换 liftover基因组版本直接的coordinate转换

其实国际三大主流生物信息学数据库运营单位都出了自己的基因组坐标转换,它们分别是 (UCSC liftOver, NCBI Remap, Ensembl API)

Ensembl’s Assembly Converter.是基于crossmap的,我觉得挺好用的,就介绍给大家!!!

This online tool currently uses CrossMap, which supports a limited number of formats (see our online documentation for details of the individual data formats listed below). CrossMap also discards metadata in files, so track definitions, etc, will be lost on conversion.

Important note: CrossMap converts WIG files to BedGraph internally for efficiency, and also outputs them in BedGraph format.

但是不知道为什么UCSC的liftover最出名,我也写过它的教程,(http://www.bio-info-trainee.com/?p=990

还算好用,比较真正基因组坐标转换的需求很少,但是略有一些限制,所以我用起来了Ensembl的crossmap软件。

该软件说明书就在主页,很简单的:http://crossmap.sourceforge.net/

一、程序安装

http://iweb.dl.sourceforge.net/project/crossmap/test.hg19.zip

https://sourceforge.net/projects/crossmap/files/CrossMap-0.2.2.tar.gz/

直接wget到linux系统即可,解压之后就能看到该程序,是基于Python的,所以安装方式跟Python模块一样!

http://www.bio-info-trainee.com/?p=581

软件官网也有安装说明书

我的安装命令是:

python setup.py install –user

安装之后添加到环境变量

# setup PYTHONPATH. Skip this step if CrossMap was installed to default location.

$ export PYTHONPATH=/home/jmzeng/.local/lib/python2.7/site-packages:$PYTHONPATH.

二、输入数据准备

该软件提供了测试数据,也提供了基因组转换的模板chain格式文件:

它要求的待转换文件很广泛:

  1. BAM or SAM format.
  2. BED or BED-like format. BED file must has at least 3 columns (‘chrom’, ‘start’, ‘end’).
  3. Wiggle format. “variableStep”, “fixedStep” and “bedGraph” wiggle line are supported.
  4. BigWig format.
  5. GFF or GTF format.
  6. VCF format.

三、运行程序

找到你的crossmap安装目录,根据你的安装命令来的!

我的程序在/home/jmzeng/.local/bin 里面,直接可以使用

根据下面的例子,自己修改就可以了,程序用全路径调用,然后写明是什么文件格式,然后列出chain文件的地址,然后input和output即可

Example (run CrossMap with no output_file specified):

$ python CrossMap.py bed hg18ToHg19.over.chain.gz test.hg18.bed3

Conversion results were printed to screen directly (column1-3 are hg18 based, column5-7 are hg19 based):

chr1   142614848       142617697       ->      chr1    143903503       143906352
chr1   142617697       142623312       ->      chr1    143906355       143911970
chr1   142623313       142623350       ->      chr1    143911971       143912008
chr1   142623351       142626523       ->      chr1    143912009       143915181
chr1   142633862       142633883       ->      chr1    143922520       143922541
chr1   142633884       142636152       ->      chr1    143922542       143924810
chr1   142636152       142636326       ->      chr1    143924813       143924987
chr1   142636339       142636391       ->      chr1    143925000       143925052
chr1   142636392       142637362       ->      chr1    143925052       143926022
chr1   142637373       142639738       ->      chr1    143926033       143928398
chr1   142639739       142639760       ->      chr1    143928399       143928420
chr1   142639761       142640145       ->      chr1    143928421       143928805
chr1   142640153       142641149       ->      chr1    143928813       143929809

四、输出数据解读

输出数据没什么好解读的了,进去的是什么数据,出来的就是什么数据,只是把你的坐标进行了转换

liftover基因组版本直接的coordinate转换

Posted on 2015年9月7日

下载地址:http://hgdownload.cse.ucsc.edu/admin/exe/

我一般是使用linux版本的:wget http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/liftOver

使用方法:【从hg38转到hg19】

因为主流的基因组版本还是hg19,但是时代在进步,已经有很多信息都是以hg38的形式公布出来的了。

比如,我下载了pfam.df这个protein domain注释文件,对人的hg38基因组每个坐标都做了domain注释,数据形式如下:

查看文件内容head pfam.hg38.df ,如下:

PFAMID chr start end strand

Helicase_C_2 chr1 12190 12689 +

7tm_4 chr1 69157 69220 +

7TM_GPCR_Srsx chr1 69184 69817 +

7tm_1 chr1 69190 69931 +

7tm_4 chr1 69490 69910 +

7tm_1 chr1 450816 451557 -

7tm_4 chr1 450837 451263 -

EPV_E5 chr1 450924 450936 -

7TM_GPCR_Srsx chr1 450927 451572 -

我想把domain的起始终止坐标转换成hg19的,就必须要借助UCSC的liftover这个工具啦

这个工具需要一个坐标注释文件 http://hgdownload-test.cse.ucsc.edu/goldenPath/hg38/liftOver/

我这里需要下载的是http://hgdownload-test.cse.ucsc.edu/goldenPath/hg38/liftOver/hg38ToHg19.over.chain.gz

而且它只能对bed等符合要求的格式进行转换

http://www.ensembl.org/info/website/upload/bed.html

示例如下:

chr7  127471196  127472363  Pos1  0  +  127471196  127472363  255,0,0
chr7  127472363  127473530  Pos2  0  +  127472363  127473530  255,0,0

很简单的,把自己的文件随便凑几列信息,做成这个9列的格式即可

cat pfam.hg38.df |sed ‘s/\r//g’ |awk ‘{print $2,$3,$4,$1,0,$5,$3,$4,”255,0,0″}’ >pfam.hg38.bed

这样就有了足够的文件可以进行坐标转换啦,转换的命令非常简单!

chmod 777 liftOver

./liftOver pfam.hg38.bed hg38ToHg19.over.chain pfam.hg19.bed unmap

然后运行成功了会有 提示,报错一般是你的格式不符合标准bed格式,自己删掉注释行等等不符合的信息即可

Reading liftover chains

Mapping coordinates

转换后,稍微检查一下就可以看到坐标的确发生了变化,当然,我们只需要看前面几列信息即可

grep -w p53 *bed

pfam.hg19.bed:chr11 44956439 44959858 p53-inducible11 0 – 44956439 44959858 255,0,0

pfam.hg19.bed:chr11 44956439 44959767 p53-inducible11 0 – 44956439 44959767 255,0,0

pfam.hg19.bed:chr2 669635 675557 p53-inducible11 0 – 669635 675557 255,0,0

pfam.hg19.bed:chr22 35660826 35660982 p53-inducible11 0 + 35660826 35660982 255,0,0

仔细看看坐标是不是变化啦!

pfam.hg38.bed:chr11 44934888 44938307 p53-inducible11 0 – 44934888 44938307 255,0,0

pfam.hg38.bed:chr11 44934888 44938216 p53-inducible11 0 – 44934888 44938216 255,0,0

pfam.hg38.bed:chr2 669635 675557 p53-inducible11 0 – 669635 675557 255,0,0

pfam.hg38.bed:chr22 35264833 35264989 p53-inducible11 0 + 35264833 35264989 255,0,0

其实R里面的bioconductor系列包也可以进行坐标转换 http://www.bioconductor.org/help/workflows/liftOver/

这个可以直接接着下载pfam.df数据库来做下去。更方便一点。

我的数据如下,需要自己创建成一个GRanges对象
library(GenomicRanges)
pfam.hg38 <- GRanges(seqnames=Rle(a[,2]),
               ranges=IRanges(a[,3], a[,4]),
               strand=a[,5])

这样就OK拉,虽然这只是一个很简陋的GRanges对象,但是这个GRanges对象可以通过R的liftover方法来转换坐标啦。

library(rtracklayer)
ch = import.chain("hg38ToHg19.over.chain")

pfam.hg19 = liftOver(pfam.hg38, ch)

pfam.hg19 =unlist(pfam.hg19)

再把这个转换好的pfam.hg19 写出即可

参考:http://www.zilhua.com/906.html

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2017-01-06

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏机器之心

教程 | 只需15分钟,使用谷歌云平台运行Jupyter Notebook

选自Medium 机器之心编译 参与:路雪 近日,Amulya Aankul 在 Medium 上发表文章,描述他在谷歌云平台上运行 Jupyter Noteb...

4428
来自专栏游戏开发那些事

【Unity游戏开发】AssetBundle杂记--AssetBundle的二三事

  马三在公司大部分时间做的都是游戏业务逻辑和编辑器工具等相关工作,因此对Unity AssetBundle这块的知识点并不是很熟悉,自己也是有打算想了解并熟悉...

1012
来自专栏MasiMaro 的技术博文

WinSock 完成端口模型

之前写了关于Winsock的重叠IO模型,按理来说重叠IO模型与之前的模型相比,它的socket即是非阻塞的,也是异步的,它基本上性能非常高,但是它主要的缺点在...

771
来自专栏FreeBuf

WireShark+Winhex:流量分析的好搭档

这篇文章你将学会的知识点有 1、进阶的wireshark的流量分析、解码、追踪流、导出文件 2、利用hackbar进行base64、URL编码转换 3、利用wi...

5916
来自专栏祝威廉

如何使用MLSQL中的帮助指令学习模块的使用

MLSQL 已经实现了文章中描述的功能 如何实现语法的自解释(MLSQL易用性设计有感) 。

814
来自专栏猿人谷

对缓存的思考——提高命中率

开篇 编写高效的程序并不只在于算法的精巧,还应该考虑到计算机内部的组织结构,cpu微指令的执行,缓存的组织和工作原理等。 好的算法在实际中不见得有高效率,如果完...

2139
来自专栏北京马哥教育

爬虫大神,又出新招

1565
来自专栏QQ会员技术团队的专栏

Android 动态库压缩壳的实现

计算机软件领域所说的壳实际上是一种软件加密技术。壳主要分为两大类:加密壳和压缩壳,加密壳侧重于防止软件被篡改,而压缩壳则侧重于减小软件体积。其实,在Window...

1.5K1
来自专栏互联网杂技

内存卡存储原理,你知道吗?

1、 简介: SD卡(Secure Digital Memory Card)是一种为满足安全性、容量、性能和使用环境等各方面的需求而设计的一种新型存储器...

3706
来自专栏生信宝典

别人的电子书,你的电子书,都在bookdown

bookdown是著名R包作者谢益辉开发的,支持采用Rmarkdown (R代码可以运行)或普通markdown编写文档,然后编译成HTML, WORD, PD...

40511

扫码关注云+社区