hisat2比对软件将reads比对到参考基因组

转载：
https://biozx.top/hisat2.html

HISAT2是TopHat2/Bowti2的继任者，使用改进的BWT算法，实现了更快的速度和更少的资源占用。TopHat2已经就不再维护，作者推荐TopHat2/Bowti2和HISAT的用户转换到HISAT2。HISAT利用大量FM索引，以覆盖整个基因组。下面一段话是tohat2官网的说明：

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 官网有使用教程。Hisat2压缩包里面有个example目录，下面的操作以example为例做下示范。

第一步：建立索引

HISAT2建立基因组索引时还支持将SNP信息以及转录组信息加入到索引中，这样比对的时候就可以考虑SNP的情况和转录本情况。生成8个索引文件，以.ht2为后缀。Hisat2官网上有建好的索引可下载，并且压缩包中有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：索引文件名的前缀

image

其中，剪切位点文件和外显子文件可从转录本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：指输出结果到哪里，输出到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文件转为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