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

使用方法:【從hg38轉(zhuǎn)到hg19】

因?yàn)橹髁鞯幕蚪M版本還是hg19，但是時(shí)代在進(jìn)步，已經(jīng)有很多信息都是以hg38的形式公布出來的了。

比如，我下載了pfam.df這個(gè)protein domain注釋文件，對人的hg38基因組每個(gè)坐標(biāo)都做了domain注釋，數(shù)據(jù)形式如下：

查看文件內(nèi)容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的起始終止坐標(biāo)轉(zhuǎn)換成hg19的，就必須要借助UCSC的liftover這個(gè)工具啦

而且它只能對bed等符合要求的格式進(jìn)行轉(zhuǎn)換

示例如下：

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

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

很簡單的，把自己的文件隨便湊幾列信息，做成這個(gè)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

這樣就有了足夠的文件可以進(jìn)行坐標(biāo)轉(zhuǎn)換啦，轉(zhuǎn)換的命令非常簡單！

chmod 777 liftOver

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

然后運(yùn)行成功了會(huì)有 提示，報(bào)錯(cuò)一般是你的格式不符合標(biāo)準(zhǔn)bed格式，自己刪掉注釋行等等不符合的信息即可

Reading liftover chains

Mapping coordinates

轉(zhuǎn)換后，稍微檢查一下就可以看到坐標(biāo)的確發(fā)生了變化，當(dāng)然，我們只需要看前面幾列信息即可

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

仔細(xì)看看坐標(biāo)是不是變化啦！

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

這個(gè)可以直接接著下載pfam.df數(shù)據(jù)庫來做下去。更方便一點(diǎn)。

我的數(shù)據(jù)如下，需要自己創(chuàng)建成一個(gè)GRanges對象

library(GenomicRanges)

pfam.hg38

ranges=IRanges(a[,3], a[,4]),

strand=a[,5])

這樣就OK拉，雖然這只是一個(gè)很簡陋的GRanges對象，但是這個(gè)GRanges對象可以通過R的liftover方法來轉(zhuǎn)換坐標(biāo)啦。

library(rtracklayer)

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

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

pfam.hg19 =unlist(pfam.hg19)

再把這個(gè)轉(zhuǎn)換好的pfam.hg19 寫出即可

總結(jié)

以上是生活随笔為你收集整理的linux+基因组字符替换,liftover基因组版本直接的coordinate转换的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。