一、比對(duì)原理
原文：HISAT: a fast spliced aligner with low memory requirements.
摘要
1、Hisat是一種高效的RNA-seq實(shí)驗(yàn)比對(duì)工具
2、它使用了基于BWT和Ferragina-manzini (Fm) index 兩種算法的索引框架
3、使用了兩類索引去比對(duì)，一類是全基因組范圍的FM索引來錨定每一個(gè)比對(duì)，另一類是大量的局部索引對(duì)這些比對(duì)做快速的擴(kuò)展
背景
1、自2008年起，RNAseq已經(jīng)成為研究基因表達(dá)、轉(zhuǎn)錄本結(jié)構(gòu)、長鏈非編碼RNA確定以及融合轉(zhuǎn)錄本的重要手段
2、隨著測序深度的加深和read讀長的延長，給比對(duì)工作帶來很多困難
3、目前的工具如Tophat2和GSNAP等在對(duì)單個(gè)轉(zhuǎn)錄組實(shí)驗(yàn)的比對(duì)中耗時(shí)幾天，而新型的STAR雖然很快，但是會(huì)吃掉大量的內(nèi)存空間，如基于人類基因組的比對(duì)需要消耗28G左右，而hisat2是4.3G
Hisat設(shè)計(jì)原則
優(yōu)化了索引建立策略
1、hisat應(yīng)用了基于bowtie2的方法去處理很多低水平的用于構(gòu)建和查詢FM索引的操作
2、但是與其它比對(duì)器不同的是，該軟件應(yīng)用了兩類不同的索引類型：代表全基因組的全局FM索引和大量的局部小索引，每個(gè)索引代表64000bp
3、以人類基因組為例，創(chuàng)建了48000個(gè)局部索引，每一個(gè)覆蓋1024bp，最終可以覆蓋這個(gè)3 billion 的堿基的基因組。這種存在交叉（overlap）的邊界可以輕松的比對(duì)那些跨區(qū)域的read（可變剪切體）
4、盡管有很多索引，但是hisat會(huì)把他們使用合適方法壓縮，最終只占4gb左右的內(nèi)存
采用了新的比對(duì)策略
1、RNA-seq產(chǎn)生的reads可能跨長度比較大的內(nèi)含子，哺乳動(dòng)物中甚至最長能達(dá)到1MB，同時(shí)外顯子比較短，read也比較短，會(huì)有很多read（模擬數(shù)據(jù)中大概34%）跨兩個(gè)外顯子的情況
2、為了更好的比對(duì)，將跨外顯子的reads分成了三類：1）長錨定read，至少有16bp在兩個(gè)外顯子的每一個(gè)上 2）中間錨定read，有8-15bp在一個(gè)外顯子上 3）短錨定read，只有1-7bp在一個(gè)外顯子上
3、所以總的reads可以被劃分為五類：1）不跨外顯子的read 2）長錨定read 3）中間錨定read 4）短錨定read 5）跨兩個(gè)外顯子以上的read
4、在模擬的數(shù)據(jù)中，有25%左右的read是長錨定read，這種read在大多數(shù)情況下可以被唯一的定位到人的基因組上
5、5%為中間錨定read，對(duì)于這類，很多依賴于全局索引的算法就很難執(zhí)行下去（需要比對(duì)很多次），而hisat，可以先將read中的長片段實(shí)現(xiàn)唯一比對(duì)，之后再使用局部索引對(duì)剩下的小片段進(jìn)行比對(duì)（局部索引可以實(shí)現(xiàn)快速檢索）
6、4.2%為短錨定read，因?yàn)檫@些序列特別短，因此只能通過在hisat比對(duì)其它read時(shí)發(fā)現(xiàn)的剪切位點(diǎn)或者用戶自己提供的剪切位點(diǎn)來輔助比對(duì)
7、最后還有3%的是跨多個(gè)外顯子的read，比對(duì)策略在hisat的online method中有介紹，文章中沒有詳解
8、比對(duì)過程中，中間錨定read、短錨定read、跨多個(gè)外顯子read的比對(duì)占總比對(duì)時(shí)長的30%-60%，而且比對(duì)錯(cuò)誤率很高！
二、比對(duì)方案
說明書：http://ccb.jhu.edu/software/hisat2/manual.shtml
網(wǎng)頁版如何轉(zhuǎn)化為pdf：https://jingyan.baidu.com/article/6b97984df464561ca2b0bfb9.html
基本功能
建索引：hisat2-build [options]* <reference_in> <ht2_base>
#<reference_in> ：fasta文件 list，如果為list，使用逗號(hào)分開
#<ht2_base> ：索引文件的前綴名，如設(shè)為xxx，則生成的索引文件為xxx.1.ht2,xxx.2.ht2，默認(rèn)的前綴名為NAME
#option:詳見說明書
檢查索引：hisat2-inspect [options]* <ht2_base>
#輸出結(jié)果為一個(gè)fasta文件，主要用于檢查已經(jīng)構(gòu)建好的索引所用的構(gòu)建信息，感覺沒啥用
比對(duì)：hisat2 [options]* -x <hisat2-idx> {-1 <m1> -2 <m2> | -U <r> | --sra-acc <SRA accession number>} [-S <hit>]
#參數(shù)說明：
#-p ：線程數(shù)目
#--dta ：注意！?。≡谙掠问褂胹tringtie組裝的時(shí)候一定要在hisat中設(shè)置這個(gè)參數(shù)?。?！
#-x <hisat2-idx> ：參考基因組索引的basename，即前綴名
#{}：其中的內(nèi)容意思為hisat2可以接受單端測序，雙端測序，或者直接提交SRA ID號(hào)
#-1 <m1> ：雙端測序的read1 list ，若為list，使用逗號(hào)隔開，名字與2要匹配，如-1 flyA_1.fq,flyB_1.fq
#-2 <m2> ：雙端測序的read2 list ，若為list，使用逗號(hào)隔開，名字與1要匹配，如-2 flyA_2.fq,flyB_2.fq
#-U <r>：單端測序list，若為list，使用逗號(hào)隔開，-U lane1.fq,lane2.fq,lane3.fq,lane4.fq
#--sra-acc <SRA accession number> : SRAID list，若為list，使用逗號(hào)隔開，--sra-acc SRR353653,SRR353654
#-S <hit> ：SAM寫入的文件名，默認(rèn)寫入到標(biāo)準(zhǔn)輸出中
#options:這里只列出可調(diào)節(jié)的類別，至于參數(shù)調(diào)整，詳見說明書
#Input options
#Alignment options
#Scoring options
#Spliced alignment options（重要）
#Reporting options
#Paired-end options（重要）
#Output options（重要）
#SAM options
#Performance options
#Other options
report格式
若為單端測序
20000 reads; of these:
20000 (100.00%) were unpaired; of these:
1247 (6.24%) aligned 0 times
18739 (93.69%) aligned exactly 1 time
14 (0.07%) aligned >1 times
93.77% overall alignment rate
若為雙端測序
10000 reads; of these:
10000 (100.00%) were paired; of these:
650 (6.50%) aligned concordantly 0 times
8823 (88.23%) aligned concordantly exactly 1 time
527 (5.27%) aligned concordantly >1 times
----
650 pairs aligned concordantly 0 times; of these:
34 (5.23%) aligned discordantly 1 time
----
616 pairs aligned 0 times concordantly or discordantly; of these:
1232 mates make up the pairs; of these:
660 (53.57%) aligned 0 times
571 (46.35%) aligned exactly 1 time
1 (0.08%) aligned >1 times
96.70% overall alignment rate
三、具體實(shí)例
hisat2-build human.fa --snp human.snp human_hg38
#最終會(huì)生成以human_hg38為前綴，以.*.ht2為后綴的索引文件
#其中提供了snp信息
單端測序比對(duì)
hisat2 -f -x human_hg38 -U reads_1.fa -S eg1.sam
#-f:輸入文件為fasta文件
#-q：輸入文件為fastq文件
#-x：只寫前綴
#最終生成eg1.sam文件
雙端測序比對(duì)
hisat2 -f -x human_hg38 -1 reads_1.fa -2 reads_2.fa -S eg2.sam
#最后生成eg2.sam
可以使用samtools或者bcftools作為下游
samtools view -bS eg2.sam > eg2.bam
#轉(zhuǎn)化為bam文件
samtools sort eg2.bam -o eg2.sorted.bam
#bam文件轉(zhuǎn)變?yōu)閟orted bam
samtools mpileup -uf human.fa eg2.sorted.bam | bcftools view -bvcg - > eg2.raw.bcf
#生成bcf文件
bcftools view eg2.raw.bcf
#查看bcf文件
#詳見samtools工具的使用
