rust-wfa

简介

本项目是波前比对算法的实现。该算法首次在2020年发表的一篇文章中被描述： 使用波前算法进行快速间隙斜比对。

以下为摘要部分：

在这篇文章中，我们提出了波前比对算法（WFA），这是一种利用序列之间的同源区域来加速比对过程的精确间隙斜算法。与运行时间为二次方的传统动态规划算法相比，WFA的运行时间为O(ns)，与读长度n和比对得分s成正比，使用O(s²)内存。此外，我们的算法具有简单的数据依赖性，可以被现代编译器的自动功能轻松向量化，适用于不同架构，无需修改代码。我们评估了我们的算法以及其他最先进的实现。结果表明，WFA比其他方法快20-300倍，在比对短Illumina-like序列时，以及使用Oxford Nanopore Technologies等产生的长噪声读数时，速度可快10-100倍。

该算法的参考实现是文章作者用C语言编写的，可在github.com/smarco/WFA2-lib找到。

波前解释

波前基础知识

本说明假定读者熟悉Smith-Waterman-Gotoh 序列比对算法及其间隙斜版本。

波前比对算法使用三个二维矩阵来计算一对字符串 Query 和 Text 之间的最优比对：Insert_{i, j}、Matches_{i, j} 和 Deletes_{i, j}。

与SWG比对不同，i 和 j 并不对应于 Query 或 Text 中的位置。

• i 对应于一个 得分。

• j 对应于动态规划矩阵中将用于对齐 Query 和 Text 的对角线。以下图显示了这些对角线的布局。

每个波前矩阵单元格中存储的值是，对于分数i，对角线j处可以对齐的文本的字符数。

如果我们对齐了n个字符的文本，则对齐的查询字符数等于n加上当前对角线的编号。例如，对于"CAT"和"CAC"的对齐，Matches_{i, j} = 2，由于在0号对角线处，如果匹配/不匹配惩罚/奖励设置为0，则可以匹配最多2个字符("CA")。

我们将插入的方向定义为在查询中对齐一个额外的字符（在对齐矩阵中水平移动）和删除定义为相反的文本（在对齐矩阵中垂直移动）。

在最高层次上，波前对齐算法可以定义为如下（其中pens是一个包含不匹配/间隙惩罚的结构）

    let mut current_front = new_wavefront_state(query, text, pens);
    loop {
        current_front.extend();
        if current_front.is_finished() {
            break;
        }
        current_front.increment_score();
        current_front.next();
    }

匹配扩展

波前对齐的一个关键方面是匹配惩罚设置为0。回想一下矩阵中存储的值的定义

每个矩阵单元格中存储的值是，对于分数i，对角线j处可以对齐的文本的字符数。

对于给定位置，只要文本和查询在SWG矩阵的对角线上的字符匹配，分数就可以扩展。因此，我们可以获得波前扩展函数

for diag in current_diagonals {
	let mut x = matches[current_score][diag];
	while text[x] == query[x + diag] {
		x += 1;
		matches[current_score][diag] = x;
	}
}

波前递归关系

Smith-Waterman-Gotoh定义以下递归关系来构建对齐矩阵

• Deletes_{i, j} = min( Deletes_{i - 1, j} + 间隙扩展惩罚, Matches_{i - 1, j} + 间隙扩展惩罚 + 间隙开放惩罚 )
• Inserts_{i, j} = min( Inserts_{i, j - 1} + 间隙扩展惩罚, Matches_{i, j - 1} + 间隙扩展惩罚 + 间隙开放惩罚 )
• Matches_{i, j} = min( Matches_{i, j} + 匹配/不匹配文本_i到查询_j的惩罚, Inserts_{i, j}, Deletes_{i, j} )

WFA_next函数基于这些关系。WFA_extend函数等价于在DP矩阵的对角线下降，当文本_i == 查询_j。一旦达到该方程不成立的单元格，WFA_next等价于计算该单元格左侧和右侧的DP单元格。

重新表述SWG递归关系以适用于波前

