一、 基本介紹

HISAT2 是一款利用改進(jìn)的BWT算法進(jìn)行序列比對(duì)的軟件。由約翰霍普金斯大學(xué)計(jì)算生物學(xué)中心（CCB at JHU）開(kāi)發(fā)，是TopHat的升級(jí)版本，速度提高了50倍。利用 HISAT2 + StringTie 流程，可以快速地分析轉(zhuǎn)錄組測(cè)序數(shù)據(jù)，獲得每個(gè)基因和轉(zhuǎn)錄本的表達(dá)量。

• TopHat采用外顯子優(yōu)先的方法，第一步使用BWT方法將完整的匹配外顯子的reads優(yōu)先回貼到參考基因組，第二步將沒(méi)有比對(duì)到基因組的reads切割成小片段，分別比對(duì)到基因組。
• HISAT2使用了疊加的比對(duì)算法，全局比對(duì)整個(gè)基因組建立index，輔助成千上萬(wàn)個(gè)小的局部index。
• STAR采用了種子和擴(kuò)展的方法。

首先需要構(gòu)建參考基因組索引用于下一步的比對(duì)。HISAT2提供了兩個(gè)腳本用于從基因組注釋GTF文件中提取剪接位點(diǎn)和外顯子位置，基于這些特征，可以使 RNA-Seq reads 比對(duì)更加準(zhǔn)確。建議直接在HISAT2官網(wǎng)下載構(gòu)建好的index。（為什么需要索引：為了提高比對(duì)速度，需要根據(jù)參考基因組序列，經(jīng)過(guò)BWT算法轉(zhuǎn)換成index，而我們比對(duì)的序列其實(shí)是index的一個(gè)子集。因此不是直接把reads回貼到基因組上，而是把reads與index進(jìn)行比較。）

然后再進(jìn)行reads mapping。reads比對(duì)到基因組上的一個(gè)位置稱為一個(gè)alignment。軟件會(huì)對(duì)所有的alignments進(jìn)行打分和判斷，能有符合過(guò)濾條件的valid alignment才會(huì)輸出。打分機(jī)制包括：錯(cuò)配堿基罰分（--mp），reads上的gap罰分（--rdg），reference上的gap罰分（--rdg），valid alignment分?jǐn)?shù)閾值（--score—min，默認(rèn)L,0,-0.2），多個(gè)valid alignments輸出數(shù)目（-k，默認(rèn)5）。

二、 背景知識(shí)

(1) BWT索引

處理NGS測(cè)序數(shù)據(jù)最常見(jiàn)的任務(wù)（task）是讀段映射（read mapping）：在參考基因組序列（reference genomic sequence）中定位讀段（locating reads）并計(jì)算相應(yīng)的比對(duì)（alignments）。在這里，輸入的reads通常來(lái)自正在研究的有機(jī)體（例如患病個(gè)體），參考基因組是已知的具有代表性的物種基因組序列（例如參考人類基因組序列），或可能是系統(tǒng)發(fā)育相關(guān)的物種（phylogenetically related species）。一個(gè)主要的困難來(lái)自數(shù)據(jù)的大?。簲?shù)以百萬(wàn)計(jì)的reads必須與一個(gè)可以包含數(shù)十億個(gè)字母的序列進(jìn)行比對(duì)。類似BLAST的工具根本無(wú)法在合理的時(shí)間內(nèi)處理這樣的數(shù)據(jù)。

開(kāi)發(fā)出了各種各樣的read mapping軟件。最常見(jiàn)的實(shí)用映射軟件包括BWA-MEM、Bowtie2、GEM，它們比類似BLAST的工具快幾個(gè)數(shù)量級(jí)。這些算法的一個(gè)主要效率（efficiency）來(lái)源是被稱為BWT-或FM-index的索引結(jié)構(gòu)（indexing structure）。

BWT-index的根源在于詞組合（word combinatorics）：它基于Burrows-Wheeler變換（Burrows-Wheeler transform），從組合的角度來(lái)看，它與Gessel-Reutenauer雙射（bijection，一一對(duì)應(yīng)關(guān)系）有著密切的聯(lián)系。雖然BWT的主要應(yīng)用是文本壓縮（text compression），但它被應(yīng)用于構(gòu)造一個(gè)簡(jiǎn)潔的文本索引（compact text index），結(jié)果令人驚訝地強(qiáng)大。

