以前做过一期TopHat用法的推文,但是在TopHat的网站上可知,如下图,TopHat的最后一个版本TopHat 2.1.1更新于2/23/2016,并且官方团队已经明确声明,后续将大力发展HISAT2。最新的HISAT2 2.1.0版本推出于6/8/2017,相比于TopHat,HISAT2有着相同的核心功能并且更加高效,所以我们有必要学习一下HISAT2的用法。
HISAT2的使用可分成三步:建立基因组index,reads比对到基因组和sort排序。
我们首先要从网上下载人类的注释文件human.gtf 和基因组信息hg38.fa
先从human.gtf中提取外显子和剪切位点信息。
extract_splice_sites.py human.gtf > human.ss
extract_exons.py human.gtf > human.exon
然后我们需要使用hisat2-build进行建库:
hisat2-build -p 5 --ss humangenome.ss --exon humangenome.exon hg38.fa hg38_tran
除了自己建库以外,HISAT2官方网站也给出了可以直接使用的Indexes下载地址:http://ccb.jhu.edu/software/hisat2/index.shtml。由于自己建库的时间比较长,这里笔者建议大家首选官方的Indexes。
假设我要用的双端测序样本为A_1_trimmed.fq 和A_2_trimmed.fq,官方下载的index在grc目录下:
hisat2 -p 5 --dta -x ./grc/genome -1 A_1_trimmed.fq -2 A_2_trimmed.fq -S A.sam > A.txt
我们使用了5个线程运行程序,这一步的运行速度是要远快于TopHat的,并且只生成一个sam文件,我们可以将运行日志保存在A.txt文件里。
TopHat的结果是会自动排序的,这样排序后的结果可以直接用于下一步的reads组装。但是HISAT2是需要手动对sam文件进行排序的:
samtools sort -@ 5 -o A.bam A.sam
我们使用5个线程对sam文件排序,并将结果保存为bam二进制格式。
就像之前我们使用的TopHat+cufflinks一样,我们同样可以使用HISAT2+cufflinks进行转录本分析,可以用HISAT2得到的A.bam文件作为cufflinks的输入文件进行基因表达量矩阵计算和差异表达分析。另外,我们也可以用最新的套装HISAT2+StringTie做分析。
TopHat和HISAT2出于同一个团队,HISAT2在2016年更新了4次,今年到现在为止更新了一次,可以看到该团队对HISAT2的支持力度是很大的,相信如果有什么比对算法上的改进也是放在HISAT2上面,所以这里笔者建议大家用HISAT2代替tophat作为新一代的比对软件。
这期就到这里了,下期再见 ~^o^~