python - 寻找不匹配模式的效率

Question

我正在研究一个简单的生物信息学问题。我有一个可行的解决方案，但效率低得离谱。我怎样才能提高我的效率？

---

问题:

在字符串 g 中查找长度为 k 的模式，给定 k-mer 最多可以有 d 个不匹配。

这些字符串和模式都是基因组的——所以我们的可能字符集是{A, T, C, G}。

我将调用函数 FrequentWordsMismatch(g, k, d)。

所以，这里有一些有用的例子:

FrequentWordsMismatch('AAAAAAAAAA', 2, 1) → ['AA', 'CA', 'GA', 'TA', 'AC', 'AG', '在']

这里有一个更长的例子，如果你实现它并想测试:

FrequentWordsMisMatch('CACAGTAGGCGCCGGCACACACAGCCCCGGGCCCCGGGCCGCCCCGGGCCGGCGGCCGCCGGCGCCGGCACACCGGCACAGCCGTACCGGCACAGTAGTACCGGCCGGCCGGCACACCGGCACACCGGGTACACACCGGGGCGCACACACAGGCGGGCGCCGGGCCCCGGGCCGTACCGGGCCGCCGGCGGCCCACAGGCGCCGGCACAGTACCGGCACACACAGTAGCCCACACACAGGCGGGCGGTAGCCGGCGCACACACACACAGTAGGCGCACAGCCGCCCACACACACCGGCCGGCCGGCACAGGCGGGCGGGCGCACACACACCGGCACAGTAGTAGGCGGCCGGCGCACAGCC', 10, 2) → ['GCACACAGAC', 'GCGCACACAC']

使用我天真的解决方案，第二个示例很容易花费大约 60 秒，但第一个示例非常快。

---

简单的解决方案:

我的想法是，对于 g 中的每个 k 长度段，找到每个可能的“邻居”(例如，其他 k 长度段，最多 d 不匹配)并将这些邻居作为键添加到字典中。然后我计算这些邻居 kmers 中的每一个出现在字符串 g 中的次数，并将它们记录在字典中。

显然，这样做有点糟糕，因为邻居的数量随着 k 和 d 的增加而疯狂地增加，并且必须扫描每个字符串这些邻居使得这个实现非常缓慢。但是，唉，这就是我寻求帮助的原因。

我会把我的代码放在下面。肯定有很多新手错误需要解包，所以感谢您的时间和关注。

def FrequentWordsMismatch(g, k, d):
    '''
    Finds the most frequent k-mer patterns in the string g, given that those 
    patterns can mismatch amongst themselves up to d times

    g (String): Collection of {A, T, C, G} characters
    k (int): Length of desired pattern
    d (int): Number of allowed mismatches
    '''
    counts = {}
    answer = []

    for i in range(len(g) - k + 1):
        kmer = g[i:i+k]
        for neighborkmer in Neighbors(kmer, d):
            counts[neighborkmer] = Count(neighborkmer, g, d)

    maxVal = max(counts.values())

    for key in counts.keys():
        if counts[key] == maxVal:
            answer.append(key)

    return(answer)


def Neighbors(pattern, d):
    '''
    Find all strings with at most d mismatches to the given pattern

    pattern (String): Original pattern of characters
    d (int): Number of allowed mismatches
    '''
    if d == 0:
        return [pattern]

    if len(pattern) == 1:
        return ['A', 'C', 'G', 'T']

    answer = []

    suffixNeighbors = Neighbors(pattern[1:], d)

    for text in suffixNeighbors:
        if HammingDistance(pattern[1:], text) < d:
            for n in ['A', 'C', 'G', 'T']:
                answer.append(n + text)
        else:
            answer.append(pattern[0] + text)

    return(answer)


def HammingDistance(p, q):
    '''
    Find the hamming distance between two strings

    p (String): String to be compared to q
    q (String): String to be compared to p
    '''
    ham = 0 + abs(len(p)-len(q))

    for i in range(min(len(p), len(q))):
        if p[i] != q[i]:
            ham += 1

    return(ham)


def Count(pattern, g, d):
    '''
    Count the number of times that the pattern occurs in the string g, 
    allowing for up to d mismatches

    pattern (String): Pattern of characters
    g (String): String in which we're looking for pattern
    d (int): Number of allowed mismatches
    '''
    return len(MatchWithMismatch(pattern, g, d))

def MatchWithMismatch(pattern, g, d):
    '''
    Find the indicies at which the pattern occurs in the string g, 
    allowing for up to d mismatches

    pattern (String): Pattern of characters
    g (String): String in which we're looking for pattern
    d (int): Number of allowed mismatches
    '''
    answer = []
    for i in range(len(g) - len(pattern) + 1):
        if(HammingDistance(g[i:i+len(pattern)], pattern) <= d):
            answer.append(i)
    return(answer)