與以BLAST為代表的基于哈希的索引基礎(chǔ)（hash-based indexes underpinning）支持種子和擴(kuò)展策略（seed-and-extend strategies）不同，BWT索引是一個(gè)全文索引（full-text index），因?yàn)樗С炙阉魅我獯笮〉淖址╝rbitrary-size strings）。BWT-index屬于一個(gè)更普遍的簡(jiǎn)潔索引（succinct indexes）系列，其以位（bits）為單位的大小接近于對(duì)象（這里是序列）的無(wú)損表示（lossless representation of the object）所需的信息理論最小值（information-theoretic minimum）。
——Evolution of biosequence search algorithms: a brief survey

(2) SAM和BAM文件

SAM（Sequence Alignment/Mapping）格式是一種通用的比對(duì)格式，用來(lái)存儲(chǔ)reads到參考序列的比對(duì)信息（是目前高通量測(cè)序中存放比對(duì)數(shù)據(jù)的標(biāo)準(zhǔn)格式）。BAM是SAM的二進(jìn)制格式。bam文件優(yōu)點(diǎn)：bam文件為二進(jìn)制文件，占用的磁盤空間比sam文本文件??；利用bam二進(jìn)制文件的運(yùn)算速度快。

SAM分為兩部分，標(biāo)頭注釋信息（header section）和比對(duì)結(jié)果部分（alignment section）。

• 標(biāo)頭信息：標(biāo)頭信息可有可無(wú)，都是以@開(kāi)頭，用不同的tag表示不同的信息。

@HD：符合標(biāo)準(zhǔn)的版本、對(duì)比序列的排列順序
@SQ：參考序列說(shuō)明
@RG：比對(duì)上的 reads 說(shuō)明
@PG：使用說(shuō)明
@Co：任意的說(shuō)明信息

• 比對(duì)結(jié)果部分：除注釋外，每一行是一個(gè)read，包括11個(gè)必須的字段（mandatory fields）和一個(gè)可選的字段，字段之間用tab分割。必須字段的順序固定，根據(jù)字段定義，可以為0或者*。

第1列：read的名稱（read表示本條read，mate表示pair中的另一條read）。
第2列：flag數(shù)字之和。（見(jiàn)下文）
第3列：比對(duì)到參考序列上的染色體號(hào)。若無(wú)法比對(duì)則是*。
第4列：read比對(duì)到參考序列上第一個(gè)堿基所在的位置。若無(wú)法比對(duì)則是0。
第5列：比對(duì)的質(zhì)量分?jǐn)?shù)，比對(duì)錯(cuò)誤率的-10log10值，越高說(shuō)明該read比對(duì)到參考基因組上的位置越唯一。（見(jiàn)下文）
第6列：CIGAR值，read比對(duì)的具體情況。（見(jiàn)下文）
第7列：mate比對(duì)到的染色體號(hào)，若沒(méi)有mate則是*。
第8列：mate比對(duì)到參考序列上的第一個(gè)堿基位置，若無(wú)mate則為0。
第9列：來(lái)自一個(gè)fragment的兩個(gè)reads覆蓋的堿基數(shù)，若無(wú)mate則為0。
第10列：read的堿基序列，如果是比對(duì)到互補(bǔ)鏈上則是reverse completed。
第11列：read質(zhì)量的ASCII編碼。

第2列：flag

0：read是Single-read且正向比對(duì)
1：read是pair中的一條
2：pair一正一負(fù)完美的比對(duì)上
4：這條read沒(méi)有比對(duì)上
8：mate沒(méi)有比對(duì)上
16：這條read反向比對(duì)
32：mate反向比對(duì)
64：這條read是read1
128：這條read是read2
256：secondary比對(duì)
512：比對(duì)質(zhì)量不合格
1024：read是PCR或光學(xué)copy產(chǎn)生
2048：輔助比對(duì)結(jié)果

第6列：CIGAR值

