搜索允许在字符串的任何位置出现一个不匹配的字符串

Question

我正在处理长度为 25 的 DNA 序列（参见下面的示例）。我有一个 230,000 的列表，需要在整个基因组（弓形虫寄生虫）中寻找每个序列。我不确定基因组有多大，但比 230,000 个序列长得多。

我需要查找每个 25 个字符的序列，例如，( AGCCTCCCATGATTGAACAGATCAT )。

基因组被格式化为一个连续的字符串，即（ CATGGGAGGCTTGCGGAGCCTGAGGGCGGAGCCTGAGGTGGGAGGCTTGCGGAGTGCGGAGCCTGAGCCTGAGGGCGGAGCCTGAGGTGGGAGGCTT.... ）

我不在乎在哪里找到它或找到它多少次，只关心它是否存在。

我认为这很简单——

 str.find(AGCCTCCCATGATTGAACAGATCAT)

但是我也什么找到一个在任何位置定义为错误（不匹配）的紧密匹配，但只有一个位置，并记录序列中的位置。我不确定该怎么做。我唯一能想到的就是使用通配符并在每个位置使用通配符执行搜索。即搜索25次。

例如，

 AGCCTCCCATGATTGAACAGATCAT
AGCCTCCCATGATAGAACAGATCAT

位置 13 不匹配的近距离匹配。

速度不是大问题，因为我只做了 3 次，不过如果速度快点就好了。

有些程序可以执行此操作——查找匹配项和部分匹配项——但我正在寻找一种无法通过这些应用程序发现的部分匹配项。

这是 perl 的类似帖子，尽管它们只是比较序列而不是搜索连续的字符串：

相关帖子

原文由 Vincent 发布，翻译遵循 CC BY-SA 4.0 许可协议

Answer 1

在继续阅读之前，您 是否了解过 biopython ？

您似乎想要找到具有一个替换错误和零插入/删除错误（即汉明距离为 1）的近似匹配。

如果你有汉明距离匹配函数（参见例如 Ignacio 提供的链接），你可以像这样使用它来搜索第一个匹配项：

 any(Hamming_distance(genome[x:x+25], sequence) == 1 for x in xrange(len(genome)))

但这会相当慢，因为 (1) 汉明距离函数会在第二次替换错误后继续研磨 (2) 失败后，它将光标前进一个而不是根据它看到的内容向前跳过（就像 Boyer-摩尔搜索确实如此）。

您可以使用如下函数克服 (1)：

 def Hamming_check_0_or_1(genome, posn, sequence):
    errors = 0
    for i in xrange(25):
        if genome[posn+i] != sequence[i]:
            errors += 1
            if errors >= 2:
                return errors
    return errors

注意：这不是 Pythonic，而是 Cic，因为您需要使用 C（可能通过 Cython）来获得合理的速度。

Navarro 和 Raffinot（谷歌“Navarro Raffinot nrgrep”）已经完成了一些关于带跳过的位并行近似 Levenshtein 搜索的工作，这可以适用于 Hamming 搜索。请注意，位并行方法对查询字符串的长度和字母大小有限制，但你的分别是 25 和 4，所以没有问题。更新：对于字母大小为 4 的情况，跳过可能没有多大帮助。

当你在谷歌上搜索汉明距离搜索时，你会注意到很多关于在硬件中实现它的东西，而不是在软件中。这是一个很大的暗示，也许你想出的任何算法都应该用 C 或其他一些编译语言来实现。

更新： 位并行方法的工作代码

我还提供了一个简单的方法来帮助进行正确性检查，并且我打包了 Paul 的重新代码的变体以进行一些比较。请注意，使用 re.finditer() 会提供非重叠结果，这可能会导致距离为 1 的匹配掩盖精确匹配；请参阅我的最后一个测试用例。

位并行方法具有以下特点： 保证线性行为 O(N) 其中 N 是文本长度。注意天真的方法是 O(NM) 和正则表达式方法一样（M 是模式长度）。 Boyer-Moore 风格的方法是最坏情况 O(NM) 和预期 O(N)。当必须缓冲输入时，也可以轻松使用位并行方法：一次可以输入一个字节或一兆字节；没有前瞻性，缓冲区边界没有问题。最大的优势是：用 C 编写的简单的按输入字节代码的速度。

