問題情形:
在自己組裝出一個基因組之后们颜,會想要看自己的序列的組裝情況和染色體朝向脑慧,或者是想根據(jù)已有的基因組注釋文件創(chuàng)建自己的基因組的注釋文件幌衣。這個時候需要的是一個chain file來進行坐標的轉(zhuǎn)換了嚎。
一般在ucsc的官網(wǎng)下會有一些不同版本的基因組的chain file葱色,有時我們需要利用自己的序列產(chǎn)生自己的chainfile,這里對liftover的使用方法進行記錄抹镊。
主要參考為:http://genomewiki.ucsc.edu/index.php/DoSameSpeciesLiftOver.pl
這里我們使用自己組裝的序列asm.fa和hg19來創(chuàng)建锉屈,為了避免不必要的bug,我們按照文檔說明構(gòu)建文件夾和操作垮耳。
一.構(gòu)建程序環(huán)境
####從官網(wǎng)下載需要的腳本和程序####
mkdir /home/jfh/projects/hic_analysis/data/liftover/data/bin
mkdir /home/jfh/projects/hic_analysis/data/liftover/data/scripts
chmod 755 /home/jfh/projects/hic_analysis/data/liftover/data/bin
chmod 755 /home/jfh/projects/hic_analysis/data/liftover/data/scripts
cd /home/jfh/projects/hic_analysis/data/liftover/data
rsync -a rsync://hgdownload.soe.ucsc.edu/genome/admin/exe/linux.x86_64/ ./bin
git archive --remote=git://genome-source.soe.ucsc.edu/kent.git \
--prefix=kent/ HEAD src/hg/utils/automation \
| tar vxf - -C ./scripts --strip-components=5 \
--exclude='kent/src/hg/utils/automation/incidentDb' \
--exclude='kent/src/hg/utils/automation/configFiles' \
--exclude='kent/src/hg/utils/automation/ensGene' \
--exclude='kent/src/hg/utils/automation/genbank' \
--exclude='kent/src/hg/utils/automation/lastz_D' \
--exclude='kent/src/hg/utils/automation/openStack'
wget -O ./bin/bedSingleCover.pl \
'http://genome-source.soe.ucsc.edu/gitlist/kent.git/raw/master/src/utils/bedSingleCover.pl'
將程序路徑加入到~/.bashrc后面
export PATH=/home/jfh/projects/hic_analysis/data/liftover/data/bin:/data/scripts:$PATH
export PATH=/home/jfh/projects/hic_analysis/data/liftover/data/bin:/home/jfh/projects/hic_analysis/data/liftover/data/bin/blat:/home/jfh/projects/hic_analysis/data/liftover/data/scripts:$PATH
另外颈渊,需要搭建parasol集群環(huán)境:
參考:http://genomewiki.ucsc.edu/index.php/Parasol_job_control_system
需要注意的是:文檔中并未提及需要能免密ssh localhost
,需要自己手動配置免密ssh本地终佛,設置好后其余的按照文檔要求做即可俊嗽。
二.準備文件
將自己的序列分別壓縮為asm.fa.gz
和hg19.fa.gz
形成
/home/jfh/projects/hic_analysis/data/liftover/data/genomes/hg19/hg19.fa.gz
/home/jfh/projects/hic_analysis/data/liftover/data/genomes/asm/asm.fa.gz
準備工作
genome=/home/hugo/storage/jifh/projects/liftover/data/genomes
query=hg19
target=asm
cd $genome/asm
faToTwoBit $genome/genbank/asm.fa.gz asm.2bit
twoBitInfo asm.2bit stdout | sort -k2,2nr > asm.chrom.sizes
cd $genome/hg19
faToTwoBit $genome/refseq/hg19.fa.gz hg19.2bit
twoBitInfo hg19.2bit stdout | sort -k2,2nr > hg19.chrom.sizes
創(chuàng)建ooc file
blat asm.2bit /dev/null /dev/null -tileSize=11 -makeOoc=asm.ooc -repMatch=1024
三、運行Liftover程序
data=/home/hugo/storage/jifh/projects/liftover/data
export target="asm"
export query="hg19"
cd $data/genomes/${target}
time (doSameSpeciesLiftOver.pl -verbose=2 -buildDir=`pwd` \
-ooc=`pwd`/${target}.ooc -fileServer=localhost -localTmp="/dev/shm" \
-bigClusterHub=localhost -dbHost=localhost -workhorse=localhost \
-target2Bit=`pwd`/${target}.2bit -targetSizes=`pwd`/${target}.chrom.sizes \
-query2Bit=$data/genomes/${query}/${query}.2bit \
-querySizes=$data/genomes/${query}/${query}.chrom.sizes ${target} ${query})
運行l(wèi)iftover的時候铃彰,可能會出現(xiàn)parasol運行停止的情況绍豁,機器提示sick machine,這時只要進入到batch所在的目錄下牙捉,運行以下命令即可:
para clearSickNodes
para try
para push
最后在自己的target文件夾下面就可以看到asmToHg19.over.chain.gz竹揍,就可以拿去用了。
Reference Link:
http://genomewiki.ucsc.edu/index.php/DoSameSpeciesLiftOver.pl
http://genomewiki.ucsc.edu/index.php/Parasol_job_control_system
https://genecats.gi.ucsc.edu/eng/parasol.html