欢迎来到《生信修行指南》！

由于测序仪机器读长的限制，在建库过程中需要先对DNA进行片段化，测序得到的序列只是基因组上的部分序列。 为了确定测序读数在基因组上的位置，需要将读数与参考基因组进行比较，这一步称为。

这样做时，需要考虑以下因素

1、硬件资源消耗

一般来说，基因组越大，占用的显存就越大。 优化内存消耗对于小型基因组（例如人类基因组）至关重要。

2、运行速度

随着测序价格的提高和数据深度挖掘的需要，测序量越来越大，海量测序reads的比对必须足够快。

3.确定性

SNP/indel、测序错误率等诱因会导致测序的reads与基因组上的原始序列存在数个bp的偏差，因此算法必须支持核苷酸错配，即gap的存在。 同时，由于测序的短序列可能与基因组中的多个位置存在同源性，因此一个read会与基因组中的多个位置进行比对。 双端测序技术在一定程度上可以校准多个位置。 由于-end reads来自同一个DNA片段，两者在基因组上的位置不会相差太远，仅靠这一点并不能解决所有的同源比。 是的，这就需要比对算法对多个位置进行判别和打分，才能给出比对结果的可靠性。

4.核糖核酸

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

目前的工具有很多，如bwa、hisat、star等，其中hisat速度最快，是软件的升级版。 使用改进后的算法，对于人类基因组来说，只需要大约4.3GB的显存。同时支持DNA和RNA数据的比对。 软件官网如下

目前最新版本。 安装过程如下

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文件获取cut site和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

支持多种格式的输入文件，常用格式如下

法斯塔

快速地

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

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

读取映射到基因组上的一个位置，我们称之为 . 软件会对所有的项目进行评分和判断，只有满足过滤条件的才称为有效，只有有效才能输出。

与blast类似，每个玩家也有相应的计分机制。 hisat评分从以下几个方面

1.错配核苷酸惩罚

不匹配核苷酸的惩罚由--mp参数指定，它的值是由冒号分隔的两个数字，第一个数字是最大惩罚，第二个数字是最小惩罚

2. reads的gap

空位罚分为两部分，第一空位罚分和空位延长罚分。 读取的间隙惩罚由 --rdg 参数指定，其值是由冒号分隔的两个数字。 一个数字是空位第一位置的罚分，第二个数字是空位延伸的罚分。

3.差距处罚

上面的gap 由--rdg参数指定，它的值是两个用冒号隔开的数字。 第一个数字是空位第一位置的罚分dnastar序列比对，第二个数字是空位延伸的罚分。

经过一系列惩罚机制后，每个人都会有一个对应的分数，然后通过一个阈值来判断这个分数是否满足有效要求。

hisat 通过--score--min 参数指定阈值。 指定方法是与读取程度相关的函数。 默认值为L、0、-0.2dnastar序列比对，对应函数为

f(x) = 0 - 0.2 * x

根据reads的宽度，可以估计score的阈值，小于阈值的就认为是有效的，可以输出。 L代表线性函数。 此外，还支持其他类型的函数，例如常数、自然对数等，更多选项请参考官方文档。

一次读取可能有多个有效值。 输出时，不会全部输出，只输出-k参数指定的N个。 -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

·结尾·

—喜欢就分享给您的同事吧—

扫一扫关注微信，更多精彩内容等你来！

如有侵权请联系删除！