计算生物学与生物信息学漫谈-5-mapping算法

之前的文章我们介绍了参考基因组，也介绍了一些基本概念，具体可以看之前的博客：

计算生物学与生物信息学漫谈-4-参考基因组与Mapping准备_基因组的map-CSDN博客

这次我们讲如何将read map到基因组上所用到的各种算法：

1.1 简述

1.2 字典树（Trie）

1.3 后缀树

1.4 后缀数组（SA）

1.5 Burrows-Wheeler变换

1.6 FM-索引

1.1 简述

序列比对是生物信息学中最重要的任务之一。在最近的基因组学革命之前，它已经被使用了数十年，这场革命是在高通量测序技术的发现之后爆发的。它已被用于众多的生物信息学应用中。基本的比对算法包括使用Needleman-Wunsch算法的全局序列比对和使用Smith-Waterman（SW）算法的局部序列比对。这两种传统的算法都使用了动态规划，这是数学家理查德·贝尔曼在1950年代初期发明的。

序列比对中的动态规划与计算机编程无关，它是一种尝试通过将整个比对划分为更小的子比对来找到两个序列最优比对的算法。

全局序列比对是从端到端比对两个序列，因此，可以引入空位以使序列长度相等。很明显，这种类型的比对不适用于将大量短序列比对到可能跨越数百万个碱基的基因组序列。第二种类型是局部序列比对，它仅将两个序列比对到相似的区域。很明显，这种序列比对适用于将读段比对到参考基因组。局部序列比对的原始形式是穷尽的，也就是说，如果你有一个单一的序列，并且想要将其比对到多个序列，算法将穷尽所有可能的比对，然后报告结果。在大多数情况下，这种方法对于比对数百万个读段是不切实际的，而且计算成本高昂，尤其是面对大量的基因组序列时。

基本的局部比对搜索工具或BLAST是SW局部序列比对算法的启发式版本，它比穷尽方法更节省时间，因为它只搜索序列中最显著的模式并据此进行序列比对。

BLAST速度快，但由于其启发式的性质，可能会跳过一些比对。现在，我们可以问一个问题：BLAST是否适合将数百万个短读段序列比对到可能跨越数百万个碱基的参考基因组？为了实际回答这个问题，假设BLAST可以在一秒内将一个读段比对到参考基因组；如果我们有2000万个读段，BLAST将需要2000万秒（约231.5天）来将它们比对到参考基因组，忽略可能的断电、网络中断或其他自然灾害的可能性。总之，BLAST不适合读段序列比对。

考虑到参考基因组序列的大小和读段数量之大，执行搜索和精确定位映射区域的最大挑战之一就是索引。需要高效的索引算法来有效地组织参考基因组的序列和短读段在内存中的存储，并促进快速搜索模式。计算机科学家在开发允许计算机程序有效存储和处理数据的数据结构方面发挥着关键作用，同时减少计算成本。高效的数据结构将有助于对参考基因组的序列进行索引，并在索引序列中快速搜索读段，同时高效地使用内存。

高效的数据结构将有助于对参考基因组的序列进行索引，并在索引序列中快速搜索读段，同时高效地使用内存。已经提出了几种用于索引的数据结构，并在读段序列比对的软件包中实现了它们。

最常用的索引数据结构包括后缀树、后缀数组、Burrows-Wheeler变换（BWT）和全文本最小空间（FM-index）。比对器实现各种搜索算法以确定短读段在索引参考基因组序列中的来源位置。读段与基因组位置的比对取决于序列的相似性。然而，一个好的比对器应该预期到由于测序错误和参考基因组与被测个体之间的遗传变异而导致的不匹配。接下来，我们将讨论这些常用的索引数据结构和流行的搜索方法。

1.2 字典树（Trie）

字典树或前缀树是一种用于在大型数据集（如查找测序读段）上快速检索的数据结构。

这个名字是由E. Fredkin于1960年从“信息检索”这一短语中衍生出来的；然而，这个想法最初是由挪威数学家Axel Thue在1912年描述的。字典树是一个有序树，可以用来表示一组字符串（读段），这些字符串是基于有限字母表集合的，对于DNA序列来说，这个集合是A、C、G和T。

它允许具有共同前缀的读段使用那个前缀，并且只存储读段的末端以表示不同的序列。字典树由节点和边表示。根节点是空的，然后字典树中的每个节点代表一个单独的字符。在DNA的情况下，一个节点可能拥有的最大子节点数为四个。节点的子节点按字母顺序排列。叶节点是代表读段最后一个字符的节点。一个读段由从根节点到叶节点的路径表示。字典树数据结构本身并不适合用于索引基因组序列，因为它存储的是一组字符串。然而，字典树的广义概念被用于后缀树中，后者广泛用于索引参考基因组序列。

1.3 后缀树

