BWA

Burrows-Wheeler Aligner作为一个十分出名的alignment软件，在各个生信领域的比对中都发挥了重要的作用。但是如果只是会使用这个软件而不能理解其中的原理的话，那就十分的遗憾了。

今天这一篇主要讲后缀树的相关知识，网上的文章很多，我除了使用自己的理解以外，也会参考别的文章的内容，相互佐证以增加可读性。但是网上的文章中使用到python进行后缀树实现的并不多，一来是应用的实际作用不大，二是本就有不少的后缀树相关的python模块，相信直接使用那些模块会比自己写好很多。但是由于自己从零实现对算法的理解帮助很大，我最后还是自己从零开始写了一个简单的后缀树基础模块。

后缀树的前前身---Trie树

以下这张图就是Trie树。由{bear, bell, bid, bull, buy, sell, stock, stop}几个单词生成的一个Trie树，Trie树本身就具有十分显著的优点，当在一个构建好的树中进行查询的时候，可以看出仅需要线性复杂度的时间。当然缺点也是十分明显的，空间复杂度是其比较大的问题。构建Trie树并不会消耗太多的时间。

特点

1. 根节点无字符
2. 从根节点到各个叶节点的唯一路径的组合，即构建时的字符串
3. 任意节点的子节点所包含的字符都不相同。
4. 每一个节点只包含一个字符 (该特点并不强制，看后续压缩Trie树可知)

可以想象得到，如何在构建好的Trie树中的搜索，也可以直观的从肉眼看出规律。

Trie的搜索过程

1. 从根节点出发。(定义父节点为Node1，搜寻字符串为Sstr)
2. 寻找搜寻字符第一个字符对应的子节点 (寻找Node1的子节点中与Sstr[0]相等的节点，并将其赋值给Node1)
3. 从该子节点继续往下，搜寻字符第二个字符对应的子节点的子节点(继续寻找Node1的子节点与Sstr[1]相等的节点...并...)
4. 重复至找不到对应节点或者到达Sstr的末尾。

Trie树

后缀树的前身---压缩Trie树

从Trie树继续靠近后缀树的过程中，还存在一个压缩Trie树。顾名思义，压缩Trie树是在Trie树的基础上进行压缩。那么何为压缩，并且压缩的意义在哪？

Trie树的压缩
除根节点外，若任意节点的父节点只有其一个子节点。(即 若该节点为独生子女)，那儿我们将其和父节点压缩为一个节点

压缩的过程是一个精简Trie树的分支的过程，但也由此带来了一些代码上实现的困难。例如任意一个节点不一定为单个字符，即任意一个节点均可能含有多个字符了(此处强调的是节点含多字符，父节点含单字符，父父节点仍可以是多字符。)

压缩Trie树

后缀树的最后一步---什么是后缀

后缀其实指的是一个字符的后缀集合。
例如，

对于abcabxabcd，则可以生成这样的后缀集合
abcabxabcd
bcabxabcd
cabxabcd
abxabcd
bxabcd
xabcd
abcd
bcd
cd
d

若将这些后缀集合左对齐的写在每一行，就可以看到这样的一个倒三角的结构。每次删掉第一个字符，从而生成一个含等同于 字符串长度 数量的字符的后缀集合(即含有10行，该字符串长度为10)

后缀树的定义

在看完上面的几个前要后，我们可以来提出后缀树的基本定义了。

后缀树是一个由单个字符串的后缀集合所构建的压缩Trie树。

后缀树的构建

后缀树本身与压缩Trie树是没什么大的区别的，完全可以使用一个（压缩后缀树的构建算法+生成后缀集合）的算法去生成一个后缀树，但这样的话，时间复杂度太高。
后来，在1995年，Ukkenon发表了论文《On-line construction of suffix trees》。这样，一下子将构建后缀树的时间复杂度降低到了线性的尺度上。

接下来我们就重点在描述Ukkenon的后缀树构建算法，并且在后文提供Python的代码实现。

Ukkenon算法

由于其中部分与我python代码实现相关的地方，可能与原论文的部分表述有所不同。请倍加注意。

一开始为了构建Suffix Tree，我们先初始化一些具有重要作用的存储对象。

active_point = (root, '', 0) # 分别称为 (父节点，指示子节点，位置索引)
remainder = 0 # 称为 剩余后缀

现在无法理解很正常，之后再讲述算法的过程中会渐渐讲解其意义和作用。
我们接下来从构建字符串abcabxabcd的后缀树的实际过程中讲解这个算法

