在 Perl 中高效地匹配子字符串

Question

我正在寻找一种高效的解决方案，以在主字符串中容忍n个不匹配项的情况下找到可能最长的子字符串。

例如： 主字符串

1. AGACGTAC TACTCTACT AGATGCA*TACTCTAC* 2. AGACGTAC TACTCTACT AGATGCA*TACTCTAC* 3. AGACGTAC TACTCTACA AGATGCA*TACTCTAC* 4. AGACGTAC TACTTTACA AGATGCA*TACTCTAC*

搜索字符串： TACTCTACT：应被视为与上述所有主字符串匹配。

另外，有可能出现子字符串的一部分位于主字符串的末尾，我也想将其包括在内。

如果您能提供一些指导意见，我将不胜感激。

注：我将有一个搜索字符串和大约1亿个主字符串来搜索子字符串。

谢谢！ - Abhi

- Abhi

---

问题。1）“主字符串”中的数字是一个“主字符串”中的单独行，还是四个单独的“主字符串”（我假设是后者，但想确认一下）？2）您能解释一下为什么搜索字符串会被认为与＃3和＃4匹配吗？前两个明显匹配，但我不清楚为什么最后两个会匹配，这听起来可能是一个关键的细节需要理解。3）您能解释一下您所说的子字符串的一部分可能在主字符串的末尾的含义是什么吗？ - Brian Gerard

你提到了“最长”，但是你的例子并没有给出任何关于“最长”的线索。你是否指的是最少不匹配的？ - ikegami

@Brian Gerard：1）这些是数百万个主字符串中的4个。2）因为他容忍不匹配。＃3：A对于T，＃4，T对于C，A对于T。3）"... TACT"应该与结束行匹配搜索字符串TACTCTACT，因为它可能后跟“CTACT”。 - ikegami

在您的第四个示例中，“TACTCTAC*”比“TACTTTACA”更匹配。 - ikegami

@ikegami：非常感谢您迄今为止的帮助。我有一个后续问题。以前，我只有一个查询字符串，一次比较数百万个搜索字符串的子字符串出现情况。现在，如果我想将搜索空间扩大到包括多个查询字符串，我可以再次运行相同的逻辑，但我可以使用其他数据结构来构建类似于查询字符串的后缀树，并将其用于在搜索字符串中搜索子字符串吗？谢谢！ - Abhi

Answer 1

我不确定这是否符合你的要求，但是BioPerl有一个名为Bio::Grep::Backend::Agrep的近似匹配工具：

Bio::Grep::Backend::Agrep使用agrep查找查询。

而agrep是“近似grep”。据我所知，它会建立一个数据库，然后使用该数据库使您的搜索速度更快，因此听起来这就是您想要的。

看起来您正在处理DNA序列，因此您可能已经安装了BioPerl。

您还可以尝试直接使用String::Approx：

用于近似匹配（模糊匹配）的Perl扩展。

尽管如此，我认为Bio::Grep::Backend::Agrep会更快、更适合您的任务。

- mu is too short

---

我会尝试这个并让大家知道它是否有效和高效。 - Abhi

Answer 2

use strict;
use warnings;
use feature qw( say );

sub match {
   my ($s, $t, $max_x) = @_;

   my $m = my @s = unpack('(a)*', $s);
   my $n = my @t = unpack('(a)*', $t);

   my @length_at_k     = ( 0 ) x ($m+$n);
   my @mismatches_at_k = ( 0 ) x ($m+$n);
   my $offset = $m;

   my $best_length = 0;
   my @solutions;
   for my $i (0..$m-1) {
      --$offset;

      for my $j (0..$n-1) {
         my $k = $j + $offset;

         if ($s[$i] eq $t[$j]) {
            ++$length_at_k[$k];
         }
         elsif ($length_at_k[$k] > 0 && $mismatches_at_k[$k] < $max_x) {
            ++$length_at_k[$k];
            ++$mismatches_at_k[$k];
         }
         else {
            $length_at_k[$k] = 0;
            $mismatches_at_k[$k] = 0;
         }

         my $length = $length_at_k[$k] + $max_x - $mismatches_at_k[$k];
         $length = $i+1 if $length > $i+1;

         if ($length >= $best_length) {
            if ($length > $best_length) {
               $best_length = $length;
               @solutions = ();
            }

            push @solutions, $i-$length+1;
         }
      }
   }

   return map { substr($s, $_, $best_length) } @solutions;
}

say for match('AABBCC', 'DDBBEE', 2);

输出：

AABB
ABBC
BBCC

- ikegami

---

你的代码完全正常。只是想知道在一定时间内检查的序列数量是否可行，或者是否有方法可以扩展它。我正在搜索一个包含36个字母的子字符串（最多5个不匹配）在9000万个100个字母长的搜索序列中。这个搜索大约需要13个小时。根据你的经验，这听起来正常吗？我还没有对我的代码进行时间/内存分析。谢谢！ - Abhi

@Abhi，我不是算法方面的专家，但我尽可能地采取了每一个捷径。如果你用C语言实现它，你可以让它快上一两个命令。 - ikegami

@Mariya，这是基于“动态规划”解决最长公共子串问题的解决方案。在此处查看其基础解释：http://en.wikipedia.org/wiki/Longest_common_substring_problem#Dynamic_programming。 - ikegami