RNAseq测序reads定位
获得RNA-seq的原始数据后,首先需要将所有测序读段通过序列映射(mapping)定位到参考基因组上,这是所有后续处理和分析的基础.在读段定位之前,有时还需要根据测序数据情况对其做某些基本的预处理.
例如,过滤掉测序质量较差的读段,对miRNA测序读段数据去除接头序列等.
高通量测序的海量数据对计算机算法的运行时间提出了很高的要求.针对诸如Illumina/Solexa等测序平台得到的读段一般较短、且插入删除错误较少等特点,人们开发了一些短序列定位算法.这些算法主要采用空位种子索引法(spaced-seedindexing)或Burrows-Wheeler转换(Burrows-WheelerTransform,BWT)技术来实现.空位种子索引法首先将读段切分,并选取其中一段或几段作为种子建立搜索索引,再通过查找索引、延展匹配来实现读段定位,通过轮换种子考虑允许出现错配(mismatch)的各种可能的位置组合.BWT
方法通过B-W转换将基因组序列按一定规则压缩并建立索引,再通过查找和回溯来定位读段,在查找时
可通过碱基替代来实现允许的错配.表1列出了目前可免费下载使用的部分短序列定位软件.其中采用空位种子片段索引法的代表是Maq,而采用Burrows-Wheeler转换的代表是Bowtie.总的来说,采用BWT的定位算法在时间效率上要优于空位种子片段索引法.随着读长的增加,允许读段序列中存在插入删除(indel)的定位变得可行而重要.由于以上两类方法对序列中插入删除的处理较为困难,近来人们开发了一些基于改进的Smith-Waterman动态规划算法的序列比对工具,如BFAST、SHRiMP、Mosaik(http://bioinformatics.bc.edu/marthlab/Mosaik)等,但算法速度较慢,大多需采用计算机并行编程技术来解决运行时间的问题. 下面是 MAQ, Bowtie, BWA, ZOOM, ELAND, SOAP2, RazerS, Novoalign, SHRiMP, BFAST, 以及 Mosaik等mapping软件格式算法的比对信息。
在RNA测序数据的基因组定位中,一个特殊的问题是跨越两个外显子接合区的读段(junctionreads)的定位.在真核生物中,成熟的mRNA是经过由mRNA前体中的外显子经过剪接形成的.如果一个读段跨越了两个外显子,那么就无法将这个读段完整地定位到基因组序列上.而同时,这种跨两个外显子的读段在分析转录本的剪接形式和研究选择性剪接中有重要的作用.为了解决这一问题,人们采取两种典型的策略来进行接合区读段的定位:一是根据已知的基因外显子注释,构建所有可能的外显子接合区序列,与基因组序列一并作为定位的参考基因组;二是不依赖基因注释,而是先利用能完整定位到基因组的读段得到粗略的外显子区域,并结合剪接位点序列构建出可能的剪接位点,然后将不能完整定位的读段分段定位到两个外显子可能的结合区域.Illumina/Solexa平台提供的RNA-seq软件分析包GApipeline采用了第一种策略.采用第二种策略的软件有Tophat和G-Mo.R-Se等,最新的Tophat软件增加了利用已知外显子边界注释信息的选项.
不论是哪种测序平台,测序中都不可避免地存在一定的错误,基因组中又存在单核苷酸多态性等引起的序列变化,所以在读段定位时通常允许一定数量的错配,可以根据不同应用调节允许错配的程度.另一方面,由于基因组中重复序列和高相似度序列的影响,某些读段会出现定位到基因组多个位置的情况.这些因素影响了各个读段到基因组的定位质量,在一些新的读段定位算法中,同时给出每个读段与基因组匹配质量.通常在后续处理前,人们将多定位的读段都过滤掉,也有人尝试用适当的策略把多定位读段“分配”到其中某些位置上.
读段定位到基因组后推荐采用SAM(SequenceAlignment/Map)格式或其二进制版本BAM格式来存储.二进制版本可大大节省存储空间,但不能直接用普通文本编辑工具显示.关于SAM格式的详细介绍,可查阅(http://samtools.sourceforge.net/SAM1.pdf)