0

AATGGGTTCCA 在多个序列中寻找字符串序列。但希望允许 N 编辑距离（替换或插入字符串）

所以在一个长字符串中它可能匹配 AATG* C *GTTCCA （替换）或 AATGGTTCCA （删除）或 AATGGG* T *TTCCA （插入）

处理许多序列时最快的算法是什么。

编辑：因为你可以匹配无限的字符串..让我们假设最多 5 个位置的 1 个核苷酸插入、删除或替换 BLAST 也是一个选项

4

1 回答 1

5

一种可能的策略是预处理您正在搜索的子序列，然后将结果模式应用于所有序列。

预处理将产生最大Levenshtein 距离为N的所有模式。缺点是这需要相当大的空间：长度为s的字符串将产生多达s N+1 个可能的模式。要从可能性中获取正则表达式，您可以$re = join "|", @possibilities. 由于此处将使用 trie 优化，因此生成的正则表达式应该非常快。

获得可能性的例子

因为大小可以非常快速地增长，所以这里有一个示例来获取所有具有 Levenshtein 距离 1 到 的字符串AC。我们可以通过

my@actg=qw/A C T G/;

sub ld{
    my $distance = shift;
    ($distance and @_) or return [@_];
    my $car = shift;
    my @unchanged   = map [$car, @$_], ld($distance, @_);
    my @inserted    = map { my $ins = $_; map [$ins, $car, @$_], ld($distance-1, @_) } @actg;
    my @substituted = map { my $sub = $_; map [      $sub, @$_], ld($distance-1, @_) } @actg;
    my @deleted     = ld($distance-1, @_);
    return @unchanged, @inserted, @substituted, @deleted;
}

如您所见，此代码尚未优化。这可以通过记忆和分解常见ld($distance-1, @_)调用来改进。它还会产生不必要的重复。

然后我们可以打印出所有独特的可能性，比如

my %uniq;
$uniq{$_} = undef for map {join "", @$_} ld(1, split//, "AC");
my @possibilities = keys %uniq;
say for sort @possibilities;

输出：

A
AA
AAC
AC
ACC
AG
AGC
AT
ATC
C
CAC
CC
GAC
GC
TAC
TC

编辑：优化

我对发布上述如此低性能的代码感到内疚。这是一个优化的版本，它更喜欢字符串而不是数组，并采用了记忆：

use constant ACTG => qw/A C T G/;
use List::MoreUtils 'uniq';

sub ld2 {
    my ($d, $t) = @_;
    ($d and length $t) or return $t;
    state $cache = {};
    my $loc = $t;
    my $result = $cache->{$loc}[$d] //= do {
        my $c = substr $t, 0, 1, '';
        my @unchanged   = map $c . $_, length $t ? ld2($d, $t) : $t;
        my @changed     = map {
            my $s = $_;
            map {; ($s . $c . $_), ($s . $_),  ($_) } length $t ? ld2($d-1, $t) : $t;
        } ACTG;
        [ uniq @unchanged, @changed ];
    };
    return wantarray ? @$result : $result;
}

像join "|", ld2($N, $string). 它运行得更快，并且应该使用更少的内存。

于 2013-06-28T14:27:46.910 回答