声明:本文部分内容和部分图片来源于网络。本文为生信小白学习笔记,不能保证专业名词和内容全部正确或权威。
下图为某一条RNAseq从数据预处理,序列回帖到数据可视化的工作流程,包含了较多的软件(Linux环境运行)和若干个包(R语言环境运行),本系列将按下图,对每一个步骤进行学习和理解。
HISAT2
简介
HISAT2是将下一代测序读段结果基于图比对到一组基因组(graph-based alignment of next generation sequencing reads to a population of genomes)。
HISAT2是一种快速而灵敏的比对程序,可用于将下一代测序数据(包括DNA和RNA)比对到人类基因组和单个参考基因组上。基于图的BWT扩展,创造性地设计并完成了一个图FM索引(GFM)。除了使用一个代表全人类基因组的全球GFM索引,HISAT2使用大量小的GFM索引,这些索引共同覆盖了全基因组。这些小的索引(也被称为局部索引),与集中比对方式结合在一起,能够实现快速和准确的序列比对。这个新的索引方案被称为层次图片FM索引(HGFM)。
HISAT2工作原理
1. HISAT2应用了基于bowtie2的方法处理很多低水平的用于构建和查询FM索引的操作。(*)
2. 与其他比对器相比,HISAT2应用了两类不同的索引类型,代表全基因组的全局FM索引和大量的局部小索引,每个索引代表64000bp。
3. 以人类基因组为例,创建了48000个局部索引,每一个覆盖1024bp,最终可以覆盖这个3 billion碱基的基因组。这种存在交叉(overlap)的边界可以轻松的比对那些跨区域的read(可变剪切体)。
4. 尽管有很多索引,但是HISAT2可以把他们使用合适的方式进行压缩,最终只占4GB左右的内存。
模式
报告模式
报告模式管理HISAT2寻找多少个比对以及如何报告它们。
通常,当我们说一个读段有一个比对,是指它有一个有效比对。当我们说一个读段有多个比对时,是指它有多个有效且彼此不同的比对方式。
默认情况下,HISAT2会对5‘和3’端进行温和地剪切。
比对总结
当HISAT2完成运行,会输出运行结果。这些信息将输入到‘标准错误’(stderr)文件中。对包含未匹配读段地数据文件,HISAT2总结可能如下所示:
针对包含已匹配读段的数据文件,HISAT2总结如下所示:
Alignment rate越高表示HISAT2对该文件比对成功率越高。
索引大小
hisat2-build能够索引任何尺寸的参考基因组。对小于40亿个核苷酸长度的基因组,hisat2-build使用32位数字在索引的不同位置建立一个‘小’索引。当基因组更长,hisat2-build能够使用64位数字建立较大的索引。小索引保存在.ht2文件中,而大索引会保存在.ht21文件中。使用者无需担心特定的索引的尺寸,HISAT2中的包装脚本将自动生成并使用合适的索引。
性能调试
如果运行的电脑有多线程或多核,可以使用 -p
-p选项可以使HISAT2启动一定数量的并行搜索线程。每一个线程运行在一个不同的中央处理器或核中,而所有的线程并行地查找比对,将比对量提高了大概并行线程的倍数(虽然在现实中,加速有时比线性较差)。
HISAT2使用
主要参数
⚪ hisat2 [options]* -x <hisat2-idx> {-1 <m1> -2 <m2> | -U <r> | --sra-acc <SRA accession number>} [-S <hit>]
1. -x <hisat2-idx>
参考基因组索引的名字。该名称是任何索引文件的名称。HISAT2会首先寻找在现有文件中特定的索引,然后再在HISAT2_INDEXES指定地环境变量的目录中搜索。
2. -1 <m1>
逗号分隔的文件列表包括了双端测序的文件1,例如,-1 flyA_1.fq,flyB_1.fq。使用此命令指定的文件-文件的顺序必须与<m2>读取-读取的顺序相一致。
3. -2 <m2>
逗号分隔的文件列表包括了双端测序的文件2,例如,-2 flyA_2.fq,flyB_2.fq。对文件顺序的要求同上。
4. -U <r>
逗号分隔的文件列表包含待比对的未成功匹配(unpaired)读段,例如,lane1.fq,lane2.fq,lane3.fq,lane4.fq
5. --sra-acc <SRA accession number>
逗号分隔的SRA登录号文件列表,例如,--sra-acc SRR353653,SRR353654
6. -s <hit>
写入SAM比对结果的文件。
选项
输入选项
比对选项
计分选项
拼接对齐选项
报告选项
双端测序选项
输出选项
SAM选项
性能选项
其他选项
具体选项见链接。
HISAT2比对操作
HISAT2提供了一些示例文件,这些示例文件的结果并不具有科学意义,这些文件只供运行HISAT2和相应的下游分析。
首先是获取和安装HISAT2,并设置相应的环境变量到包含hisat2, hisat2-build和hisat2-inspect的HISAT2目录中。
比对实例读段
从HISAT2网站获取待分析物种参考基因组,下一步将待分析读段比对到参考基因组上。命令如下:
$HISAT2_HOME/hisat2 -f -x $HISAT2_HOME/example/index/22_20-21M_snp -U $HISAT2_HOME/example/reads/reads_1.fa -S eg1.sam
本例使用的是使用hisat2-build构建的索引文件(22_20-21M_snp)。这行命令将一组未配对的读段数据比对到索引上。比对结果被写入进eg1.sam文件中,同时,一段简短的比对总结被写入进console。
可使用下列语句查看SAM文件的前几行。
head eg1.sam
可能会得到下图类似的结果。
上图前几行(以@开始)是SAM文件表头行,其他行是SAM比对结果,每读段或每对读段一行。
双端测序比对
为了使用HISAT2比对双端测序数据,首先,需要需要进入相同更多目录然后运行以下命令:
$HISAT2_HOME/hisat2 -f -x $HISAT2_HOME/example/index/22_20-21M_snp -1 $HISAT2_HOME/example/reads/reads_1.fa -2 $HISAT2_HOME/example/reads/reads_2.fa -S eg2.sam
SAMtools转换文件格式
SAMtools是管理和分析SAM和BAM比对文件的一组工具,提供了一个可以方便转换SAM和BAM文件格式。在HISAT2软件进行序列比对后,可用SAMtools将SAM文件转换为BAM文件,命令如下:
samtools view -bS eg2.sam > eg2.bam
同时,SAMtools也可以转换为BAM文件的同时进行排序(版本需要1.2或更高)。命令如下:
samtools sort eg2.bam -o eg2.sorted.bam
对BAM进行排序时非常有用的,因为比对通常是压缩的,这对于长期存储是很方便的,同时,排序的BAM文件也有助于突变的发现。