1. Bitap Algorithm
论文第 2.3 节 Bitap Algorithm 介绍了 GenASM 所基于的核心近似字符串匹配算法——Bitap。下面根据论文内容,用中文解释该算法的原理与流程。
Bitap 算法的目标
Bitap 算法用于在给定参考文本(如参考基因组)和查询模式(如 read)之间,计算允许最多 个错误的最小编辑距离。
- 当 时,算法退化为精确字符串匹配。
- 当 时,允许替换、插入、删除三类错误,最终找到匹配位置及对应的最小编辑距离。
Bitap 的核心思想是位并行:用二进制位向量表示模式与文本的部分匹配状态,通过简单的位运算(移位、与、或)同时处理多个位置,因此非常适合硬件加速。
算法流程(对应 Algorithm 1)
1. 预处理:生成模式位掩码(Pattern Bitmask, PM)
- 将查询模式 转换为 位的位掩码,其中 是模式长度。
- 对字母表中的每个字符 (例如 DNA 中的 A、C、G、T),生成一个位掩码 。
- 规则:0 表示匹配。因此,当 时,令 ;否则为 1。
- 这样,每个字符的位掩码就编码了该字符在模式中出现的位置。
2. 初始化状态位向量
- 定义 个状态位向量 ,其中 ,表示允许最多 个错误时的部分匹配信息。
- 初始时,所有 的每一位都置为 1(表示尚未匹配)。
- 同时,用 保存上一轮文本迭代的状态,以便在计算中考虑部分匹配。
3. 遍历文本,逐字符计算
算法从文本的最后一个字符开始,向前遍历每个文本字符(算法中 for i in (n-1):-1:0),每次处理一个字符。
对于当前文本字符 curChar:
-
取出该字符对应的模式位掩码 。
-
计算精确匹配状态 :
这里左移一位表示消耗一个模式字符,或上 表示当前文本字符与模式字符的匹配情况。
-
对于每个编辑距离 ,计算四种中间位向量: 1.) 代表上一个文本位置(因为算法是从文本末尾向前遍历, 是文本位置 的状态)。使用 意味着消耗了一个文本字符。
2). 代表当前文本位置(文本位置 的状态)。使用 意味着没有消耗文本字符(文本指针停在原地)。
- 删除(Deletion, D):当前模式字符缺失,只消耗文本字符,不消耗模式字符,错误数加 1。
文本位置 的个错误状态,到达状态,需由当前文本位置 的 转换而来
- 替换(Substitution, S):当前模式字符与文本字符不匹配,同时消耗文本和模式字符,错误数加 1。
文本位置 的个错误状态, 消耗模式字符,到达状态,需由当前文本位置 的 转换而来
- 插入(Insertion, I):当前文本字符缺失,只消耗模式字符,不消耗文本字符,错误数加 1。
文本位置 的个错误状态, 消耗模式字符,到达状态,需由当前文本位置 的 转换而来
- 匹配(Match, M):当前模式字符与文本字符匹配,同时消耗文本和模式字符,错误数不变。
文本位置 的个错误状态, 消耗模式字符, 因错误数不变,所以是第个错误状态
然后,将这四个中间位向量按位与(AND),得到当前编辑距离 的状态位向量:
按位与的作用是保留所有可能的部分匹配位置(即所有位中为 0 的位置),表示在允许 个错误时,从当前位置开始可能存在的匹配。
- 检测匹配:
- 检查每个 的最高有效位(MSB)。如果某一位为 0,则说明在文本位置 处找到了一个匹配,其编辑距离为 。
- 记录起始位置 和编辑距离 。
- 算法继续遍历,直到所有文本位置都被检查完。
关键点说明
- 位向量的含义: 中的每一位对应模式中的一个位置。某位为 0 表示从该位对应的模式位置开始,到当前文本位置为止,已经找到了一个允许 个错误的匹配。
- 从后往前遍历:Bitap 通常从文本末尾向前扫描,这样左移操作对应模式字符的消耗,最终 MSB 为 0 表示整个模式匹配完成。
- 并行性:所有位运算都是按位操作,一条指令可以同时处理 个位置,因此计算效率高。
- 复杂度:时间复杂度为 次位运算,其中 是文本长度, 是最大编辑距离。但原始 Bitap 受限于机器字长,模式长度不能超过一个字(如 64 位)。
原始 Bitap 的局限性(论文 3.1 节)
论文指出原始 Bitap 存在五个主要限制,导致其难以直接用于基因组分析:
- 不支持长读长:位向量长度受机器字长限制。
- 迭代间数据依赖: 依赖上一轮文本迭代和同一轮内前一个 的结果,无法直接并行。
- 不支持回溯:只能找到匹配位置和编辑距离,无法生成 CIGAR 字符串等最优比对结果。
- 计算并行度有限:CPU 核心数有限。
- 内存带宽有限:GPU 实现受限于每个线程的可用带宽。
GenASM 正是针对这些限制,对 Bitap 进行了算法修改(支持多字、消除循环依赖、增加文本级并行)并设计了专用硬件加速器。