下載地址: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 寫出即可

總結

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

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。