一、Hisat2介绍
Hisat是一种高效的RNA-seq实验比对工具。它使用了基于BWT和Ferragina-manzini (Fm) index 两种算法的索引框架。使用了两类索引去比对，一类是全基因组范围的FM索引来锚定每一个比对，另一类是大量的局部索引对这些比对做快速的扩展。比对原理可阅读文献原文：HISAT: a fast spliced aligner with low memory requirements.
二、Hisat2设计原则
1.优化了索引建立策略
hisat2应用了基于bowtie2的方法去处理很多低水平的用于构建和查询FM索引的操作。但是与其它比对器不同的是，该软件应用了两类不同的索引类型：代表全基因组的全局FM索引和大量的局部小索引，每个索引代表64000bp的区域。以人类基因组为例，创建了约48000个局部索引，每一个索引与其相邻索引重叠1024bp，最终可以覆盖这个30亿个碱基的基因组。这种存在交叉（overlap）的边界可以轻松的比对那些跨区域的read（可变剪切体）。尽管有很多索引，但是hisat会把他们使用合适方法压缩，最终只占4GB左右的内存。
2.采用了新的比对策略
RNA-seq产生的reads可能跨越长度比较大的内含子，在哺乳动物中甚至最长能达到1MB，同时外显子比较短，read也比较短。当使用100-bp reads时，会有很多read（模拟数据中大概34%）跨两个外显子的情况。为了更好的比对，将跨外显子的reads分成了三类：1）长锚定read，两个外显子中的每个都至少有16 bp；2）中间锚定read，一个外显子中具有8-15bp；3）短锚定read，仅与其中一个外显子比对上1-7bp。所以总的reads可以被划分为五类：1）不跨外显子的read 2）长锚定read 3）中间锚定read 4）短锚定read 5）跨两个外显子以上的read。在模拟的数据中，有25%左右的read是长锚定read，这种read在大多数情况下可以被唯一的定位到人的基因组上。5%为中间锚定read，对于这类，很多依赖于全局索引的算法就很难执行下去（需要比对很多次），而hisat，可以先将read中的长片段实现唯一比对，之后再使用局部索引对剩下的小片段进行比对（局部索引可以实现快速检索）。4.2%为短锚定read，因为这些序列特别短，因此只能通过在hisat比对其它read时发现的剪切位点或者用户自己提供的剪切位点来辅助比对。最后还有3%的是跨多个外显子的read，比对策略在hisat的online method中有介绍，文章中没有详解。比对过程中，中间锚定read、短锚定read、跨多个外显子read的比对占总比对时长的30%-60%，而且比对错误率很高！
三、软件安装
conda安装
代码语言：javascript复制conda install hisat2
四、hisat2 index建立
1.直接下载
直接在网站http://daehwankimlab.github.io/hisat2/download/#index下载已建好的index
代码语言：javascript复制wget -c https://genome-idx.s3.amazonaws.com/hisat/grch38_genome.tar.gz
2.自己构建index
hisat2可以自己建立index文件，如果没有现成的index的话，那就得自己去建立了，这个就比较麻烦而且耗内存和时间
代码语言：javascript复制# 其实hisat2-buld在运行的时候也会自己寻找exons和splice_sites，但是先做的目的是为了提高运行效率
extract_exons.py gencode.v26lift37.annotation.sorted.gtf > hg19.exons.gtf &
extract_splice_sites.py gencode.v26lift37.annotation.gtf > hg19.splice_sites.gtf &
# 建立index， 必须选项是基因组所在文件路径和输出的前缀
hisat2-build --ss hg19.splice_sites.gtf --exon hg19.exons.gtf genome/hg19/hg19.fa hg19
五、Hisat2的用法
安装完成以后，可以使用hisat2 -h来查看软件的帮助文档。
- 软件用法：
2.常用参数：
代码语言：javascript复制-x ：索引的前缀,也可以添加路径，例如 ~/genome/homo_GRCh38.dna.primary
-1 ：双端测序的第一个文件，若有多组，则用逗号隔开，名字与2要匹配，如-1 flyA_1.fq,flyB_1.fq
-2 ：类比上，名字与1要匹配，如-2 flyA_2.fq,flyB_2.fq
-U ：单端数据文件
-S ：保存到的sam文件，默认输出到标准输出
-p ：线程数
--dta ：在下游使用stringtie组装的时候量身定做的参数。使用此选项，
：HISAT2需要更长的锚长度才能从头发现拼接位点，这样可以减少与短锚的对齐，
： 这有助于转录汇编程序显着提高计算和内存使用率。
： 当然下游不使用stringtie也可以使用此参数减少短锚定read
--dta-cufflinks ：下游使用cufflinks则需要添加此参数
-q ：指定读取的测序文件是fastq文件（含有质量信息）
-f ：指定读取的测序文件是fa文件
-t ：打印加载索引文件和对齐读数所需的时钟时间
--phred33 ：质量表示方法，如果使用fq文件，建议选择，同理有小众的--phred64
--min-intronlen ：设置最小内含子的长度，默认值20
--max-intronlen ：设置最大内含子长度，默认500000
--known-splicesite-infile ：提供已知的剪接位点列表，HISAT2利用这些列表比对,可以用 hisat2_extract_splice_sites.py和gtf生成信息
--novel-splicesite-outfile ：在结果中生成一个剪接位点的列表
--novel-splicesite-infile ： 可以使用novel-splicesite-outfile所生成的剪接列表作为
： 新剪接列表，应该可以自己定义。为特定基因。
三、软件运行命令
hisat2 输出文件是sam格式，可通过管道符与Samtools工具连用，直接生成bam，并对bam文件进行sorted以方便后续数据处理
代码语言：javascript复制hisat2 -p 2 -x ${index} -1 $fq1 -2 $fq2 | samtools sort -@ 6 -o ./${sample}.bam -
批量运行脚本
代码语言：javascript复制mkdir 04.mapping
ls /home/data/lihe/lncRNA_project/03.trim/*_1.fq.gz > 1
ls /home/data/lihe/lncRNA_project/03.trim/*_2.fq.gz > 2
ls /home/data/lihe/lncRNA_project/03.trim/*_1.fq.gz | cut -d "/" -f 7 | cut -d "_" -f 1 > 0
paste 0 1 2 > config
cd 04.mapping
cat > 04.hg38_mapping.sh
index=/home/data/server/reference/index/hisat/hg38/genome
config=$1
number1=$2
number2=$3
cat $1 | while read id
do
if((i%$number1==$number2))
then
arr=(${id})
fq1=${arr[1]}
fq2=${arr[2]}
sample=${arr[0]}
hisat2 -p 2 -x ${index} -1 $fq1 -2 $fq2 | samtools sort -@ 6 -o ./${sample}.bam -
fi ## end for number1
i=$((i 1))
done
for i in {0..6}
do
(nohup bash 04.hg38_mapping.sh config 7 $i 1>log.$i.txt 2>&1 & )
done
四、结果解读
代码语言：javascript复制17968900 reads; of these: #reads数目
17968900 (100.00%) were paired; of these: #能配对的reads，由于用的是trimmomatic后的paired clean data，故100%
4434269 (24.68%) aligned concordantly 0 times #不能比对到基因组上
10748474 (59.82%) aligned concordantly exactly 1 time #比对到基因组上，且唯一匹配，即unique map
2786157 (15.51%) aligned concordantly >1 times #比对到基因组上，但结果>1次，即不唯一匹配
----
4434269 pairs aligned concordantly 0 times; of these: #上述不能比对上的reads中
269467 (6.08%) aligned discordantly 1 time #这些虽然比对上了，但是不合理，例如方向、R1 R2之间的距离等
----
4164802 pairs aligned 0 times concordantly or discordantly; of these: #剩余的不匹配reads中
8329604 mates make up the pairs; of these:
5311882 (63.77%) aligned 0 times
2073153 (24.89%) aligned exactly 1 time #作为单端来比对，这些可以唯一比对
944569 (11.34%) aligned >1 times #这些不唯一比对
85.22% overall alignment rate #结果是唯一比对 不唯一比对 不能比对中可以单端比对的
五、导入IGV可视化
使用samtools 建立bam文件的index，然后可把bam文件导入igv可视化
代码语言：javascript复制ls /home/data/lihe/lncRNA_project/04.mapping/*.bam > 1
ls /home/data/lihe/lncRNA_project/04.mapping/*.bam | cut -d "/" -f 7 | cut -d "_" -f 1 | cut -d "." -f 1 > 0
paste 0 1 > config
cat > index.sh
config=$1
number1=$2
number2=$3
cat $1 | while read id
do
if((i%$number1==$number2))
then
arr=(${id})
input=${arr[1]}
sample=${arr[0]}
samtools index -@ 4 ${input} ./${sample}.bam.bai
fi ## end for number1
i=$((i 1))
done
for i in {0..4}
do
(nohup bash index.sh config 5 $i 1>log.$i.txt 2>&1 & )
done