高层次的来看构建过程

1. 从左到右，每次只向后缀树加入一个字符。 (注意：用的是加入，不是插入。) （加入的是一个字符，但实际上操作的是每一个以该字符开头的后缀。)
2. 加入树前，会先确认 该字符是否在某已有字符串的 前缀中 （前缀与后缀类似，均为一个对应字符串的集合，确认的目的是，若在，则说明该字符开头的后缀需要与别的后缀共享一些节点。）

第一个字符，a

由于本身这个是个空的后缀树，那么直接在根节点后插入一个节点，a即可。

active_point = (root, '', 0),remainder = 0

1. 在插入前，remainder= remainder + 1，表明，我们现在准备要插入一个后缀。
2. 插入a节点到父节点root上。
3. 插入后，由于我们新增了一个叶节点，且完成，所以remainder= remainder - 1
active_point = (root, '', 0),remainder = 0

第二个字符，b

1. 在插入前，remainder= remainder + 1，表明，我们现在准备要插入一个后缀。
2. 扩展每一个叶节点，所以a节点-->ab节点
3. 查询父节点下，(所有已插入字符串的前缀中不包含b)，所以，b需要独立建立新的一个节点，所以向父节点插入b节点 ，注意的是，父节点下查询时只有ab节点，但ab节点的前缀为[a,ab]。
4. 插入后，由于我们新增了一个叶节点，且完成，所以remainder= remainder - 1
active_point = (root, '', 0),remainder = 0

第三个字符，c

1. 在插入前，remainder= remainder + 1，表明，我们现在准备要插入一个后缀。
2. 扩展每一个叶节点，所以ab节点-->abc节点，b节点-->bc节点。
3. 查询父节点下，(所有已插入字符串的前缀中不包含c)，所以，c需要独立建立新的一个节点，所以向父节点插入c节点
4. 插入后，由于我们新增了一个叶节点，且完成，所以remainder= remainder - 1
active_point = (root, '', 0),remainder = 0

第四个字符，a

1. 在插入前，remainder= remainder + 1，表明，我们现在准备要插入一个后缀。
2. 扩展每一个叶节点，所以abc节点-->abca节点，bc节点-->bca节点，c节点-->ca节点。
3. 查询父节点下，所有已插入字符串的前缀中是否包含a？发现的确是有的，且“最后一个”节点为，就是abca节点。

特殊操作1
active_point，由于该节点abca的父节点还是root，所以第一个不变。
active_point[1] = 'a'，指示子节点则等于，当前字符a。
active_point[2] = '1'，位置索引，因为a在abca的第一个位置，所以等于1。

3. 插入后，由于我们无新增了一个叶节点，且完成，所以remainder保持不变。
   active_point = (root, 'a', 1),remainder = 1

为什么我们这里要这样操作？

因为一旦我们发现了已插入字符串的前缀就是我们要插入的后缀，那么我们要找到“最后一个”节点。（为啥要最后一个？因为当字符串更长时，会出现你找到节点有很多。这时我们需要最后一个，可见第九个字符步骤）。那么说明这个将要插入的后缀与这“最后一个”节点的是共享一些字符串的。那么就会产生一次分叉(这分叉可能已经存在)，那么这个分叉还不能直接分叉，因为还不能确定是在这个地方分叉，一定要到不一样的地方才是分叉点，如果一直都一样，就是完全一个重复的后缀。

所以这个地方我们插入了我们准备插入的后缀吗？
没有，现在整个后缀树仍然只有后缀集合里的三个后缀。

第五个字符，b

1. 在插入前，remainder= remainder + 1，表明，我们现在准备要插入一个后缀。
2. 扩展每一个叶节点，所以abca节点-->abcab节点，bca节点-->bcab节点，ca节点-->cab节点。
3. 查询父节点下，所有已插入字符串的前缀中是否包含ab？发现的确是有的，且“最后一个”节点为，就是abcab节点。（由于上面那个后缀我们没插入，现在它也被扩展为ab了。）

特殊操作1
active_point，由于该节点abcab的父节点还是root，所以第一个不变。
active_point[1] = 'ab'，指示子节点则等于，扩展后的字符ab。
active_point[2] = '2'，位置索引，因为b在abcab的第一个位置，所以等于1。

