LiftOver chain 文件创建流程

当我们所研究的物种有多个基因组组装版本的时候,我们往往需要将两个版本的位置对应起来。比如某个版本注释了某些转录因子,而另外一个版本没有注释,我们需要知道这些转录因子在另一个版本的位置。比如早期的QTL研究是基于老版本的参考序列,我们需要知道这些QTL在新版本参考序列中的位置等等。

最常用的工具是UCSC的LiftOver,如果你所研究的物种刚好是LiftOver支持的物种,那很幸运,可以直接使用其在线工具,或者下载chain文件后,通过pyliftoverCrossMap等工具转换。如果没有做好的chain文件,那需要自己生成。UCSC对于自行生成chain文件有详细的pipeline,所以这篇博客不对每一步作详细说明,只提供主要步骤说明,并避开UCSC流程中的一些坑。

  • BLAT比对

既然要对应两个序列,那最直接的方式就是比对了,liftover原理主要基于此,并把位置关系存在chain文件中,用于后续查询。参考序列相近,比如同个物种的两个不同组装版本,可以用blat比对,如果是不同物种,则需要用blastz进行比对。

我对UCSC的pipeline进行了一些修改,命名为chain_step*.sh,中间需要修改或注意的地方加了注释。在运行之前,需要 1) 下载并编译好UCSC的kentUtils,liftover主要用到其中一系列工具。2) 将需要对应的两个参考序列转为2bit格式,并放在同一个目录下,可以用kentUtils的faToTwoBit进行转换。3)修改chain_step1.sh中的必要参数。主要需要修改的是kentUtils路径,工作路径,两个2bit格式参考序列前缀。需要指出的是,targetDb为已知位置的序列,queryDb为查询序列,即将targetDb的位置转为queryDb的位置。

运行该步后,会在run.blat目录下生成一个jobList,如果任务少,可以直接sh jobList运行,如果任务多,可以在集群上用SGE投递任务。运行需要blatJob.csh脚本,注意修改blatJob.csh中blat和kentUtils路径。

最大的坑出现了!!!

jobList里面的任务运行完之后,并不能直接运行流程接下来的部分,这是UCSC流程最大的坑。在run.blat/psl/*/目录下,生成了一系列的psl格式的比对结果,需要对该结果进行修改,主要是因为对targetDb序列进行了切割,psl文件里面记录的染色体名和长度以及比对的位置都需要进行调整。该调整可以通过脚本change_psl.py实现,如需要对所有psl文件进行调整,则在psl目录下运行ls */*.psl|while read a;do python change_psl.py $a target.fa.fai ;done,其中target.fa.fai是待转换参考序列fa格式的fai文件,通过samtools faidx target.fa生成。

  • 生成chain文件

第一步完成后,接下来就比较简单了。在工作目录下运行chain_step2.sh后,会在run.chain目录下生成jobList,直接运行或拆分后用SGE投递。
jobList任务运行结束后,在工作目录运行chain_step3.sh,会在run.chain目录下生成doNet.csh,直接运行该脚本。注意修改step2和step3脚本中的路径。

运行结束后,会在工作目录生成后缀为.over.chain.gz的压缩文件,该文件就是我们需要的chain文件。

至此,大功告成!

利用该文件,就可以使用CrossMap对bed/gtf/gff/bam/vcf等文件中的参考序列位置进行转换。

如果我的博客能帮你节约时间,我会非常高兴你请我喝杯咖啡。