---

更多测试

FrequentWordsMismatch('ACGTTGCATGTCGCATGATGCATGAGAGCT', 4, 1) → ['ATGC', 'ATGT', 'GATG']

FrequentWordsMismatch('AGTCAGTC', 4, 2) → ['TCTC', 'CGGC', 'AAGC', 'TGTG', 'GGCC', 'AGGT', 'ATCC', 'ACTG', 'ACAC', 'AGAG', 'ATTA', 'TGAC', 'AATT', 'CGTT', 'GTTC', 'GGTA', 'AGCA', 'CATC']

FrequentWordsMismatch('AATTAATTGGTAGGTAGGTA', 4, 0) → ["GGTA"]

FrequentWordsMismatch('ATA', 3, 1) → ['GTA', 'ACA', 'AAA', 'ATC', 'ATA', 'AGA', 'ATT', 'CTA', 'TTA', 'ATG']

FrequentWordsMismatch('AAT', 3, 0) → ['AAT']

FrequentWordsMismatch('TAGCG', 2, 1)  → ['GG', 'TG']

Answer 1

问题描述在几个方面是模棱两可的，所以我将通过示例进行说明。你似乎想要所有k -来自字母表的长度字符串 (A, C, G, T}这样与 g 的连续子串的匹配数是最大的 - 其中“匹配”表示逐个字符最多等于 d字符不等式。

我忽略了你的 HammingDistance()即使输入具有不同的长度，函数也会产生一些东西，主要是因为它对我来说没有多大意义;-)，但部分是因为在你给出的任何例子中不需要它来获得你想要的结果。

以下代码在所有示例中生成您想要的结果，即生成您提供的输出列表的排列。如果您想要规范的输出，我建议在返回之前对输出列表进行排序。

该算法非常简单，但依赖于 itertools “以 C 速度”进行重型组合提升。所有示例的运行速度都低于第二个总数。

对于每个长度- k g 的连续子串, 考虑所有 combinations(k, d)套d不同的索引位置。有4**d用 {A, C, G, T} 中的字母填充这些索引位置的方法，并且每一种这样的方式都是“一个模式”，最多匹配子字符串 d差异。通过记住已经生成的模式来剔除重复项；这比一开始就只生成独特的模式做出巨大的努力要快。

所以，总的来说，时间要求是O(len(g) * k**d * 4**d) = O(len(g) * (4*k)**d , 其中k**d是，对于 k 的相当小的值和 d ，二项式系数的夸大代表 combinations(k, d) .需要注意的重要一点是 - 毫不奇怪 - 它在 d 中呈指数增长。 .

def fwm(g, k, d):
    from itertools import product, combinations
    from collections import defaultdict

    all_subs = list(product("ACGT", repeat=d))
    all_ixs = list(combinations(range(k), d))
    patcount = defaultdict(int)

    for starti in range(len(g)):
        base = g[starti : starti + k]
        if len(base) < k:
            break
        patcount[base] += 1
        seen = set([base])
        basea = list(base)
        for ixs in all_ixs:
            saved = [basea[i] for i in ixs]
            for newchars in all_subs:
                for i, newchar in zip(ixs, newchars):
                    basea[i] = newchar
                candidate = "".join(basea)
                if candidate not in seen:
                    seen.add(candidate)
                    patcount[candidate] += 1
            for i, ch in zip(ixs, saved):
                basea[i] = ch

    maxcount = max(patcount.values())
    return [p for p, c in patcount.items() if c == maxcount]

编辑:唯一地生成模式

与其保留一组到目前为止看到的重复项来剔除重复项，不如直接从一开始就防止生成重复项。事实上，下面的代码更短更简单，尽管有些微妙。作为减少冗余工作的返回，有对 inner() 的递归调用层。功能。哪种方式更快似乎取决于具体的输入。

def fwm(g, k, d):
    from collections import defaultdict

    patcount = defaultdict(int)
    alphabet = "ACGT"
    allbut = {ch: tuple(c for c in alphabet if c != ch)
              for ch in alphabet}

    def inner(i, rd):
        if not rd or i == k:
            patcount["".join(base)] += 1
            return
        inner(i+1, rd)
        orig = base[i]
        for base[i] in allbut[orig]:
            inner(i+1, rd-1)
        base[i] = orig

    for i in range(len(g) - k + 1):
        base = list(g[i : i + k])
        inner(0, d)

    maxcount = max(patcount.values())
    return [p for p, c in patcount.items() if c == maxcount]