使用liftover创建注释Chain文件（基因组坐标转换）

问题情形：

在自己组装出一个基因组之后，会想要看自己的序列的组装情况和染色体朝向，或者是想根据已有的基因组注释文件创建自己的基因组的注释文件。这个时候需要的是一个chain file来进行坐标的转换。
一般在ucsc的官网下会有一些不同版本的基因组的chain file，有时我们需要利用自己的序列产生自己的chainfile，这里对liftover的使用方法进行记录。

主要参考为：http://genomewiki.ucsc.edu/index.php/DoSameSpeciesLiftOver.pl

这里我们使用自己组装的序列asm.fa和hg19来创建，为了避免不必要的bug，我们按照文档说明构建文件夹和操作。

一.构建程序环境

####从官网下载需要的脚本和程序####
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集群环境：
参考：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

创建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})

运行liftover的时候，可能会出现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

最后编辑于：2019.03.28 16:12:43