3. 插入后，由于我们无新增了一个叶节点，且完成，所以remainder保持不变。
   active_point = (root, 'ab', 2),remainder = 2

第六个字符，x

1. 在插入前，remainder= remainder + 1，表明，我们现在准备要插入一个后缀。
2. 扩展每一个叶节点，所以abcab节点-->abcabx节点，bcab节点-->bcabx节点，cab节点-->cabx节点。
3. 查询父节点下，所有已插入字符串的前缀中是否包含abx？发现没有。所以这是一个分叉点（分叉点的解释可见第四个字符的解释文字。）
   

特殊操作2
active_point[0]，从父节点开始找，其子节点们。
active_point[1] = 'abx'，在子节点们中找其前缀中含abx的节点。(其实就是第五个字符时找到的abcab，但现在为abcabx)
active_point[2] = '2'，使用位置索引来寻找到将分叉的位置，即上图中绿色箭头所指向的位置

1. 则我们插入abx这个节点，即如下图

   
   
2. 由于新增了一个叶节点，所以remainder= remainder - 1，但是由于remainder还>1，则我们继续。
3. 插入bx这个节点(为什么还要插入bx？，因为我们中间漏了bx的后缀，所以得补上。)

一样的，也是从active_point的父节点开始寻找，active_point[1]=bx的节点，且active_point[2]由于父节点不变，active_point[1]变化了。所以我们得-1。同样为新增了叶节点，所以remainder= remainder - 1

4. 插入x这个节点，此时，active_point = (root, 'x', 0),remainder = 1

此时，由于active_point[2]已经为0，说明，我们要在父节点的右侧直接插入一个节点。插入后，同样为新增了叶节点，所以remainder= remainder - 1。

自此，遇到分叉点的情况已经解决。
active_point = (root, '', 0),remainder = 0

Q：新增的两个分叉点有关系吗？
A：即ab和b节点，创建时就是因为存在abx到bx的前后关系。即如果我们再遇到abk之类的，我们插入完abk后，必然也要插入bk节点，所以一定会从ab节点到b节点。如果我们每次都要重新再找一次b节点，会使速度变慢。
Q：能利用其关系来加速吗？
A：能。加入后缀链接，即图中的绿色虚线箭头(有向)

第七个字符，a

同第四个字符。

active_point = (root, 'a', 1),remainder = 1

第八个字符，b

同第五个字符

active_point = (root, 'ab', 2),remainder = 2

第九个字符，c

此时遇到一个比较神奇的事情，我们需要寻找已插入字符串时，发现abc存在。但“最后一个”节点不能cover整个字符串。那我们先不管，继续下面的步骤。

1. 在插入前，remainder= remainder + 1，表明，我们现在准备要插入一个后缀。
2. 扩展每一个叶节点。
3. 查询abc，找到“最后一个”节点是c。

特殊操作1
active_point，由于该节点c的父节点是ab，所以active_point[0] = 节点ab
active_point[1] = 'abc'，指示子节点则等于，扩展后的字符abc。
active_point[2] = '1'，位置索引，因为c在cabxabc的第一个位置，所以等于1。

3. 插入后，由于我们无新增了一个叶节点，且完成，所以remainder保持不变。
   active_point = (节点ab, 'abc', 1),remainder = 3

第十个字符，d

那么就到了最后一个字符了，这个字符中也会使用到前文说的后缀连接。

1. 在插入前，remainder= remainder + 1，表明，我们现在准备要插入一个后缀。
2. 扩展每一个叶节点
3. 查询父节点下，所有已插入字符串的前缀中是否包含abcd？发现没有。所以这是一个分叉点（分叉点的解释可见第四个字符的解释文字。）
   

特殊操作2
active_point[0]，从父节点开始找，即节点ab，找其子节点们。
active_point[1] = 'abcd'，在子节点们对应字符串找其前缀中含abc的节点。(其实就是第九个字符时找到的cabxabc，但现在为cabxabcd)(虽然cabxabc自己不含有abc的前缀，但你要考虑到其对应字符串为abcabxabc。)
active_point[2] = '1'，使用位置索引来寻找到将分叉的位置，即上图中绿色箭头所指向的位置

1. 则我们插入abcd这个节点，即分裂cabxabcd，分成c -->[abxabcd,d]
2. 由于新增了一个叶节点，所以remainder= remainder - 1，但是由于remainder还>1，则我们继续。
3. 由于active_point[0]存在后缀连接，所以我们要沿着后缀连接的方向改变，active_point[0] = 节点b，active_point[1] ='bcd'，但由于父节点的变化，不涉及根节点，所以active_point[2]不变
4. 插入bcd这个节点

