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

发布时间:2021-11-22 19:33:04 作者:柒染
来源:亿速云 阅读:482

这篇文章将为大家详细讲解有关如何解析对比基因组工具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就分享到这里了,希望以上内容可以对大家有一定的帮助,可以学到更多知识。如果觉得文章不错,可以把它分享出去让更多的人看到。

推荐阅读:
  1. nali工具解析ip来源
  2. python3中深浅copy对比解析

免责声明:本站发布的内容(图片、视频和文字)以原创、转载和分享为主,文章观点不代表本网站立场,如果涉及侵权请联系站长邮箱:is@yisu.com进行举报,并提供相关证据,一经查实,将立刻删除涉嫌侵权内容。

上一篇:如何理解IMGT和HLA数据库

下一篇:c语言怎么实现含递归清场版扫雷游戏

相关阅读

您好,登录后才能下订单哦!

密码登录
登录注册
其他方式登录
点击 登录注册 即表示同意《亿速云用户服务条款》