linux+基因组字符替换,liftover基因组版本直接的coordinate转换

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

使用方法:【从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这个工具啦

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

示例如下：

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

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

我的数据如下，需要自己创建成一个GRanges对象

library(GenomicRanges)

pfam.hg38

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 写出即可