(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格式文件:
它要求的待转换文件很广泛:
- BAM or SAM format.
- BED or BED-like format. BED file must has at least 3 columns (‘chrom’, ‘start’, ‘end’).
- Wiggle format. “variableStep”, “fixedStep” and “bedGraph” wiggle line are supported.
- BigWig format.
- GFF or GTF format.
- 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
- JavaScript 教程
- JavaScript 编辑工具
- JavaScript 与HTML
- JavaScript 与Java
- JavaScript 数据结构
- JavaScript 基本数据类型
- JavaScript 特殊数据类型
- JavaScript 运算符
- JavaScript typeof 运算符
- JavaScript 表达式
- JavaScript 类型转换
- JavaScript 基本语法
- JavaScript 注释
- Javascript 基本处理流程
- Javascript 选择结构
- Javascript if 语句
- Javascript if 语句的嵌套
- Javascript switch 语句
- Javascript 循环结构
- Javascript 循环结构实例
- Javascript 跳转语句
- Javascript 控制语句总结
- Javascript 函数介绍
- Javascript 函数的定义
- Javascript 函数调用
- Javascript 几种特殊的函数
- JavaScript 内置函数简介
- Javascript eval() 函数
- Javascript isFinite() 函数
- Javascript isNaN() 函数
- parseInt() 与 parseFloat()
- escape() 与 unescape()
- Javascript 字符串介绍
- Javascript length属性
- javascript 字符串函数
- Javascript 日期对象简介
- Javascript 日期对象用途
- Date 对象属性和方法
- Javascript 数组是什么
- Javascript 创建数组
- Javascript 数组赋值与取值
- Javascript 数组属性和方法
- Spring Boot Actuators
- 用图机器学习探索 A 股个股相关性变化
- Python 技术篇-连接oracle数据库并执行sql语句实例演示,python连接oracle数据库oci详细配置方法
- MySQL 技术篇-mysql数据库的安装、配置与使用实例演示
- JavaScript 技术篇-js代码获取当前操作系统信息、浏览器版本信息实例演示,windows NT版本对照表
- Oracle 数据库impdp导入数据库版本和dmp数据库文件版本不匹配问题解决方法,ORA-39142版本号不兼容、ORA-39000转储文件说明错误解决方法
- 实践总结:基于Kbone使用React同构开发小程序
- BAT 批处理命令 - 实现输出当前文件夹下的所有文件夹名的功能实例演示
- Python+Selenium 自动化-指定chrome驱动运行selenium实例演示,运行指定位置下的浏览器驱动
- Linux 命令查找指定文件夹下符合查询条件的文件和文件夹实例演示
- 用Python实现一个最新QQ办公版(TIM)的登录界面
- Oracle 数据库直接执行本地sql文件、sql脚本实例演示
- Oracle 数据库利用回收站恢复删除的表实例演示
- Linux 命令利用scp实现从服务器共享地址上传下载文件、文件夹实例演示,scp命令的参数详解
- Oracle 数据库利用sql语句判断某个表是否是临时表实例演示,达梦数据库查询出所有临时表