2.3.1 BWT-FM-Index
下面用一个具体的 DNA 序列例子,解释线性 BWT 和 FM-Index 是怎么工作的。参考序列选 GATTACA,末尾加结束符 $,变成 GATTACA$。
1. BWT 的构建
BWT(Burrows-Wheeler Transform)的核心是:把所有循环旋转排序,然后取最后一列。
1.1 生成所有循环旋转
原字符串:G A T T A C A $,索引 0 到 7。
循环旋转就是每次把最后一个字符移到最前面,得到 8 个字符串:
| 原始索引 | 旋转字符串 |
|---|---|
| 0 | G A T T A C A $ |
| 1 | A T T A C A $ G |
| 2 | T T A C A $ G A |
| 3 | T A C A $ G A T |
| 4 | A C A $ G A T T |
| 5 | C A $ G A T T A |
| 6 | A $ G A T T A C |
| 7 | $ G A T T A C A |
1.2 按字典序排序
字符顺序:$ < A < C < G < T。排序后得到:
| 行号 | 旋转字符串 | 原始索引 |
|---|---|---|
| 0 | $ G A T T A C A | 7 |
| 1 | A $ G A T T A C | 6 |
| 2 | A C A $ G A T T | 4 |
| 3 | A T T A C A $ G | 1 |
| 4 | C A $ G A T T A | 5 |
| 5 | G A T T A C A $ | 0 |
| 6 | T A C A $ G A T | 3 |
| 7 | T T A C A $ G A | 2 |
1.3 取最后一列,得到 BWT
把每一行的最后一个字符拿出来:
| 行号 | 最后一列字符 |
|---|---|
| 0 | A |
| 1 | C |
| 2 | T |
| 3 | G |
| 4 | A |
| 5 | $ |
| 6 | T |
| 7 | A |
所以 BWT = A C T G A $ T A,写成一串就是 ACTGA$TA。
BWT 的特点是:相同字符会聚集在一起,比如 A 集中在几个位置,便于压缩和索引。
2. FM-Index 的两个关键表
FM-Index 在 BWT 基础上,增加两个表:C 表 和 Occ 表。
2.1 第一列 F
排序后每一行的第一个字符组成第一列 F:
| 行号 | F 列 |
|---|---|
| 0 | $ |
| 1 | A |
| 2 | A |
| 3 | A |
| 4 | C |
| 5 | G |
| 6 | T |
| 7 | T |
2.2 C 表
C 表记录:每个字符在 F 列中第一次出现的位置。
$第一次出现在行 0,所以C[$] = 0A第一次出现在行 1,所以C[A] = 1C第一次出现在行 4,所以C[C] = 4G第一次出现在行 5,所以C[G] = 5T第一次出现在行 6,所以C[T] = 6
2.3 Occ 表
Occ 表记录:BWT 中到某个位置为止,某个字符出现了多少次。
BWT = A C T G A $ T A,索引 0 到 7:
| 索引 | BWT 字符 | 到此处 A 的次数 | C 的次数 | G 的次数 | T 的次数 | $ 的次数 |
|---|---|---|---|---|---|---|
| 0 | A | 1 | 0 | 0 | 0 | 0 |
| 1 | C | 1 | 1 | 0 | 0 | 0 |
| 2 | T | 1 | 1 | 0 | 1 | 0 |
| 3 | G | 1 | 1 | 1 | 1 | 0 |
| 4 | A | 2 | 1 | 1 | 1 | 0 |
| 5 | $ | 2 | 1 | 1 | 1 | 1 |
| 6 | T | 2 | 1 | 1 | 2 | 1 |
| 7 | A | 3 | 1 | 1 | 2 | 1 |
3. 用 FM-Index 搜索模式 ACA
假设我们要在参考序列 GATTACA 中找模式 ACA。
FM-Index 的搜索方式是:从模式最后一个字符开始,反向向前搜索。
3.1 第一步:找最后一个字符 A
在 F 列中,A 出现在行 1、2、3。所以初始范围是:
l = 1, r = 3
这表示所有以 A 开头的旋转,行 1 到行 3。
3.2 第二步:向前找 C
当前范围是行 1 到行 3。我们看这些行的 BWT 字符(最后一列):
| 行号 | BWT 字符 |
|---|---|
| 1 | C |
| 2 | T |
| 3 | G |
我们要找前一个字符 C,只有行 1 的 BWT 是 C。
所以行 1 匹配。现在需要通过 LF 映射,找到这个 C 在 F 列中的位置。
F 列中 C 出现在行 4,所以新范围变成:
l = 4, r = 4
3.3 第三步:向前找 A
当前范围是行 4。行 4 的 BWT 字符是 A。
我们要找 A。F 列中 A 出现在行 1、2、3。所以新范围是:
l = 1, r = 3
3.4 结果
最终范围是行 1 到行 3?等等,这里需要更精确的 LF 映射计算。实际上,用标准公式:
- 初始:
l = 1, r = 3(F 列中 A 的范围) - 找
C:l_new = C[C] + Occ(C, l-1) = 4 + Occ(C, 0) = 4 + 0 = 4r_new = C[C] + Occ(C, r) - 1 = 4 + Occ(C, 3) - 1 = 4 + 1 - 1 = 4- 得到
l = 4, r = 4
- 找
A:l_new = C[A] + Occ(A, l-1) = 1 + Occ(A, 3) = 1 + 1 = 2r_new = C[A] + Occ(A, r) - 1 = 1 + Occ(A, 4) - 1 = 1 + 2 - 1 = 2- 得到
l = 2, r = 2
最终范围是行 2 到行 2,也就是行 2。
行 2 对应的旋转是:
A C A $ G A T T
它的原始起始索引是 4。在原始序列 G A T T A C A 中,从索引 4 开始正好是 A C A。
所以 FM-Index 成功找到了模式 ACA,位置在原始序列的索引 4。
4. 总结
- BWT:把参考序列的所有循环旋转排序,取最后一列。它让相同字符聚集,便于压缩和索引。
- FM-Index:在 BWT 上建立 C 表和 Occ 表,支持反向搜索模式。
- 搜索过程:从模式末尾字符开始,在 F 列中找范围,然后利用 LF 映射逐步向前,直到模式开头。
- 例子:参考
GATTACA$,BWT =ACTGA$TA,搜索ACA,最终定位到原始序列索引 4。
线性 BWT 和 FM-Index 是线性参考比对的核心索引结构;而 GBWT 就是把这种思想推广到泛基因组图的单倍型路径上。