后缀树作为字典树数据结构的一种泛化，主要用于模式匹配和在给定的字符字符串中查找子串。它构建为键值对（如Python字典），其中参考序列的所有可能后缀是键，而它们在序列中的位置（索引）是值。

1995年，Esko Ukkonen 提出了一种称为Ukkonen算法的线性时间算法，该算法以线性时间复杂度构建后缀树，即索引所需的时间随着节点数量的增加而线性增加O(n)。为了简单起见，我们将尝试向您展示如何为参考序列构建后缀树以及如何进行读取映射。假设我们的参考序列由10个碱基组成，即“CTTGGCTGGA$”，其中碱基的位置分别是0、1、2、…、9，而$是一个空的尾部位置。

构建后缀树的第一步是从序列中形成后缀（键）和索引（值）对。我们从序列的最后一个字符开始，即位于序列第10位的空字符“$”。然后，我们形成后缀“A$”，其第一个字符位于序列的第9位。我们继续这样做，直到创建所有可能的后缀，如下所示：

请注意，每一行都包含一个后缀（键）和一个位置（值）。然后，从这些键值对中，我们可以构建一个由节点和边组成的树。位置（值）将成为节点，而后缀（键）将成为树的边。如图所示构建后缀树，从顶部的第一个后缀开始，并向下移动以创建分支节点和边，以避免重复常见的字符。

这样，我们将构建一个由节点和边组成的后缀树。整个参考基因组可以被划分为后缀，并以这种方式存储，同时在未分枝的节点中包含后缀和索引位置，因此查找模式或读取在参考基因组中的位置将变得容易。一旦使用后缀树对参考序列进行索引，就可以使用多种搜索算法之一来找到读取映射的位置。例如，要在图2.5中查找“TGG”，我们将从根节点开始搜索“T”，然后从下一个节点开始搜索“GG”；因此，由于有两个叶节点的索引为2和6，这意味着“TGG”与参考序列在位置2和6对齐。

1.4 后缀数组（SA）

后缀数组（SA）与后缀树在模式匹配和查找子串（读取）方面类似，也可以从后缀树构建。它基本上是给定序列所有后缀的排序数组。有多种算法可以用于构建后缀数组，这些算法已由软件包实现。在其最简单的形式中，后缀数组是为一个序列或字符串构建的。对于长度为N的字符串（序列）S，其字符（碱基）索引为1、2、3、…、N，我们可以构建包含序列所有后缀或子串的数组，如S[1…N]、S[2...N]、…、S[N...N]，然后按字典顺序对后缀进行排序。例如，对于序列S = “CTTGGCTGGA$”，我们可以构建后缀数组并按图所示进行排序。

上图展示了如何从排序后的后缀构建后缀数组。数字是序列中的位置。

后缀的排序使得以相同字符串开头的后缀连续出现，这在我们尝试查找读取（子串）的确切匹配时允许快速查找。例如，可以通过跳转到以“TGG”开头的排序后缀来找到读取“TGG”的确切匹配，并且我们可以快速定位位置7和3（坐标从1开始，而不是像后缀树那样从0开始）。当使用后缀数组对参考基因组进行索引时，查找模式或读取的位置将具有线性时间复杂度。使用后缀数组的一个主要缺点是它们需要大量的内存存储，这取决于所使用的基因组的大小。STAR是使用后缀数组的比对器的一个例子。

1.5 Burrows-Wheeler变换

Burrows-Wheeler变换（简称BWT）是一种数据结构算法，它将一个字符串（序列）转换成可压缩的形式，从而实现快速搜索。

BWT被BWA、Bowtie等比对工具所使用。

此外，BWT算法也被Unix和Linux系统用于数据压缩，例如作为bzip2压缩工具。设S是一个长度为n的字符序列或字符串。对于基因组序列而言，它由A、C、G、T或N（表示不确定的碱基）组成。另设$为一个特殊的单个字符，作为序列S的尾部空结束符（例如，S = “CTTGGCTGGA$”）。 Burrows-Wheeler变换的结果是通过对序列生成循环移位得到的，如图的第一列所示。

然后将这些循环移位按字典序排序，形成一个矩阵，称为BWM（Burrows-Wheeler矩阵），如图的第二列所示。BWM中最后一列字符即为序列的BWT。如图所示，序列S的BWT是“AGG$GGTTCTC”。

当比对工具使用BWT时，它会将整个基因组序列转换成BWT形式。然而，使用上述朴素方法计算BWT在计算上非常昂贵，时间复杂度为O(n^2 log n)。（如果还是不理解可以看下图）

相反，我们可以使用后缀数组在O(n)的时间复杂度内计算参考基因组的BWT。

这是使用BWT的比对工具所采用的策略。从后缀数组构建BWT非常简单直接。要获取字符串S的BWT中的第i个字符（其字符索引为i=1, 2, ..., n），只需使用以下规则：

现在，我们从BWT(S)[11]开始推断BWT。