M表示 match或 mismatch
I表示 insert
D表示 deletion
N表示 skipped（跳過(guò)這段區(qū)域）
S表示 soft clipping（被剪切的序列存在于序列中）
H表示 hard clipping（被剪切的序列不存在于序列中）
P表示 padding
=表示 match
X表示 mismatch（錯(cuò)配，位置是一一對(duì)應(yīng)的）

第5列：
實(shí)際上hisat2比對(duì)的結(jié)果MAPQ沒(méi)有實(shí)際意義，不表示-10log10比對(duì)錯(cuò)誤率，它只有0、1、60三個(gè)值。

• 60：uniquely mapped read, regardless of number of mismatches / indels
• 1：multiply mapped, perfect match or few mismatches / indels
• 0：unmapped, or multiply mapped and with lots of mismatches / indels

(3) GTF格式

提供基因位置的注釋文件通常以GTF（General Transfer Format）或GFF3（General Feature Format）格式呈現(xiàn)。有GTF文件后，就可以利用注釋信息計(jì)算每個(gè)基因/轉(zhuǎn)錄本/外顯子比對(duì)了多少reads，從而獲取counts值。GTF是GFF的便于傳輸版。

GTF格式各列的內(nèi)容為：
seqname - 染色體或scaffold的名稱。
source - 生成這個(gè)特征的項(xiàng)目名稱，或數(shù)據(jù)庫(kù)來(lái)源。
feature - 特征類型名稱，如gene、transcript、exon、CDS。
start - 開(kāi)始位置，使用基于1的索引。
end - 結(jié)束位置，使用基于1的索引。
score - 組裝的轉(zhuǎn)錄本的可信度分?jǐn)?shù)（目前該字段未被StringTie使用，如果轉(zhuǎn)錄本與read alignment bundle有連接，StringTie將報(bào)告一個(gè)常量值1000）。
strand - 正鏈或負(fù)鏈+/-。
frame - 密碼子的第幾個(gè)堿基0/1/2（StringTie不使用這個(gè)字段，只記錄一個(gè)點(diǎn).）。
attribute - 附加信息。

例如：
a. 來(lái)自ensembl的果蠅的Drosophila_melanogaster.BDGP6.32.109.gtf。

feature 的內(nèi)容是gene、transcript、exon、start_codon、stop_codon、five_prime_utr、three_prime_utr、CDS、Selenocysteine。
attributes的內(nèi)容例如：gene_id "FBgn0250732"; transcript_id "FBtr0091512"; gene_name "gfzf"; gene_source "FlyBase"; gene_biotype "protein_coding"; transcript_name "gfzf-RB"; transcript_source "FlyBase"; transcript_biotype "protein_coding"; tag "Ensembl_canonical";

b. 來(lái)源UCSC的人類GTF文件，1號(hào)染色體是 chr1。

c. 來(lái)源Ensembl的Homo_sapiens.GRCh38.chr.gtf.gz，1號(hào)染色體是1。

d. Stringtie得到的sample.gtf。

GFF格式各列的內(nèi)容為：
seqid - 染色體或scaffold的名稱。
source - 生成這個(gè)特征的項(xiàng)目名稱，或數(shù)據(jù)庫(kù)來(lái)源。
feature - 特征類型名稱，來(lái)自SOFA sequence ontology。
start
end
score
strand - 正鏈或負(fù)鏈+/-。
phase - 密碼子的第幾個(gè)堿基0/1/2。
attribute - 附加信息。A semicolon-separated list of tag-value pairs。

GTF和GFF之間的區(qū)別：
數(shù)據(jù)結(jié)構(gòu)：都是由9列構(gòu)成，分別是reference sequence name; annotation source; feature type; start coordinate; end coordinate; score; strand; frame; attributes.前8列都是相同的，第9列不同。
GFF第9列：都是以鍵值對(duì)的形式，鍵值之間用“=”連接，不同屬性之間用“；”分隔，都是以ID這個(gè)屬性開(kāi)始。
GTF第9列：同樣以鍵值對(duì)的形式，鍵值之間是以空格區(qū)分，值用雙引號(hào)括起來(lái)；不同屬性之間用“；”分隔；開(kāi)頭必須是gene_id, transcript_id兩個(gè)屬性。

三、 hisat2使用方法

(1) 用法和參數(shù)

hisat2  [options]*  -x <ht2-idx>  {-1 <m1>  -2 <m2> |  -U <r>}  [-S <sam>]

