一、 基本介绍

HISAT2 是一款利用改进的BWT算法进行序列比对的软件。由约翰霍普金斯大学计算生物学中心（CCB at JHU）开发，是TopHat的升级版本，速度提高了50倍。利用 HISAT2 + StringTie 流程，可以快速地分析转录组测序数据，获得每个基因和转录本的表达量。

• TopHat采用外显子优先的方法，第一步使用BWT方法将完整的匹配外显子的reads优先回贴到参考基因组，第二步将没有比对到基因组的reads切割成小片段，分别比对到基因组。
• HISAT2使用了叠加的比对算法，全局比对整个基因组建立index，辅助成千上万个小的局部index。
• STAR采用了种子和扩展的方法。

首先需要构建参考基因组索引用于下一步的比对。HISAT2提供了两个脚本用于从基因组注释GTF文件中提取剪接位点和外显子位置，基于这些特征，可以使 RNA-Seq reads 比对更加准确。建议直接在HISAT2官网下载构建好的index。（为什么需要索引：为了提高比对速度，需要根据参考基因组序列，经过BWT算法转换成index，而我们比对的序列其实是index的一个子集。因此不是直接把reads回贴到基因组上，而是把reads与index进行比较。）

然后再进行reads mapping。reads比对到基因组上的一个位置称为一个alignment。软件会对所有的alignments进行打分和判断，能有符合过滤条件的valid alignment才会输出。打分机制包括：错配碱基罚分（--mp），reads上的gap罚分（--rdg），reference上的gap罚分（--rdg），valid alignment分数阈值（--score—min，默认L,0,-0.2），多个valid alignments输出数目（-k，默认5）。

二、 背景知识

(1) BWT索引

处理NGS测序数据最常见的任务（task）是读段映射（read mapping）：在参考基因组序列（reference genomic sequence）中定位读段（locating reads）并计算相应的比对（alignments）。在这里，输入的reads通常来自正在研究的有机体（例如患病个体），参考基因组是已知的具有代表性的物种基因组序列（例如参考人类基因组序列），或可能是系统发育相关的物种（phylogenetically related species）。一个主要的困难来自数据的大小：数以百万计的reads必须与一个可以包含数十亿个字母的序列进行比对。类似BLAST的工具根本无法在合理的时间内处理这样的数据。

开发出了各种各样的read mapping软件。最常见的实用映射软件包括BWA-MEM、Bowtie2、GEM，它们比类似BLAST的工具快几个数量级。这些算法的一个主要效率（efficiency）来源是被称为BWT-或FM-index的索引结构（indexing structure）。

BWT-index的根源在于词组合（word combinatorics）：它基于Burrows-Wheeler变换（Burrows-Wheeler transform），从组合的角度来看，它与Gessel-Reutenauer双射（bijection，一一对应关系）有着密切的联系。虽然BWT的主要应用是文本压缩（text compression），但它被应用于构造一个简洁的文本索引（compact text index），结果令人惊讶地强大。

与以BLAST为代表的基于哈希的索引基础（hash-based indexes underpinning）支持种子和扩展策略（seed-and-extend strategies）不同，BWT索引是一个全文索引（full-text index），因为它支持搜索任意大小的字符串（arbitrary-size strings）。BWT-index属于一个更普遍的简洁索引（succinct indexes）系列，其以位（bits）为单位的大小接近于对象（这里是序列）的无损表示（lossless representation of the object）所需的信息理论最小值（information-theoretic minimum）。
——Evolution of biosequence search algorithms: a brief survey

(2) SAM和BAM文件

SAM（Sequence Alignment/Mapping）格式是一种通用的比对格式，用来存储reads到参考序列的比对信息（是目前高通量测序中存放比对数据的标准格式）。BAM是SAM的二进制格式。bam文件优点：bam文件为二进制文件，占用的磁盘空间比sam文本文件小；利用bam二进制文件的运算速度快。

SAM分为两部分，标头注释信息（header section）和比对结果部分（alignment section）。

• 标头信息：标头信息可有可无，都是以@开头，用不同的tag表示不同的信息。

@HD：符合标准的版本、对比序列的排列顺序
@SQ：参考序列说明
@RG：比对上的 reads 说明
@PG：使用说明
@Co：任意的说明信息

• 比对结果部分：除注释外，每一行是一个read，包括11个必须的字段（mandatory fields）和一个可选的字段，字段之间用tab分割。必须字段的顺序固定，根据字段定义，可以为0或者*。

