序列搜索算法的进化2

Next generation sequencing era

第二代测序技术（NGS）出现后发生了重大的技术转变。测序仪产生非常短的DNA片段，称为reads（35-400nt），数量巨大（每次运行数十GB的测序数据），成本低。其应用包括：研究全基因组（genome-wide）、RNA转录本（RNA-seq）、DNA甲基化（Methyl-seq）、蛋白质-DNA相互作用（ChIP-seq）和蛋白质组学图谱（proteomic profiles）（外显子测序，exome sequencing）。

处理NGS测序数据最常见的任务是读段映射（read mapping）：在参考基因组序列中定位reads并计算相应的比对（locating reads within a reference genomic sequence and computing corresponding alignments）。数以百万计的reads必须与一个包含数十亿个字母的序列进行比对，类似BLAST的工具无法在合理的时间内处理这样的数据。

常用的read mapping软件包括：BWA-MEM、Bowtie2、GEM，它们比类似BLAST的工具快几个数量级。这些算法的一个主要效率（efficiency）来源是被称为BWT-或FM-index的索引结构（indexing structure）。

BWT-index基于Burrows-Wheeler transform，从组合学（combinatorics）的角度BWT与Gessel-Reutenauer bijection（双射，一一对应关系）有着密切的联系。虽然BWT的主要应用是文本压缩（text compression），但它也适用于构造一个简洁的文本索引（compact text index）。与以BLAST为代表的基于哈希（hash）的索引基础结构支持种子和扩展策略（seed-and-extend strategies）不同，BWT索引是一个全文索引（full-text index），因为它支持搜索任意大小的字符串。BWT-index以位（bits）为单位的大小接近于序列的无损表示所需的信息理论最小值。其他流行的全文索引（如后缀树（suffix tree）、后缀数组（suffix array）及其变体）也被用于生物信息学算法，但需要更多的空间。

BWT算法是一种数据转换算法，它将一个字符串中的相似字符放在相邻的位置，以便于后续的压缩。BWT算法可以分为编码和解码两部分。编码后，原始字符串中的相似字符会处在比较相邻的位置；解码就是将编码后的字符串重新恢复成原始字符串的过程。BWT的一个特点就是经过编码后的字符串可以完全恢复成原始字符串（无损压缩）。

Burrows和Wheeler在1994年提出了字符串匹配和压缩索引的BWT方法，应用于无损压缩。BWT通过对字符串循环移位得到一个字符矩阵，然后通过排序和变换得到一个新的字符串，只改变了字符串中字符的顺序，未改变字符本身。转换后的字符串中许多连续相同的字符会被重组在一起，更容易被压缩。解码时有序地按字典序排列BWT转换后的字符串，即可还原未经转换前的字符串。在BWT结构和技术的基础上，Ferragina和Manzini在2000年提出了一种后缀查找算法FM-index，其搜索查找的时间复杂度为O(m)（m为查找串的长度），并且与参考库中序列长度无关。

编码：输入字符串→加标记$（$比字符串里的所有字符都要小）→循环移位（重复将最后一个字符移到开头）→得到的新字符串从小到大排序→排序后的矩阵中最后一个字符构成L列。

解码：输入L列；对L列进行从小到大排序得到F列；根据L列的字符Li找到F列中的相同字符Fj，然后得到Fj所在行的最后一个字符Lj。将Lj记录下来；重复上面一步，直到Fj等于标记字符为止；按照上述步骤找到的各个Lj进行反向排列，得到字符串r。一般地，如果Li (i=1,2,…,n)这个字符在F列中出现过多次，分别是Fj, Fj+1, …；并且假设L1到Li里这i个字符中一共有k个(k<=i)字符等于Li这个字符，那么我们要为Li在F列中找的对应的字符就是Fj+k-1。

