取代bowtie+tophat进行RNA-seq比对
HISAT全称为Hierarchical Indexing for Spliced Alignment of Transcripts,由约翰霍普金斯大学开发。它取代Bowtie/TopHat程序,能够将RNA-Seq的读取与基因组进行快速比对。这项成果发表在3月9日的《Nature Methods》上。
HISAT利用大量FM索引,以覆盖整个基因组。以人类基因组为例,它需要48,000个索引,每个索引代表~64,000 bp的基因组区域。这些小的索引结合几种比对策略,实现了RNA-Seq读取的高效比对,特别是那些跨越多个外显子的读取。尽管它利用大量索引,但HISAT只需要4.3 GB的内存。这种应用程序支持任何规模的基因组,包括那些超过40亿个碱基的。
HISAT软件可从以下地址获取:http://ccb.jhu.edu/software/hisat/index.shtml。
首先,我们安装这个软件!
Wget http://ccb.jhu.edu/software/hisat/downloads/hisat-0.1.5-beta-source.zip
官网下载的是源码包,需要make一下,make之后目录下面就多了很多程序,绿色的那些都是,看起来是不是很眼熟呀!!!
哈哈,这完全就是bowtie的模拟版本!!!
也可以从github里面下载,wget https://codeload.github.com/infphilo/hisat/zip/master
下载后直接解压即可使用啦。当然这个软件本身也有着详尽的说明书
http://ccb.jhu.edu/software/hisat/manual.shtml
然后就是准备数据,它跟tophat一样的功能。就是把用RNA-seq方法测序得到的fastq文件比对到参考基因组上面,所以就准这两个文件了哦
接下来是运行程序!
说明书上面写着分成两个步骤,构建索引和比对。
这个软件包模仿bowtie自带了一个example数据,而且它的说明书也是针对于那个example来的,我也简单运行一下。
$HISAT_HOME/hisat-build $HISAT_HOME/example/reference/22_20-21M.fa 22_20-21M_hisat
构建索引的命令如上,跟bowtie一样我修改了一下
/home/jmzeng/hoston/RNA-soft/hisat-0.1.5-beta/hisat-build 22_20-21M.fa my_hisat_index
连日志都跟bowtie一模一样,哈哈,可以看到我们的这个参考fasta文件 22_20-21M.fa 就变成索引文件啦,索引还是很多的!
然后就是比对咯,还是跟bowtie一样
$HISAT_HOME/hisat -x 22_20-21M_hisat -U $HISAT_HOME/example/reads/reads_1.fq -S eg1.sam
我的命令是
/home/jmzeng/hoston/RNA-soft/hisat-0.1.5-beta/hisat -x my_hisat_index -U ../reads/reads_1.fq -S reads1.sam
1000 reads; of these:
1000 (100.00%) were unpaired; of these:
0 (0.00%) aligned 0 times
1000 (100.00%) aligned exactly 1 time
0 (0.00%) aligned >1 times
100.00% overall alignment rate
哈哈,到这里。这个软件就运行完毕啦!!!是不是非常简单,只有你会用bowtie,这个就没有问题。当然啦,软件还是有很多细节是需要调整的。我下面就简单讲一个实际的例子哈!
首先,我用了1.5小时把4.6G的小鼠基因组构建了索引
/home/jmzeng/hoston/RNA-soft/hisat-0.1.5-beta/hisat-build Mus_musculus.GRCm38.fa.fa mouse_hisat_index
然后对我的四个测序文件进行比对。
for i in *fq
do
/home/jmzeng/hoston/RNA-soft/hisat-0.1.5-beta/hisat -x /home/jmzeng/hoston/mouse/mouse_hisat_index \
-p 30 -U $i.trimmed.single -S ./hisat_out/${i%.*}.sam
done
它运行的速度的确要比tophat快好多,太可怕的速度!!!!至于是否多消耗了内存我就没有看了
4.6G的小鼠,5G的测序数据,我只用了五个核,居然十分钟就跑完了!
然后听群友说是因为没有加 --known-splicesite-infile <path>这个参数的原因,没有用gtf文件来指导我们的RNA数据的比对,这样是不对的!
需要用下面这个脚本把gtf文件处理一下,然后导入什么那个参数来指导RNA比对。
extract_splice_sites.py genes.gtf > splicesites.txt
但是我报错了,错误很奇怪,没解决,但是我换了个 extract_splice_sites.py 程序,就可以运行啦!之前是HISAT 0.1.5-beta release 2/25/2015里面的python程序,后来我换做了github里面的就可以啦!
/home/jmzeng/hoston/RNA-soft/hisat-master/extract_splice_sites.py Mus_musculus.GRCm38.79.gtf >mouse_splicesites.txt
21192819 reads; of these:
21192819 (100.00%) were unpaired; of these:
14236834 (67.18%) aligned 0 times
5437800 (25.66%) aligned exactly 1 time
1518185 (7.16%) aligned >1 times
感觉没有变化,不知道为什么?
21192819 reads; of these:
21192819 (100.00%) were unpaired; of these:
14236838 (67.18%) aligned 0 times
5437793 (25.66%) aligned exactly 1 time
1518188 (7.16%) aligned >1 times
32.82% overall alignment rate
发表这个软件的文献本身也把这个软件跟其它软件做了详尽的对比
http://www.nature.com/nmeth/journal/v12/n4/full/nmeth.3317.html
Program | Run time (min) | Memory usage (GB) |
Run times and memory usage for HISAT and other spliced aligners to align 109 million 101-bp RNA-seq reads from a lung fibroblast data set. We used three CPU cores to run the programs on a Mac Pro with a 3.7 GHz Quad-Core Intel Xeon E5 processor and 64 GB of RAM. | ||
HISATx1 | 22.7 | 4.3 |
HISATx2 | 47.7 | 4.3 |
HISAT | 26.7 | 4.3 |
STAR | 25 | 28 |
STARx2 | 50.5 | 28 |
GSNAP | 291.9 | 20.2 |
OLego | 989.5 | 3.7 |
TopHat2 | 1,170 | 4.3 |
参考:http://www.plob.org/2015/03/20/8980.html
http://nextgenseek.com/2015/03/hisat-a-fast-and-memory-lean-rna-seq-aligner/