一样的，也是从active_point的父节点，即节点b，开始寻找，active_point[1]=bcd的“最后一个”节点。所以找到了cabxabcd，所以分裂cabxabcd，，分成c -->[abxabcd,d]。同样新增了叶节点，所以remainder= remainder - 1

4. 因为此时的active_point[0]无后缀连接，所以我们把active_point[0]重置为根节点，active_point[1] = cd，active_point[2]-=1。由于父节点重置为根节点，所以-1。
5. 插入cd这个节点

一样的，也是从active_point的父节点，即节点b，开始寻找，active_point[1]=bcd的“最后一个”节点。所以找到了cabxabcd，所以分裂cabxabcd，，分成c -->[abxabcd,d]。同样新增了叶节点，所以remainder= remainder - 1

6. 插入d这个节点，此时，active_point = (root, 'd', 0),remainder = 1

此时，由于active_point[2]已经为0，说明，我们要在父节点的右侧直接插入一个节点。插入后，同样为新增了叶节点，所以remainder= remainder - 1。

最后的后缀树长这样。

image.png

回头总结

1. active_point的变化规律

active_point[0]，只有当遇到分叉点时，我们把最后能找到的那个节点的父节点作为active_point[0]。
active_point[1]，原作者跟我在这个地方有点差异，我会选择在这个位置储存 已存在的字符的。
active_point[2]，原则上等于最后一个字符(非分叉点)在能找到的那个节点上的位置。
active_point[2]的变化，当active_point[0]为根节点时，每次插入完要-1，当active_point[0]不是根节点，且其没有后缀连接时，每次-1

2. remainder则是观察有没有叶节点的变化，因为每个叶节点代表一个字符串，而且是唯一的一个字符串。若叶节点增加了，则代表加入了一个新的字符串。则remainder-=1

全文结束。。。。。。
写完大概是个废鱼了。。。对了。。。还有代码。。。还有些注意事项

注意事项

1. 后缀树的查找时存在一种情况，若上述的后缀树不以d结尾，会发生什么情况呢？

会使得很多的分支都不会出现在图上，也就是abc，bc，c这几个后缀，变成了一个隐式的字符串在后缀树中。(因为我们一直存着，但是一直遇不到分叉点。)那么这种情况只要简单的加一个结尾即可，有时会加入特殊字符。

代码在这...
algorithms_in_python/suffix_tree/trie_tree_Ukkonen_al.py

其中，另一个trie_tree是trie的实现，并且构建树的方法也完全不一样。而trie_tree_Ukkonen_al才是使用Ukkonen构建的后缀树。

里面如果有些奇怪的注释之类大家可以忽略。。。
还有一些奇怪的命名大家也可以忽略。。。

写完大概是个废鱼了，继续写点别的好了。。。
如果有人问这个东西能做啥的话。

1.查找字符串O是否在字符串S中。
方案：用S构造后缀树，按在trie中搜索字串的方法搜索O即可。
原理：若O在S中，则O必然是S的某个后缀的前缀。
例如：leconte，查找O：con是否在S中，则O（con）必然是S（leconte）的前缀。
2.指定字符串T在字符串S中的重复次数。
方案：用S+’$’构造后缀树，搜索T节点下的叶子节点数目即为重复次数
原理：如果T在S中重复了两次，则S应有两个后缀以T为前缀，重复次数自然统计出来了。
3.字符串S中的最长重复子串
方案：原理同2，具体做法是找到最深的非叶子节点。
这个深指从root所经历过的字符个数，最深非叶子节点所经历的字符串起来就是最长重复子串。为什么非要是叶子节点呢？因为既然是要重复的，当然叶子节点个数要>=2
4.两个字符串S1,S2的最长公共子串（而非以前所说的最长公共子序列，因为子序列是不连续的，而子串是连续的。）
方案：将S1#S2$作为字符串压入后缀树，找到最深的非叶子节点，且该节点的叶子节点既有#也有$.
5.最长回文子串

以上是我copy and paste来的。。。
。。。

参考

1. http://www.cnblogs.com/gaochundong/p/suffix_tree.html
2. https://en.wikipedia.org/wiki/Suffix_tree
3. http://www.cnblogs.com/aiyelinglong/archive/2012/04/09/2439777.html