read mapping的结果：得到SAM（sequence alignment/map）格式，或它的二进制（binary）BAM格式。可以上传到UCSC Genome Browser。使用SAMtools和BEDTools可以识别变异（call variants）、堆叠读段（pile reads）、将SAM/BAM转换为BED格式（浏览器扩展数据，Browser Extensible Data）。BEDTools能够发现重叠（finding overlaps）、计算覆盖度（computing coverage）。SAMTools能够对SAM文件进行排序、合并和索引（sorting, merging, and indexing），在每个位置生成比对，能够发现SNPs（single-nucleotide polymorphisms）。

无比对方法（alignment-free methods）

宏基因组学（metagenomics）任务可能是将数百万reads比对到数千计的微生物基因组（microbial genomes），以阐明正在研究的环境样本的组成。即使是专门的比对软件也可能不满足这项任务的实际时间要求。另外，宏基因组学可能不需要计算比对，通常只需要识别给定read可能起源于哪个基因组。一种方法是放弃比对，降低精度，以加快速度和节约内存。

一个常见的想法是将序列看作是出现在其中的一（多）组模式（a (multi-)set of patterns）。在最简单的情况下，这些模式是固定长度k的子字符串（substrings）（k-mers）。然后通过相应集合之间的相似性来表示序列之间的相似性。这种方法被称为无比对（alignment-free）或基于组合的比较（composition-based comparison）。大多数用于全基因组测序数据的现代宏基因组分类器（classifiers）都是基于k-mer分析的。并且可以通过将k-mers替换为spaced k-mers来提高精度（accuracy）。

k-mers是包含在生物序列中的长度为k的子序列。mer即为monomeric unit（单体单元）。长度为L的序列对于一个给定的K可以得到L-k+1个k-mers。每个位点的碱基可以为（A、T、C、G）中的任意一个，因此k-mer理论上有4^k个不同的序列。由k-mer构成的特征向量（x1, x2, …, xn），其中理论的n = 4^k。例如，如果k=4，x1=AAAA，x2=AAAT，...，x256=GGGG。

基因组组装中，二代测序技术reads可能有100多个碱基，如果直接用reads进行组装，原始reads中可能出现的错误就会累加到最终的序列中，导致结果的不准确性。通过将reads切割成以更短的k为单位的k-mer，由于测序错误具有随机性，这些包含测序错误的k-mer绝大多数都是原测序物种中不存在的k-mer，只出现了1次，将这些k-mer去掉可以使组装的结果更加可靠。切割成相同长度的较小k-mer也有助于缓解不同起始read长度的问题。另外，能够利用k-mer对物种基因组特征进行评估。例如：评估物种基因组大小、杂合度和重复序列比例、物种样品污染等。

De Bruijn graph是组装基因组最常用的算法之一（其他还有Overlap-Layout-Consensus，String Graph等）。De Bruijn graph就是通过将reads打断成k-mer后，利用k-mer之间的重复部分构建图，得到最优化路径从而拼接contig。寻找k-mer之间的重叠关系，建立De Bruijn图，即对于任意两个k-mer，如果u的后k-1个碱基序列与w的前k-1个碱基序列相同，则建立一条由u指向w的有向边（edge）。寻找欧拉路径（Eulerian path）来获得结果序列Contigs。实际生成的De Bruijn graph非常复杂，例如由于测序错误、杂合或者高重复序列产生的tips翼尖（a）和bubbles气泡（c），由于低覆盖率造成的链接low coverage links（b）和（d），微小重复造成的压缩（e）。因此需要对DBG图进行简化，移除错误链接、删除低覆盖链接、解开短重复序列等，最终输出Contigs序列。

根据overlap构图并搜寻contigs仅是第一步，完整的基因组组装流程还包含Scaffold搭建、gap修补等过程。双末端测序中，相同序列名字的read1和read2来自同一条序列（中间不一定有overlap），可以根据paired-end信息将不同的Contigs搭建成Scaffold。得到的Scaffold中间会有较多的gap，为了使组装的序列更完整，需再次利用测序的双末端数据之间的配对关系连接Contigs，并利用测序数据与已经组装的Contig之间的覆盖关系对Contig之间空隙进行补洞，延长Contigs。

