如何解析对比基因组工具hisat2

这篇文章将为大家详细讲解有关如何解析对比基因组工具hisat2，文章内容质量较高，因此小编分享给大家做个参考，希望大家阅读完这篇文章后对相关知识有一定的了解。

由于测序仪机器读长的限制，在构建文库的过程中首先需要将DNA片段化，测序得到的序列只是基因组上的部分序列。为了确定测序reads在基因组上的位置，需要将reads比对回参考基因组上，这个步骤叫做mapping。

在进行mapping时，需要考虑以下几个因素

1. 硬件资源的消耗

通常来说，基因组越大，占用的内存越大。对于大型基因组，比如人类基因组而言，优化内存消耗是很关键的一点。

2. 运行速度

随着测序价格的下降和数据深入挖掘的需求，测序量越来越大，海量测序reads的比对，要求速度上必须够快。

3. 准确性

SNP/indel, 测序错误率等因素都使得测序的reads和基因组上的原始序列会存在几个bp的误差，所以mapping的算法必须支持碱基的错配，或者是gap的存在。同时由于测序的短序列可能和基因组多个位置存在同源，一条reads会比对到基因组上多个位置。双端测序技术在一定程度上能够校正多个位置，因为双端reads 来自同一个DNA片段，二者在基因组上的位置不会相距太远，但是仅靠这一点并不能解决所有的同源比对，这就要求比对算法对多个位置进行判断和打分，给出比对结果的可靠性。

4. RNA

对于转录组数据， 真核生物可变剪切的存在，导致cDNA片段在基因组上的位置并不是连续的，中间可能存在内含子。在比对转录组数据时，就需要考虑跳过剪切位点。

目前mapping的工具有很多，比如bwa, hisat, star等。hisat 是其中速度最快的，是tophat软件的升级版本。采用了改进的FM index 算法，对于人类基因组，只需要4.3GB左右的内存。同时支持DNA和RNA数据的比对，软件官网如下

http://ccb.jhu.edu/software/hisat2/index.shtml

目前最新版为为hisat2. 安装过程如下

wget ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/downloads/hisat2-2.1.0-Linux_x86_64.zip
unzip hisat2-2.1.0-Linux_x86_64.zip

下载解压缩即可。

在进行比对前，首先需要对参考基因组建立索引, 基本用法如下

hisat2-build -p 20   hg19.fa hg19

对于转录组数据，在构建索引时，可以通过gtf文件，得到剪切位点和exon的信息，用法如下

hisat2_extract_splice_sites.py hg19.gtf > hg19.ss
hisat2_extract_exons.py hg19.gtf > hg19.exon
hisat2-build -p 20  --ss hg19.ss --exon hg19.exon  hg19.fa hg19

hisat2 支持多种格式的输入文件，常见格式有以下两种

1. fasta

2. fastq

   
   

-f参数表示输入问下格式为fasta, -q参数表示输入文件格式为fastq。输入文件可以是经过gzip压缩之后的文件，默认输入文件是fastq格式。

对于单端数据，采用-U指定输入文件；对于双端数据，采用-1和-2分别指定R1端和R2端的输入文件。

reads比对到基因组上的一个位置，我们称之为一个alignment。 软件会对所有的alignments 进行打分和判断，能够符合过滤条件的alignment 称之为valid  alignment， 只有valid alignments , 才会输出。

和blast类似，每个alignment也有对应的打分机制。hisat 从以下几个方面对alignment 进行打分

1. 错配碱基罚分

错配碱基的罚分通过--mp参数指定，其值为逗号分隔的两个数字，第一个数字为最大的罚分，第二个数字为最小的罚分

2. reads上的gap罚分

gap的罚分通过分成两个部分，第一次出现gap的罚分和gap延伸的罚分，reads上的gap罚分通过--rdg参数指定，其值为逗号分隔的两个数字，第一个数字为gap第一个位置的罚分，第二个数字为gap延伸的罚分。

3. reference上的gap罚分

reference上的gap罚分通过--rdg参数指定，其值为逗号分隔的两个数字，第一个数字为gap第一个位置的罚分，第二个数字为gap延伸的罚分。

经过一系列的罚分机制，每个alignment会有一个对应的得分，然后会根据一个阈值，来判断这个得分是否满足valid  alignment的要求。

hisat通过--score--min参数指定该阈值，指定方式是一个和reads程度相关的函数，默认值为L,0,-0.2, 对应函数为

f(x) = 0 - 0.2 * x

根据reads长度，可以计算出得分的阈值，大于该阈值的alignment 被认为是valid alignment , 才可能被输出。L代表线性函数，此外，也支持其他类型的函数，比如常量，自然对数等，更多选择请参考官方文档。

一条reads可能会拥有多个valid  alignments, 在输出时，并不会输出所有的alignments, 而是只输出-k参数指定的N个alignments，-k参数的默认值为5。
输出结果以SAM格式保存，默认输出到屏幕上，可以通过-S参数指定输出文件。

通常情况下，默认参数就能够满足我们的需求了。单端数据比对的用法如下

hisat -x hg19 -p 20 -U reads.fq -S align.sam

双端数据用法如下

hisat -x hg19 -p 20 -1 R1.fq -2 R2.fq -S align.sam

关于如何解析对比基因组工具hisat2就分享到这里了，希望以上内容可以对大家有一定的帮助，可以学到更多知识。如果觉得文章不错，可以把它分享出去让更多的人看到。