第1列：read的名称（read表示本条read，mate表示pair中的另一条read）。
第2列：flag数字之和。（见下文）
第3列：比对到参考序列上的染色体号。若无法比对则是*。
第4列：read比对到参考序列上第一个碱基所在的位置。若无法比对则是0。
第5列：比对的质量分数，比对错误率的-10log10值，越高说明该read比对到参考基因组上的位置越唯一。（见下文）
第6列：CIGAR值，read比对的具体情况。（见下文）
第7列：mate比对到的染色体号，若没有mate则是*。
第8列：mate比对到参考序列上的第一个碱基位置，若无mate则为0。
第9列：来自一个fragment的两个reads覆盖的碱基数，若无mate则为0。
第10列：read的碱基序列，如果是比对到互补链上则是reverse completed。
第11列：read质量的ASCII编码。

第2列：flag

0：read是Single-read且正向比对
1：read是pair中的一条
2：pair一正一负完美的比对上
4：这条read没有比对上
8：mate没有比对上
16：这条read反向比对
32：mate反向比对
64：这条read是read1
128：这条read是read2
256：secondary比对
512：比对质量不合格
1024：read是PCR或光学copy产生
2048：辅助比对结果

第6列：CIGAR值

M表示 match或 mismatch
I表示 insert
D表示 deletion
N表示 skipped（跳过这段区域）
S表示 soft clipping（被剪切的序列存在于序列中）
H表示 hard clipping（被剪切的序列不存在于序列中）
P表示 padding
=表示 match
X表示 mismatch（错配，位置是一一对应的）

第5列：
实际上hisat2比对的结果MAPQ没有实际意义，不表示-10log10比对错误率，它只有0、1、60三个值。

• 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）格式呈现。有GTF文件后，就可以利用注释信息计算每个基因/转录本/外显子比对了多少reads，从而获取counts值。GTF是GFF的便于传输版。

GTF格式各列的内容为：
seqname - 染色体或scaffold的名称。
source - 生成这个特征的项目名称，或数据库来源。
feature - 特征类型名称，如gene、transcript、exon、CDS。
start - 开始位置，使用基于1的索引。
end - 结束位置，使用基于1的索引。
score - 组装的转录本的可信度分数（目前该字段未被StringTie使用，如果转录本与read alignment bundle有连接，StringTie将报告一个常量值1000）。
strand - 正链或负链+/-。
frame - 密码子的第几个碱基0/1/2（StringTie不使用这个字段，只记录一个点.）。
attribute - 附加信息。

例如：
a. 来自ensembl的果蝇的Drosophila_melanogaster.BDGP6.32.109.gtf。

feature 的内容是gene、transcript、exon、start_codon、stop_codon、five_prime_utr、three_prime_utr、CDS、Selenocysteine。
attributes的内容例如：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. 来源UCSC的人类GTF文件，1号染色体是 chr1。

c. 来源Ensembl的Homo_sapiens.GRCh38.chr.gtf.gz，1号染色体是1。

d. Stringtie得到的sample.gtf。

GFF格式各列的内容为：
seqid - 染色体或scaffold的名称。
source - 生成这个特征的项目名称，或数据库来源。
feature - 特征类型名称，来自SOFA sequence ontology。
start
end
score
strand - 正链或负链+/-。
phase - 密码子的第几个碱基0/1/2。
attribute - 附加信息。A semicolon-separated list of tag-value pairs。

GTF和GFF之间的区别：
数据结构：都是由9列构成，分别是reference sequence name; annotation source; feature type; start coordinate; end coordinate; score; strand; frame; attributes.前8列都是相同的，第9列不同。
GFF第9列：都是以键值对的形式，键值之间用“=”连接，不同属性之间用“；”分隔，都是以ID这个属性开始。
GTF第9列：同样以键值对的形式，键值之间是以空格区分，值用双引号括起来；不同属性之间用“；”分隔；开头必须是gene_id, transcript_id两个属性。

三、 hisat2使用方法

(1) 用法和参数

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

#提取外显子和剪接位点
extract_splice_sites.py  dmel.annotation.gtf  >  dmel.ss 
extract_exons.py  dmel.annotation.gtf  >  dmel.exon 
extract_snps.py  snp142Common.txt  >  genome.snp