#提取外顯子和剪接位點(diǎn)
extract_splice_sites.py  dmel.annotation.gtf  >  dmel.ss 
extract_exons.py  dmel.annotation.gtf  >  dmel.exon 
extract_snps.py  snp142Common.txt  >  genome.snp

#構(gòu)建參考基因組索引
hisat2-build  --ss  dmel.ss  --exon  dmel.exon  (--snp  genome.snp)  dmel.genome.fa  dmel  #得到8個(gè)索引文件，以dmel為共同的前綴，dmel.1~8.ht2

# Reads mapping
hisat2  -p 4  --dta  -x dmel  -1 sample_1.fq.gz  -2 sample_2.fq.gz  -S sample.sam

{-1 <m1> -2 <m2> | -U <r>}：對(duì)于單端數(shù)據(jù)，采用-U指定輸入文件（-U <r>）；對(duì)于雙端數(shù)據(jù)，采用-1和-2分別指定R1端和R2端的輸入文件（-1 <m1> -2 <m2>）。
-x：指定基因組索引的文件前綴（前綴.1~8.ht2）
-S：指定輸出的SAM文件<sam>。默認(rèn)輸出到標(biāo)準(zhǔn)輸出，比對(duì)結(jié)束后統(tǒng)計(jì)結(jié)果輸出到標(biāo)準(zhǔn)錯(cuò)誤輸出。
-p：線程數(shù)。
--dta：輸出專門為轉(zhuǎn)錄本組裝的比對(duì)結(jié)果。Reports alignments tailored for transcript assemblers.（如果比對(duì)結(jié)果后續(xù)需要使用StringTie進(jìn)行轉(zhuǎn)錄本組裝，則一定要加入--dta選項(xiàng)。）
--rna-strandness：默認(rèn)unstranded，如果是鏈特異性文庫(kù)要加上--rna-strandness RF。

-f：表示輸入文件格式為fasta（默認(rèn)），-q表示輸入文件格式為fastq?？梢允莋zip壓縮的文件。
-phred33：輸入的FASTQ文件堿基質(zhì)量值編碼標(biāo)準(zhǔn)為phred33（默認(rèn)）。
-k <int>：report up to <int> alignments per read; MAPQ not meaningful.

(2) 舉個(gè)例子

#提取外顯子和剪接位點(diǎn)
extract_splice_sites.py  annotation/genome/Drosophila_melanogaster.BDGP6.32.105.gtf  >  annotation/genome/dmel.ss
extract_exons.py  annotation/genome/Drosophila_melanogaster.BDGP6.32.105.gtf  >  annotation/genome/dmel.exon

#構(gòu)建參考基因組索引
hisat2-build  --ss  annotation/genome/dmel.ss  --exon  annotation/genome/dmel.exon  annotation/genome/Drosophila_melanogaster.BDGP6.32.dna.toplevel.fa  ./annotation/index/dmel  2>>log/hisat2_index_log.txt  &  #在annotation/index/文件夾下生成了8個(gè)索引文件

# Reads mapping
hisat2  -p  4  --dta  --rna-strandness RF  -x  ./annotation/index/dmel  -1  02_cleandata/WTF4.Input_1_val_1.fq.gz  -2  02_cleandata/WTF4.Input_2_val_2.fq.gz  -S 04_bam/WTF4.Input.sam  2>>log/hisat2_log.txt  #雙端，鏈特異性文庫(kù)

# SAM轉(zhuǎn)為BAM
samtools sort  -@ 4  -o 04_bam/WTF4.Input.sorted.bam  04_bam/WTF4.Input.sam

(3) 批量處理

#以雙端為例：
for i in `ls 02_cleandata/*_1_val_1.fq.gz `
do

    i=`basename $i`
    i=${i%*_1_val_1.fq.gz}
    echo $i >>log/hisat2_log.txt

hisat2  -p  4  --dta  --rna-strandness RF  -x  ./annotation/index/dmel  -1  02_cleandata/$i\_1_val_1.fq.gz  -2  02_cleandata/$i\_2_val_2.fq.gz  2>>log/hisat2_log.txt | samtools sort  -@  4  -o  04_bam/$i.bam   2>>log/hisat2_log.txt

