轉(zhuǎn)載:
https://biozx.top/hisat2.html
HISAT2是TopHat2/Bowti2的繼任者,使用改進的BWT算法,實現(xiàn)了更快的速度和更少的資源占用。TopHat2已經(jīng)就不再維護,作者推薦TopHat2/Bowti2和HISAT的用戶轉(zhuǎn)換到HISAT2。HISAT利用大量FM索引,以覆蓋整個基因組。下面一段話是tohat2官網(wǎng)的說明:
HISAT2 is a successor to both HISAT and TopHat2\. We recommend that the HISAT and TopHat2 users switch to HISAT2.
Hisat2安裝
如果是下載source code的話還需要make編譯才能用,如下下載的是編譯好的,直接可以用。
wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/downloads/hisat2-2.1.0-Linux_x86_64.zip
unzip hisat2-2.1.0-Linux_x86_64.zip
echo 'PATH=$PATH:/your/hisat2/path' >> ~/.bashrc
source ~/.bashrc
Hisat2的使用
Hisat2 官網(wǎng)有使用教程。Hisat2壓縮包里面有個example目錄,下面的操作以example為例做下示范。
第一步:建立索引
HISAT2建立基因組索引時還支持將SNP信息以及轉(zhuǎn)錄組信息加入到索引中,這樣比對的時候就可以考慮SNP的情況和轉(zhuǎn)錄本情況。生成8個索引文件,以.ht2為后綴。Hisat2官網(wǎng)上有建好的索引可下載,并且壓縮包中有scripts/目錄記錄了是如何建立索引的。命令如下:
hisat2-build -p 4 genome.fa –snp genome.snp –ss genome.ss –exon genome.exon genome_snp_tran
- -p 4 : number of threads。
4個進程 - genome.fa: 基因組參考序列是必須的,多個參考序列文件請用逗號分隔。
- --snp: SNP 信息文件。
- --ss:
可變剪切位點信息文件。Splice site file name - --exon:
外顯子信息文件。Exon file name - genome_snp_tran:
索引文件名的前綴
其中,剪切位點文件和外顯子文件可從轉(zhuǎn)錄本gtf信息文件提取,Hisat2提供了提取腳本,提取如下:
extract_exons.py Homo_sapiens.GRCh38.83.chr.gtf > genome.exon
extract_splice_sites.py Homo_sapiens.GRCh38.83.chr.gtf > genome.ss
提取snp位點信息文件的腳本如下:
extract_snps.py snp142Common.txt > genome.snp
第二步、將reads比對到參考序列上
當建好索引以后,就可以將reads比對到參考序列上。這里又分為Paired-end和unpaired 兩種:
unpaired reads比對
[bio@ubuntumyindex]$ hisat2 -f -x 22_20-21M_snp -U ../reads/reads_1.fa -S eg1.sam
- -f:指query input files 文件類型是
fasta文件類型 - -x: 指
索引文件的前綴(看上一步建立索引的前綴) - -U:指input reads文件
- -S:指
輸出結(jié)果到哪里,輸出到SAM文件
Paired-end
[bio@ubuntumyindex]$ hisat2 -f -x 22_20-21M_snp -1 ../reads/reads_1.fa -2 ../reads/reads_2.fa -S eg2.sam
第三步、用samtools做下游分析
詳盡教程 http://samtools.sourceforge.net/mpileup.shtml
1.將SAM文件轉(zhuǎn)為BAM文件
samtools view -bS eg2.sam > eg2.bam
2.將BAM文件排序
samtools sort eg2.bam -o eg2.sorted.bam
3.生成VCF文件
samtools mpileup -uf ../reference/22_20-21M.fa eg2.sorted.bam|bcftools view - >eg2.raw.bcf