#构建参考基因组索引
hisat2-build  --ss  dmel.ss  --exon  dmel.exon  (--snp  genome.snp)  dmel.genome.fa  dmel  #得到8个索引文件，以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>}：对于单端数据，采用-U指定输入文件（-U <r>）；对于双端数据，采用-1和-2分别指定R1端和R2端的输入文件（-1 <m1> -2 <m2>）。
-x：指定基因组索引的文件前缀（前缀.1~8.ht2）
-S：指定输出的SAM文件<sam>。默认输出到标准输出，比对结束后统计结果输出到标准错误输出。
-p：线程数。
--dta：输出专门为转录本组装的比对结果。Reports alignments tailored for transcript assemblers.（如果比对结果后续需要使用StringTie进行转录本组装，则一定要加入--dta选项。）
--rna-strandness：默认unstranded，如果是链特异性文库要加上--rna-strandness RF。

-f：表示输入文件格式为fasta（默认），-q表示输入文件格式为fastq。可以是gzip压缩的文件。
-phred33：输入的FASTQ文件碱基质量值编码标准为phred33（默认）。
-k <int>：report up to <int> alignments per read; MAPQ not meaningful.

(2) 举个例子

#提取外显子和剪接位点
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

#构建参考基因组索引
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个索引文件

# 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  #双端，链特异性文库

# SAM转为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) 质控和建立索引

#批量对bam文件进行QC
ls 04_bam/*.bam | while read id ; do (samtools flagstat $id > 04_bam/$(basename $id ".bam").stat) ; done  &
#每一个bam文件都生成了一个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)

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

(5) 结果解读

通过2>>log将比对率等信息的标准错误输出重定向到log文件。

查看SAM结果：less 04_bam/WTF4.Input.sam

查看BAM结果：samtools view 04_bam/WTF4.Input.sorted.bam | less

四、 samtools使用方法

(1) 将sam文件转化为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格式，默认输出SAM格式。
-S：出入为sam格式，默认输入文件是bam文件。新版本自动检测，可忽略该参数。
-o：输出文件的名字。
-h：设定输出sam文件时带header信息，默认输出sam不带header信息。
-H：仅输出头文件信息。
-@：线程数。
-F：过滤掉的flag。例如，数字4代表该序列没有比对到参考序列上；数字8代表该序列的mate序列没有比对到参考序列上。-F 4表示比对到参考序列上的结果，-F 12表示双端都比对到参考序列的结果。
-f：需要的flag。-f 4表示没有比对到参考序列上的比对结果。
-q：have mapping quality >= INT
region：提取比对到指定区域上的比对结果

(2) 将bam文件进行排序

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

#版本1.4后可简化直接由sam到sorted.bam
samtools sort  -@ 4  -o sample.sorted.bam  sample.sam

-o：输出文件的名字。
-O：指定输出格式 sam/bam/cram。默认基于-o文件的扩展名选择一个格式，如.bam。如果格式无法推测，则需要指定-O。
-n/-t：默认按照染色体位置（最左边的坐标）进行排序。-n是根据read名（QNAME）进行排序，-t 根据TAG进行排序。
-@：线程数。
-m：内存参数，默认是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文件对应的.bai索引文件。

(4) 过滤

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已经过时了，新版使用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: 除去重复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个或2个以上的已经sort了的bam文件融合成一个bam文件。
-h：指定文件内的头文件复制并替换输出文件的头文件，默认从第一个输入文件复制过来。
-R：合并输出文件的指定区域。

(6) faidx

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

对fasta文件建立索引，生成的索引文件以.fai后缀结尾。由于有索引文件，可以很快从基因组中提取到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文件的比对结果：总reads数，总reads匹配率，多少paired reads，reads1中的reads数，reads2中的reads数，完美匹配的reads数，paired reads中两条都比对到参考序列上的reads数，单独一条匹配到参考序列上的reads数（和上一个相加是总的匹配上的reads数），paired reads中两条分别比对到两条不同的参考序列的reads数。

(8) depth

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

得到每个碱基位点的测序深度，并输出到标准输出。需要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进行SNP和Indel的分析。
mpileup不使用-u或-g参数时，不生成二进制的bcf文件，而生成一个文本文件，输出到标准输出。该文本文件统计了参考序列中每个碱基位点的比对情况，每一行代表了参考序列中某一个碱基位点的比对结果。
mpileup生成的结果包含6行：参考序列名；位置；参考碱基；比对上的reads数；比对情况；比对上的碱基的质量。
-g：计算基因型的似然值和输出文件格式为BCF。
-f：输入有索引文件的fasta参考序列。
-D：输出每个样本的reads深度。