done

(4) 質(zhì)控和建立索引

#批量對(duì)bam文件進(jìn)行QC
ls 04_bam/*.bam | while read id ; do (samtools flagstat $id > 04_bam/$(basename $id ".bam").stat) ; done  &
#每一個(gè)bam文件都生成了一個(gè)stat文件
274 + 0 in total (QC-passed reads + QC-failed reads)
50 + 0 primary
224 + 0 secondary
0 + 0 supplementary
0 + 0 duplicates
0 + 0 primary duplicates
270 + 0 mapped (98.54% : N/A)
46 + 0 primary mapped (92.00% : N/A)
50 + 0 paired in sequencing
25 + 0 read1
25 + 0 read2
46 + 0 properly paired (92.00% : N/A)
46 + 0 with itself and mate mapped
0 + 0 singletons (0.00% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

#批量對(duì)bam文件建立索引
for i in `ls 04_bam/*.bam` ; do (samtools index $i) ; done  &
#每一個(gè)bam文件都生成了一個(gè)bam.bai文件

(5) 結(jié)果解讀

通過(guò)2>>log將比對(duì)率等信息的標(biāo)準(zhǔn)錯(cuò)誤輸出重定向到log文件。

查看SAM結(jié)果：less 04_bam/WTF4.Input.sam

查看BAM結(jié)果：samtools view 04_bam/WTF4.Input.sorted.bam | less

四、 samtools使用方法

(1) 將sam文件轉(zhuǎn)化為bam文件

samtools view  [options]  <in.bam>|<in.sam>  [region ...]
samtools view  -bS  sample.sam  >  sample_unsorted.bam

samtools view  -bF 4  sample.bam  >  sample.F4.bam
samtools view  -b  -F 8  -f 4  sample.sam  >  only.read2.mapped.bam
samtools view  -@ 4  -h  -b  -q 30  -F 4  -F 256  example.bam  >  example.filtered.bam
samtools view  example.bam  scaffold1  >  scaffold1.sam
samtools view  example.bam  scaffold1:30000-100000  >  scaffold1_30k-100k.sam

-b：輸出為bam格式，默認(rèn)輸出SAM格式。
-S：出入為sam格式，默認(rèn)輸入文件是bam文件。新版本自動(dòng)檢測(cè)，可忽略該參數(shù)。
-o：輸出文件的名字。
-h：設(shè)定輸出sam文件時(shí)帶header信息，默認(rèn)輸出sam不帶header信息。
-H：僅輸出頭文件信息。
-@：線程數(shù)。
-F：過(guò)濾掉的flag。例如，數(shù)字4代表該序列沒(méi)有比對(duì)到參考序列上；數(shù)字8代表該序列的mate序列沒(méi)有比對(duì)到參考序列上。-F 4表示比對(duì)到參考序列上的結(jié)果，-F 12表示雙端都比對(duì)到參考序列的結(jié)果。
-f：需要的flag。-f 4表示沒(méi)有比對(duì)到參考序列上的比對(duì)結(jié)果。
-q：have mapping quality >= INT
region：提取比對(duì)到指定區(qū)域上的比對(duì)結(jié)果

(2) 將bam文件進(jìn)行排序

samtools sort  [options...]  [in.bam]
samtools sort  -@ 4  sample_unsorted.bam  -o sample.sorted.bam

#版本1.4后可簡(jiǎn)化直接由sam到sorted.bam
samtools sort  -@ 4  -o sample.sorted.bam  sample.sam

-o：輸出文件的名字。
-O：指定輸出格式 sam/bam/cram。默認(rèn)基于-o文件的擴(kuò)展名選擇一個(gè)格式，如.bam。如果格式無(wú)法推測(cè)，則需要指定-O。
-n/-t：默認(rèn)按照染色體位置（最左邊的坐標(biāo)）進(jìn)行排序。-n是根據(jù)read名（QNAME）進(jìn)行排序，-t 根據(jù)TAG進(jìn)行排序。
-@：線程數(shù)。
-m：內(nèi)存參數(shù)，默認(rèn)是500,000,000，即500M。

(3) 索引

samtools index  <in.bam>  [out.index]
samtools index  sample.sorted.bam
for i in `ls 04_bam/*.bam` ; do (samtools index $i) ; done  &

得到sample.sorted.bam.bai索引文件。IGV需要bam文件對(duì)應(yīng)的.bai索引文件。

(4) 過(guò)濾

samtools rmdup  [-sS]  <input.srt.bam>  <output.bam>
samtools rmdup  sample.sorted.bam  sample.rmdup.bam

-s：rmdup for SE reads
-S：treat PE reads as SE in rmdup (force -s)

• rmdup已經(jīng)過(guò)時(shí)了，新版使用markdup，效果和picard類似。
  按read name排序：samtools sort -n xxx.sort.bam -o xxx.sortname.bam
  然后：samtools fixmate -m xxx.sortname.bam xxx.fixmate.bam
  按position排序：samtools sort xxx.fixmate.bam -o xxx.sortposition.bam
  最后：samtools markdup -r xxx.sortposition.bam xxx.markdup.bam (-r: 除去重復(fù)reads)

(5) 合并

samtools merge  [options]  -o <out.bam>  [options]  <in1.bam> ... <inN.bam>
samtools merge  [options]  <out.bam>  <in1.bam> ... <inN.bam>
samtools merge  -h test.sam  -R chr7  test_L1_L2.bam  test_L1.sorted.bam  test_L2.sroted.bam

可以將2個(gè)或2個(gè)以上的已經(jīng)sort了的bam文件融合成一個(gè)bam文件。
-h：指定文件內(nèi)的頭文件復(fù)制并替換輸出文件的頭文件，默認(rèn)從第一個(gè)輸入文件復(fù)制過(guò)來(lái)。
-R：合并輸出文件的指定區(qū)域。

(6) faidx

samtools faidx  <file.fa|file.fa.gz>  [<reg> [...]] 
samtools faidx  genome.fasta
samtools faidx  genome.fasta  scffold_10  >  scaffold_10.fasta

對(duì)fasta文件建立索引，生成的索引文件以.fai后綴結(jié)尾。由于有索引文件，可以很快從基因組中提取到fasta格式的子序列。

(7) flagstat

samtools flagstat  [options]  <in.bam>
samtools flagstat  example.bam
ls 04_bam/*.bam | while read id ; do (samtools flagstat $id > 04_bam/$(basename $id ".bam").stat) ; done  &

給出BAM文件的比對(duì)結(jié)果：總reads數(shù)，總reads匹配率，多少paired reads，reads1中的reads數(shù)，reads2中的reads數(shù)，完美匹配的reads數(shù)，paired reads中兩條都比對(duì)到參考序列上的reads數(shù)，單獨(dú)一條匹配到參考序列上的reads數(shù)（和上一個(gè)相加是總的匹配上的reads數(shù)），paired reads中兩條分別比對(duì)到兩條不同的參考序列的reads數(shù)。

(8) depth

samtools depth [options] in.bam [in.bam ...]
samtools depth  example.bam  >  depth

得到每個(gè)堿基位點(diǎn)的測(cè)序深度，并輸出到標(biāo)準(zhǔn)輸出。需要index。

(9) mpileup

samtools mpileup [options] in1.bam [in2.bam [...]]
samtools mpileup  -f  genome.fasta  example.bam  >  example.txt
samtools mpileup  -gSDf  genome.fasta  example.bam  >  example.bcf

該命令用于生成bcf文件，再使用bcftools進(jìn)行SNP和Indel的分析。
mpileup不使用-u或-g參數(shù)時(shí)，不生成二進(jìn)制的bcf文件，而生成一個(gè)文本文件，輸出到標(biāo)準(zhǔn)輸出。該文本文件統(tǒng)計(jì)了參考序列中每個(gè)堿基位點(diǎn)的比對(duì)情況，每一行代表了參考序列中某一個(gè)堿基位點(diǎn)的比對(duì)結(jié)果。
mpileup生成的結(jié)果包含6行：參考序列名；位置；參考?jí)A基；比對(duì)上的reads數(shù)；比對(duì)情況；比對(duì)上的堿基的質(zhì)量。
-g：計(jì)算基因型的似然值和輸出文件格式為BCF。
-f：輸入有索引文件的fasta參考序列。
-D：輸出每個(gè)樣本的reads深度。