k-mer大小的选择对基因组组装有多种影响。较小的k-mer将减少图中存储的edges数，节约内存；增加所有k-mer重叠的机会；但是会面临多顶点通向单个k-mer的风险，信息也会丢失；无法解决DNA中出现小微卫星或重复序列问题。较大的k-mer会增加存储DNA序列所需的内存；顶点的数目减少有助于基因组的构建，因为路径变少了；较大的k-mer也会有较高的风险，即没有从每个k-mer出发的向外顶点；导致大量较小的contigs。

一般来说，选用17-mer来估算基因组大小，因为ATCG四种不同的碱基组成长度为17的核苷酸有4^17=17G，足以覆盖一般大小的基因组；如果选择15则只有1G，对正常的基因组来说可能覆盖度不够，导致估计不准。越长的k-mer片段越具很强的物种特异性。当基因组中有较多的重复序列时，可以使用较大的k-mer来跨过高重复的区域，从而获得更加准确及完整的基因组草图。由于reads上的碱基错误率的存在，选择较长的k-mer会带来较高的错误率，但这也可以加大测序深度来弥补。用奇数的k-mer进行分析，是为了防止组装时正反链混淆，奇数的k-mer已经被证明是不能够匹配其反义互补链的。

k-mer常用的工具：JELLYFISH、KAT(The K-mer Analysis Toolkit)、KmerGenie、KMC、Gerbil、Velvet、ABySS、AllPath-LG、SOAPdenovo、EulerUSR、IDBA、SPAdes、Ray。

除了BWT-index外，还有几种数据结构建立在Bloom filters的基础上。Bloom filters已被应用于表示de Bruijn graphs。Bloom filters的扩展包括Cascading Bloom Filter、Bloom Filter Trie、Sequence Bloom tree及其变体。Bloom filters的替代方案有节省空间（space-efficient）和缓存友好（cache-friendly）的无损哈希方案（lossless hashing schemes），例如counting quotient filters。

Bloom过滤器是一种空间高效的（space-efficient）数据结构，具有快速的查找时间和产生假阳性的可控风险。它最初是由Burton Bloom在20世纪70年代开发的，目的是减少容纳哈希编码信息所需的空间。Bloom filter可以用于检索一个元素是否在一个集合中，它的优点是空间效率和查询时间都远远超过一般的算法，缺点是有一定的误识别率和删除困难（有假阳性，无假阴性）。

Bloom Filter 原理：当一个元素被加入集合时，通过K个散列函数将这个元素映射成一个位数组（Bit Array）中的K个点，把它们置为1（它实际上是一个很长的二进制向量和一系列随机映射函数）。检索时，我们只要看看这些点是不是都是1就（大约）知道集合中有没有它了：如果这些点有任何一个0，则被检元素一定不在；如果都是1，则被检元素很可能在。

Hash（散列、哈希）函数，就是把任意长度的输入通过算法变成固定长度的输出。直观解释起来，就是对一串数据m进行杂糅，输出另一段固定长度的数据h，作为这段数据的特征（指纹）。例如，可以简单地为字符串中的每个字符添加代码值，并将结果mod（取余数）返回给定值。散列函数总是从相同的数据生成相同的散列值，但实际上可以为两个不同的数据值生成相同的散列值。也就是说，哈希值对于给定的数据项不是唯一的，并且无法反转哈希函数来获取数据值。Bloom Filter与单哈希函数Bit-Map不同之处在于：Bloom Filter使用了k个哈希函数，每个字符串跟k个bit对应。从而降低了冲突的概率。

所谓的BitMap就是用一个bit位来标记某个元素所对应的value，而key即是该元素，由于BitMap使用了bit位来存储数据，因此可以大大节省存储空间。一个byte是占8个bit，每一个bit的值是有或者没有，也就是二进制的1或者0。