转录组分析 | 使用Hisat2进行序列比对
我前面已经对数据进行了质控：
转录组分析 | 使用trim-galore去除低质量的reads和adaptor
接下来我们进行序列比对，利用的软件是Hisat2。
一.index文件下载
index文件直接去官网下载
(http://daehwankimlab.github.io/hisat2/download/)，我测序的组织来自小鼠，所以我这里下载的是小鼠的。
下载完后上传到Linux服务器，解压后就可以直接用了。
我上传到了： /data/mouse_genome/ ，就是mm10_genome.tar.gz这个文件。
解压文件，解压过程中会在当前文件夹下创建mm10文件，解压后的文件就在mm10文件夹中。
tar -zxvf /data/mouse_genome/mm10_genome_tar.gz查看解压后的文件。
ll -h ./mm10总共8个ht2格式文件，一个sh格式文件。
二.hisat2介绍
Hisat是一种高效的RNA-seq实验比对工具。它使用了基于BWT和Ferragina-manzini (Fm) index 两种算法的索引框架。使用了两类索引去比对，一类是全基因组范围的FM索引来锚定每一个比对，另一类是大量的局部索引对这些比对做快速的扩展。比对原理可阅读文献原文：HISAT: a fast spliced aligner with low memory requirements.
（一）背景
自2008年起，RNAseq已经成为研究基因表达、转录本结构、长链非编码RNA确定以及融合转录本的重要手段。随着测序深度的加深和read读长的延长，给比对工作带来很多困难。目前的工具如Tophat2和GSNAP等在对单个转录组实验的比对中耗时几天，而新型的STAR虽然很快，但是会吃掉大量的内存空间，如基于人类基因组的比对需要消耗28G左右，而hisat2是4.3G，小鼠3.8G左右。
（二）Hisat设计原则
1.优化了索引建立策略
hisat应用了基于bowtie2的方法去处理很多低水平的用于构建和查询FM索引的操作。但是与其它比对器不同的是，该软件应用了两类不同的索引类型：代表全基因组的全局FM索引和大量的局部小索引，每个索引代表64000bp。以人类基因组为例，创建了48000个局部索引，每一个覆盖1024bp，最终可以覆盖这个3 billion 的碱基的基因组。这种存在交叉（overlap）的边界可以轻松的比对那些跨区域的read（可变剪切体）。尽管有很多索引，但是hisat会把他们使用合适方法压缩，最终只占4gb左右的内存。
2.采用了新的比对策略
RNA-seq产生的reads可能跨长度比较大的内含子，哺乳动物中甚至最长能达到1MB，同时外显子比较短，read也比较短，会有很多read（模拟数据中大概34%）跨两个外显子的情况。为了更好的比对，将跨外显子的reads分成了三类：1）长锚定read，至少有16bp在两个外显子的每一个上 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%，而且比对错误率很高！【引用：1】
官方手册：https://daehwankimlab.github.io/hisat2/manual/
（三）.HiSat2进行比对的参数设置 【引用：2】
HISAT2 version 2.1.0 by Daehwan Kim (infphilo@gmail.com, www.ccb.jhu.edu/people/infphilo)
Usage:
hisat2 [options]* -x <ht2-idx> {-1 <m1> -2 <m2> | -U <r>} [-S <sam>]
<ht2-idx> Index filename prefix (minus trailing .X.ht2).
索引文件的前缀
<m1> Files with #1 mates, paired with files in <m2>.
Could be gzip'ed (extension: .gz) or bzip2'ed (extension: .bz2).
Paired-end测序的第一个文件，可以是压缩格式的（`.gz`或者`.bz2`）
<m2> Files with #2 mates, paired with files in <m1>.
Could be gzip'ed (extension: .gz) or bzip2'ed (extension: .bz2).
Paired-end测序的第二个文件
<r> Files with unpaired reads.
Could be gzip'ed (extension: .gz) or bzip2'ed (extension: .bz2).
single-end测序的文件
<sam> File for SAM output (default: stdout)
输出比对的结果（`.sam`）
<m1>, <m2>, <r> can be comma-separated lists (no whitespace) and can be
specified many times. E.g. '-U file1.fq,file2.fq -U file3.fq'.
可以是多个序列文件用逗号分隔
Options (defaults in parentheses):
Input:
-q query input files are FASTQ .fq/.fastq (default)
输入检索序列为FASTQG格式
--qseq query input files are in Illumina's qseq format
输入检索序列是Illumina qseq格式
-f query input files are (multi-)FASTA .fa/.mfa
输入文件是多个FASTA文件
-r query input files are raw one-sequence-per-line
输入检索序列是每行一个原始序列
-c <m1>, <m2>, <r> are sequences themselves, not files
直接输入序列，而不是文件
-s/--skip <int> skip the first <int> reads/pairs in the input (none)
忽略输入序列的前`<int>`个
-u/--upto <int> stop after first <int> reads/pairs (no limit)
只分析输入序列的前`<int>`个
-5/--trim5 <int> trim <int> bases from 5'/left end of reads (0)
剪去reads长度为`<int>`的5'端序列
-3/--trim3 <int> trim <int> bases from 3'/right end of reads (0)
剪去reads长度为`<int>`的3'端序列
--phred33 qualities are Phred+33 (default)
序列质量打分为Phred+33
--phred64 qualities are Phred+64
序列质量打分为Phred+64
--int-quals qualities encoded as space-delimited integers
序列质量打分为用空格分隔的整数
Alignment:
--n-ceil <func> func for max # non-A/C/G/Ts permitted in aln (L,0,0.15)
允许在比对结果中出现的最多非ACTG的计算函数
--ignore-quals treat all quality values as 30 on Phred scale (off)
将所有质量分值设为30
--nofw do not align forward (original) version of read (off)
不比对前向的read
--norc do not align reverse-complement version of read (off)
不比对后向的read
Spliced Alignment:
--pen-cansplice <int> penalty for a canonical splice site (0)
对已知剪接位点的罚分
--pen-noncansplice <int> penalty for a non-canonical splice site (12)
对新的剪接位点的罚分
--pen-canintronlen <func> penalty for long introns (G,-8,1) with canonical splice sites
对已知剪接位点中长内含子的罚分
--pen-noncanintronlen <func> penalty for long introns (G,-8,1) with noncanonical splice sites
对未知剪接位点中长内含子的罚分
--min-intronlen <int> minimum intron length (20)
最小允许的内含子长度
--max-intronlen <int> maximum intron length (500000)
最大允许的内含子长度
--known-splicesite-infile <path> provide a list of known splice sites
提供一个已知的剪接位点列表
--novel-splicesite-outfile <path> report a list of splice sites
报告新的剪接位点
--novel-splicesite-infile <path> provide a list of novel splice sites
提供新的剪接位点
--no-temp-splicesite disable the use of splice sites found
不使用发现的新剪接位点
--no-spliced-alignment disable spliced alignment
不允许剪接比对
--rna-strandness <string> specify strand-specific information (unstranded)
指定链特异性的信息
--tmo reports only those alignments within known transcriptome
报告只包含已知转录组信息的比对结果
--dta reports alignments tailored for transcript assemblers
结果适合用于stringtie进行转录本的拼接
--dta-cufflinks reports alignments tailored specifically for cufflinks
结果适合于cufflinks的应用
--avoid-pseudogene tries to avoid aligning reads to pseudogenes (experimental option)
避免与假基因的比对
--no-templatelen-adjustment disables template length adjustment for RNA-seq reads
不允许对reads的模板长度进行调整
Scoring:
--mp <int>,<int> max and min penalties for mismatch; lower qual = lower penalty <6,2>
不匹配的最大与最小罚分值，低质量对应低罚分
--sp <int>,<int> max and min penalties for soft-clipping; lower qual = lower penalty <2,1>
soft-clipping的最大与最小罚分
--no-softclip no soft-clipping
不考虑Soft-clipping
--np <int> penalty for non-A/C/G/Ts in read/ref (1)
对read或参考序列中出现的非ATCG的罚分
--rdg <int>,<int> read gap open, extend penalties (5,3)
对read的gap-open, gap-extension的罚分
--rfg <int>,<int> reference gap open, extend penalties (5,3)
对参考序列的gap-open,gap-extension的罚分
--score-min <func> min acceptable alignment score w/r/t read length
(L,0.0,-0.2)
相对read长度的最小接受比对分值
Reporting:
-k <int> (default: 5) report up to <int> alns per read
对每个read只报告`<int>`个比对结果
Paired-end:
-I/--minins <int> minimum fragment length (0), only valid with --no-spliced-alignment
paired-end的最小片段长度
-X/--maxins <int> maximum fragment length (500), only valid with --no-spliced-alignment
最大片段长度
--fr/--rf/--ff -1, -2 mates align fw/rev, rev/fw, fw/fw (--fr)
`--fr==fw/rev`，`--rf==rev/fw`, `--ff==fw/fw`
--no-mixed suppress unpaired alignments for paired reads
不输出paired reads的不能匹配的比对结果
--no-discordant suppress discordant alignments for paired reads
不输出paired-reads的不统一的比对结果
Output:
-t/--time print wall-clock time taken by search phases
输出搜索阶段所花费的wall-clock时间
--un <path> write unpaired reads that didn't align to <path>
输出没有与`<path>`比对的unpaired reads
--al <path> write unpaired reads that aligned at least once to <path>
输出至少与`<path>`有一次比对的unpaired reads
--un-conc <path> write pairs that didn't align concordantly to <path>
--al-conc <path> write pairs that aligned concordantly at least once to <path>
(Note: for --un, --al, --un-conc, or --al-conc, add '-gz' to the option name, e.g.
--un-gz <path>, to gzip compress output, or add '-bz2' to bzip2 compress output.)
将输出结果进行`gz`或者`bz2`压缩
--summary-file print alignment summary to this file.
将比对的总结输出到文件中
--new-summary print alignment summary in a new style, which is more machine-friendly.
将比对总结以新的形式输出
--quiet print nothing to stderr except serious errors
除了严重错误以外，不将任何结果输出到stderr
--met-file <path> send metrics to file at <path> (off)
将度量值i输出到文件
--met-stderr send metrics to stderr (off)
将度量值输出到stderr
--met <int> report internal counters & metrics every <int> secs (1)
每隔`<int>`秒输出当前的内部计数器和度量值
--no-head supppress header lines, i.e. lines starting with @
不输出header行，也就是所有以`@`开始的行
--no-sq supppress @SQ header lines
不输出`@SQ`行
--rg-id <text> set read group id, reflected in @RG line and RG:Z: opt field
设置`@RG`
--rg <text> add <text> ("lab:value") to @RG line of SAM header.
Note: @RG line only printed when --rg-id is set.
将`<text>`加到SAM文件的`@RG`行
--omit-sec-seq put '*' in SEQ and QUAL fields for secondary alignments.
将"*"置于SEQ和QUAL两列中，进行二次比对
Performance:
-o/--offrate <int> override offrate of index; must be >= index's offrate
重置索引的`offrate`，该值必须大于原有的`offrate`
-p/--threads <int> number of alignment threads to launch (1)
比对所用的线程数
--reorder force SAM output order to match order of input reads
强制要求SAM输出顺序比对与输入reads一致
--mm use memory-mapped I/O for index; many 'hisat2's can share
使用内存映射的I/O作为索引
Other:
--qc-filter filter out reads that are bad according to QSEQ filter
根据QSEQ filter过滤输入的reads
--seed <int> seed for random number generator (0)
设置RNG seeds的值
--non-deterministic seed rand. gen. arbitrarily instead of using read attributes
不根据reads的属性产生随机数
--remove-chrname remove 'chr' from reference names in alignment
比对结果中的参考序列名称中去除`chr`
--add-chrname add 'chr' to reference names in alignment
在比对结果的参考序列名称中加入`chr`
--version print version information and quit
输出版本信息
-h/--help print this usage message
输出帮助信息hisat2可以自己建立index文件，如果没有现成的index的话，那就得自己去建立了，这个就比较麻烦而且耗内存和时间，下面是建立index的一些参数介绍。因为我们的index文件是直接下载的，所以我后面不需要建立index文件，关于这部分，后面有机会介绍。
HISAT2 version 2.1.0 by Daehwan Kim (infphilo@gmail.com, http://www.ccb.jhu.edu/people/infphilo)
Usage: hisat2-build [options]* <reference_in> <ht2_index_base>
reference_in comma-separated list of files with ref sequences
以逗号分隔的参考序列文件的列表
hisat2_index_base write ht2 data to files with this dir/basename
输出的索引文件的`basename`
Options:
-c reference sequences given on cmd line (as
<reference_in>)
输入的参考序列
--large-index force generated index to be 'large', even if ref
has fewer than 4 billion nucleotides
强制要求产生的索引为`large`
-a/--noauto disable automatic -p/--bmax/--dcv memory-fitting
关闭-p/--bmax/--dcv等自动内存匹配
-p number of threads
线程数
--bmax <int> max bucket sz for blockwise suffix-array builder
逐块后缀数组建立采用的最大bucket size
--bmaxdivn <int> max bucket sz as divisor of ref len (default: 4)
--dcv <int> diff-cover period for blockwise (default: 1024)
--nodc disable diff-cover (algorithm becomes quadratic)
-r/--noref don't build .3/.4.ht2 (packed reference) portion
-3/--justref just build .3/.4.ht2 (packed reference) portion
-o/--offrate <int> SA is sampled every 2^offRate BWT chars (default: 5)
-t/--ftabchars <int> # of chars consumed in initial lookup (default: 10)
--localoffrate <int> SA (local) is sampled every 2^offRate BWT chars (default: 3)
--localftabchars <int> # of chars consumed in initial lookup in a local index (default: 6)
--snp <path> SNP file name
--haplotype <path> haplotype file name
--ss <path> Splice site file name
--exon <path> Exon file name
--seed <int> seed for random number generator
-q/--quiet verbose output (for debugging)
-h/--help print detailed description of tool and its options
--usage print this usage message
--version print version information and quit三. 使用Hisat2进行序列比对
创建输出数据的文件夹
mkdir cleandata/hisat2_mm10data因为比对这一过程很耗内存，所以样本多话，计算机内存不够大，需要分批比对，我就先介绍一个样本比对的命令：
hisat2 --dta -t -x /data/RNAseq/mm10/genome -1 /data/RNAseq/cleandata/trim_galoredata/CK-4_1_val_1.fq.gz -2 /data/RNAseq/cleandata/trim_galoredata/CK-4_2_val_2.fq.gz -S cleandata/hisat2_mm10data/CK4.sam -1和-2分别表示双端测序的1个文件，后面跟的是文件路径，一定要注意 /data/RNAseq/mm10/genome文件的目录，genome这个不是文件夹，是index文件的前缀，我的mm10文件下并没有这个文件，如果不加genome就会发生如下报错：
## Could not locate a HISAT2 index corresponding to basename "/data/RNAseq/mm10/genome/" ##另外一定要加上 --dta，后续用Stringtie处理会更容易一些，这是StringTie使用说明里面的一句话：【NOTE: be sure to run HISAT2 with the --dta
option for alignment, or your results will suffer.】后续再介绍。-S后面cleandata/hisat2_mm10data/CK4.sam 是输出文件路径，hisat2_mm10data这个文件夹如果没有建立，需要事先建立。
命令没有问题的话，出现下面提示表示程序正在执行，就去干其他的，等执行结束再往下：
Time loading forward index: 00:02:51
Time loading reference: 00:00:16我这里，一个样本43分钟。
计算机内存足够大的话，我们可以像前文【转录组分析 | 使用trim-galore去除低质量的reads和adaptor】一样，通过一个脚本一次执行。
创建脚本文件：hisat2_mm10_batch.sh
vim hisat2_mm10_batch.sh输入下面内容：
#!/bin/bash
# This is for hisat2 batch aligne
for i in CK-7 CK-8 HGJ-10 HGJ-6 HGJ-9
do
hisat2 --dta -t -p 8 -x /data/RNAseq/mm10/genome \
-1 /data/RNAseq/cleandata/trim_galoredata/"$i"_1_val_1.fq.gz \
-2 /data/RNAseq/cleandata/trim_galoredata/"$i"_2_val_2.fq.gz \
-S /data/RNAseq/cleandata/hisat2_mm10data/"$i".sam; done我这里由于CK-4这个样本前面单独处理了，批量处理就不需要处理。
保存脚本，然后执行脚本。
bash hisat2_mm10_batch.sh脚本运行时间和计算机配制还有数据量的大小有关。运行结束后我们查看结果：
ll -h ./cleandata/hisat2_mm10data/原本22G的数据，运行结束后155G，所以自己用本地计算机跑的话，数据多注意一下内存。
四.结果解读
我们比对后得到的是sam文件，关于sam文件是一个什么样的文件，参考文章：生信中常见的数据文件格式。
为了快速查看本公众号文章，可阅读文章：公众号文章目录
致谢：本文还参考了以下内容：
【1】https://www.jianshu.com/p/ce3f4afb9b60。
【2】https://github.com/ricket-sjtu/bioinformatics/wiki/HiSat2%E5%BA%94%E7%94%A8%E4%BB%8B%E7%BB%8D
对原作者表示感谢！
腾讯云开发者