应用上述规则，BWT(S)[11] = S[A[11] - 1] = S[10] = A。

同样地，BWT(S)[10] = S[A[10] - 1] = S[9] = G；

BWT(S)[6] = S[A[6] - 1] = S[5] = G；

但BWT(S)[1]，A[1]小于1；因此，BWT(S)[1] = $。

对于BWT(S)[i]对应于A=9, 5, 8, 4, 7, 3, 2，我们将继续使用BWT(S)[i] = S[A[i] - 1]。

BWT如图所示。

接下来的问题是，为什么我们需要将参考基因组的序列转换为BWT？

BWT有两个目的：首先，BWT将序列的字符分组，使得单个字符在一行中多次出现，因为列是按字母顺序排序的，这可以用于序列压缩，从而减少内存存储。在这个意义上，BWT是可压缩和可逆的。第二个目的是BWT是一种可用于对参考基因组的序列进行索引的数据结构，以便快速找到读取位置。 BWT的特性是我们可以使用最后一个到第一个列映射属性（简称LF映射）将其逆转以获得BWM，其中L是排序后的后缀中对应的BWT的最后一列，而F是排序后的旋转中的第一列。

L和F列如图所示分开。我们只需要第一列（F）和最后一列（L）来将BWT逆转回原始的排序后缀。BWT可以使用LF映射恢复。由于F必须以“$”开始，我们知道这个字符是原始字符串中的最后一个。

"$"前面的第一个字符是L中的“A”；因此，我们现在可以推断出原始序列中的最后两个字符是“A$”。为了推断出“A”之前的字符，让我们找到F中的第一个“A”，它在第二行，所以L中的第二行字符“G”是“A$”之前的字符；因此，我们现在可以推断出原始字符串的最后三个字符是“GA$”。

由于这个“G”是L中的第一个“G”，它指向F中的第一个“G”（第5行）；因此，“GA$”之前的字符是L中的第5行字符，即“G”。所以现在推断出的最后四个字符是“GGA$”。

由于这个“G”是L中的第三个“G”，它指向F中的第三个“G”，即第7行，而L中的第7行字符“T”是“GGA$”之前的字符，因此推断出的最后五个字符是“TGGA$”。这个“T”的排名是L中的第一个“T”，它指向F中的第9行的第一个“T”，所以“TGGA$”之前的字符是L中的第9行的“C”。

因此，推断出的原始字符串的最后6个字符是“CTGGA$”。我们可以继续使用这种L到F的映射，直到恢复整个序列“CTTGGCTGGA$”，如图中的箭头所示。

这个逆转过程被称为反向搜索机制，因为它从字符串的末端开始并向前移动。通常，计算机算法执行此任务。我们可以使用相同的LF映射属性以不同的方式恢复原始字符串。给定LF矩阵，字符串S从字符串的末尾开始重建，即$，这也是F列的第一个字符。$之前的字符将是L列的第一个字符。一般来说，可以使用图所示的L和F的关系推断出之前的字符。

首先，我们通过从L列和F列创建一个两列矩阵并排序来恢复BWM的前两列，如图a所示。第三列通过从L列和已恢复的前两列创建一个三列矩阵并排序来恢复（b）。同样的方法用于恢复排序旋转的第四列，如图c所示。我们将继续这样做，直到恢复所有列，BWM将如图d所示。原始字符串是以BWM中字符“$”结尾的字符串，如图中阴影行所示。

图中的步骤展示了BWT如何可逆。

然而，有多种软件实现的算法可以使用LF映射来逆转BWT。BWT作为数据结构，通过将整个基因组转换为BWT来索引基因组（压缩或未压缩）。一旦基因组被转换，就有多种搜索算法和方法可以在参考基因组上找到子串（读取）的位置。

1.6 FM-索引

FM-索引（全文本分钟空间索引）在BWT的基础上执行反向搜索机制，以线性计算时间复杂度O(N)查找字符串上的精确模式匹配，无论字符串的长度如何。此外，它可以实现高压缩比，使整个人类基因组的索引存储在小于2.0G的内存中。

FN-索引搜索是通过使用LF映射矩阵来完成的。使用FM-索引搜索读取（子串）采用反向匹配，其中反复使用LF映射，直到找到子串的匹配项。FM-索引引入了排名表和查找表。排名表是从BWT（列L）生成的，使用字符的排名，定义为字符在前缀L[1…k]中出现的次数（即，列出每个唯一字符的出现次数和顺序）。排名表充当函数B(c, k)，解释如下：

如果给定一个字符c和长度为k的前缀，那么我们可以推断出该字符在前缀中出现的次数。查找表列出了从排序矩阵的第一列（列F）中每个字符c第一次出现的索引。查找表中字符c的第一次出现可以描述为C。下图显示了Burrows-Wheeler矩阵（a）、排名表（b）和查找表（c）。