缺点：模式长度实际上受限于快速寄存器中的位数，例如 32 或 64。在这种情况下，模式长度为 25；没问题。它使用与模式中不同字符的数量成比例的额外内存（S_table）。在这种情况下，“字母大小”只有 4；没问题。

详细信息来自 这份技术报告。该算法用于使用 Levenshtein 距离进行近似搜索。为了转换为使用汉明距离，我简单地（！）删除了处理插入和删除的语句 2.1 的片段。您会注意到很多带有“d”上标的“R”参考。 “d”是距离。我们只需要 0 和 1。这些“R”成为下面代码中的 R0 和 R1 变量。

 # coding: ascii

from collections import defaultdict
import re

_DEBUG = 0

# "Fast Text Searching with Errors" by Sun Wu and Udi Manber
# TR 91-11, Dept of Computer Science, University of Arizona, June 1991.
# http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.20.8854

def WM_approx_Ham1_search(pattern, text):
    """Generate (Hamming_dist, start_offset)
    for matches with distance 0 or 1"""
    m = len(pattern)
    S_table = defaultdict(int)
    for i, c in enumerate(pattern):
        S_table[c] |= 1 << i
    R0 = 0
    R1 = 0
    mask = 1 << (m - 1)
    for j, c in enumerate(text):
        S = S_table[c]
        shR0 = (R0 << 1) | 1
        R0 = shR0 & S
        R1 = ((R1 << 1) | 1) & S | shR0
        if _DEBUG:
            print "j= %2d msk=%s S=%s R0=%s R1=%s" \
                % tuple([j] + map(bitstr, [mask, S, R0, R1]))
        if R0 & mask: # exact match
            yield 0, j - m + 1
        elif R1 & mask: # match with one substitution
            yield 1, j - m + 1

if _DEBUG:

    def bitstr(num, mlen=8):
       wstr = ""
       for i in xrange(mlen):
          if num & 1:
             wstr = "1" + wstr
          else:
             wstr = "0" + wstr
          num >>= 1
       return wstr

def Ham_dist(s1, s2):
    """Calculate Hamming distance between 2 sequences."""
    assert len(s1) == len(s2)
    return sum(c1 != c2 for c1, c2 in zip(s1, s2))

def long_check(pattern, text):
    """Naively and understandably generate (Hamming_dist, start_offset)
    for matches with distance 0 or 1"""
    m = len(pattern)
    for i in xrange(len(text) - m + 1):
        d = Ham_dist(pattern, text[i:i+m])
        if d < 2:
            yield d, i

def Paul_McGuire_regex(pattern, text):
    searchSeqREStr = (
        '('
        + pattern
        + ')|('
        + ')|('.join(
            pattern[:i]
            + "[ACTGN]".replace(c,'')
            + pattern[i+1:]
            for i,c in enumerate(pattern)
            )
        + ')'
        )
    searchSeqRE = re.compile(searchSeqREStr)
    for match in searchSeqRE.finditer(text):
        locn = match.start()
        dist = int(bool(match.lastindex - 1))
        yield dist, locn

if __name__ == "__main__":

    genome1 = "TTTACGTAAACTAAACTGTAA"
    #         01234567890123456789012345
    #                   1         2

    tests = [
        (genome1, "ACGT ATGT ACTA ATCG TTTT ATTA TTTA"),
        ("T" * 10, "TTTT"),
        ("ACGTCGTAAAA", "TCGT"), # partial match can shadow an exact match
        ]

    nfailed = 0
    for genome, patterns in tests:
        print "genome:", genome
        for pattern in patterns.split():
            print pattern
            a1 = list(WM_approx_Ham1_search(pattern, genome))
            a2 = list(long_check(pattern, genome))
            a3 = list(Paul_McGuire_regex(pattern, genome))
            print a1
            print a2
            print a3
            print a1 == a2, a2 == a3
            nfailed += (a1 != a2 or a2 != a3)
    print "***", nfailed

原文由 John Machin 发布，翻译遵循 CC BY-SA 2.5 许可协议

Answer 2

Python 正则表达式 库支持 模糊正则表达式匹配。与 TRE 相比的一个优势是它允许在文本中查找正则表达式的所有匹配项（也支持重叠匹配项）。

 import regex
m=regex.findall("AA", "CAG")
>>> []
m=regex.findall("(AA){e<=1}", "CAAG") # means allow up to 1 error
m
>>> ['CA', 'AG']

原文由 sefakilic 发布，翻译遵循 CC BY-SA 4.0 许可协议