• 对于删除矩阵，在分数i和对角线j处可以匹配的字符数是Deletes _{i - 间隙扩展惩罚, j + 1} 和Matches_{i - 间隙扩展惩罚 - 间隙开放惩罚, j + 1} + 1的最大值。'+'1是因为每个单元格中存储的值的定义。
• 这与插入矩阵相同，只是值将来自对角线j - 1（因为插入在SWG矩阵中是向右移动）。这次，我们不添加1，因为每个单元格中存储的是匹配的文本字符数，而插入是查询中的额外字符。
• 匹配是Deletes_{i, j}、Inserts_{i, j}和Matches_{i - 不匹配惩罚, j}中的最大值。

检查算法是否结束。

最终单元格位于特定的对角线：对角线编号为查询长度减去文本长度。在每次循环中，我们将检查对角线上的匹配字符数是否等于文本的长度。

回溯波前以构建对齐

回溯算法在原始文章中没有指定，我自行推导，因为它并不复杂。一个重要细节是，如果我们处于匹配状态，我们需要撤销波前的扩展。

Rust实现

验证我的实现

验证WFA算法给出的分数与SWG对齐相同。

我在reference.rs中实现了一个朴素、未经优化的SWG对齐版本。我编写了validate.rs（一个二进制目标，编译后生成独立的可执行二进制文件）。

./target/release/validate -h
rust_wfa 0.1.0

USAGE:
    validate [OPTIONS] --min-length <MIN_LENGTH> --max-length <MAX_LENGTH> --min-error <MIN_ERROR> --max-error <MAX_ERROR>

OPTIONS:
    -h, --help                       Print help information
        --max-error <MAX_ERROR>      
        --max-length <MAX_LENGTH>    
        --min-error <MIN_ERROR>      
        --min-length <MIN_LENGTH>    
    -p, --parallel                   
    -V, --version  

该程序生成用户指定的长度区间内的随机字符串，以及第二个经过错误率（百分比）区间变异的字符串版本。然后它使用WFA和SWG算法对这两个字符串进行对齐，并检查它们的分数是否相同（由于最优对齐分数可能有多个对齐方式，所以不对对齐本身进行比较）。它还可以并行运行，在检测到的每个CPU核心上并行执行不同的文本/查询字符串对。

使用此可执行文件修复我算法中的剩余错误后，我现在能够比较成千上万个字符串的对齐，这两种算法之间的对齐分数没有差异，这让我确信我的实现是正确的。

验证对齐及其分数是否匹配。

以类似的方式，我编写了validate_score_matches_alignment。此程序编译成具有相同名称的二进制文件，可以运行以生成字符串对、对它们进行对齐，并验证对齐是否与其分数匹配：分数从对齐中重新计算，然后进行比较。

./target/release/validate_score_matches_alignment -h
rust_wfa 0.1.0
USAGE:
    validate_score_matches_alignment [OPTIONS] --min-length <MIN_LENGTH> --max-length <MAX_LENGTH> --min-error <MIN_ERROR> --max-error <MAX_ERROR>

OPTIONS:
    -h, --help                       Print help information
        --max-error <MAX_ERROR>      
        --max-length <MAX_LENGTH>    
        --min-error <MIN_ERROR>      
        --min-length <MIN_LENGTH>    
    -p, --parallel                   
    -V, --version                    Print version information

这个程序使我能够识别另一个错误（一个错误的变量），它产生了不正确的对齐，分数是正确的（因此没有被之前的验证程序检测到）。我现在可以验证数千个对齐，而不会产生错误。

基准测试

我首先在我的朴素实现上运行了基准测试。结果如下表所示。n = 序列的长度，d = 每个序列之间元素差异的速率。

正如预期的那样，SWG对齐不依赖于错误率。我的朴素实现比原始实现慢得多，并且只有在高度相似的序列中才比SWG对齐快。

这是一个WFA算法的初始、朴素实现。例如，我将波前存储为Vec，这并不高效。

在下一个版本中，我将我的实现重写为使用更高效的1D Vec。

运行时间减少了约90%。我的实现现在与参考实现具有竞争力，并且比高效的SWG更高效。