Try/Except 和循环的 DNA 基序枚举 - Python3

Question

过去几天我一直在努力解决一个特定的生物信息学问题，我想知道是否有人可以找到我的逻辑中的错误或我的代码中的错误（或两者）。该函数应该从所有 DNA 串中找到汉明距离最多为 d 的所有 k 聚体。这就是我想要做的：

从所有可能的 k-mers 的迭代开始，并将它们与每个字符串进行比较
这意味着我还需要另一个遍历 DNA 串的循环
我为 c+k <= len(DNA[0])-1 创建了一个 while 循环。 c+k 是我的 window，大小为 k，我想在每个 DNA 串中找到至少一个 window，其中我的组合与该串的汉明距离等于或小于任意 d。如果汉明距离满足标准，则 while 循环中断，允许比较下一个字符串。如果没有，则更改window，如果c+k==len(DNA[0])-1，并且汉明距离仍然不符合标准，我创建一个名称错误int(a)并异常终止inner_loop.

但是，我的函数 returns 除了 set() 什么都没有，我不明白。

import itertools
def combination(k):
    bases=['A','T','G','C']
    combo=[''.join(p) for p in itertools.product(bases, repeat=k)]
    return combo

def hammingDistance(Pattern, seq):
        if Pattern == seq:
               return 0
        else:
                dist=0
                for i in range(len(seq)):
                        if Pattern[i] != seq[i]:
                                dist += 1
        return dist

def motif_enumeration(k, d, DNA):
    combos = combination(k)
    global pattern
    for combo in combos:
        try:
            inner_loop(k, d, DNA, combo)
        except:
            continue

    return set(pattern)

def inner_loop(k, d, DNA, combo):
    global pattern
    for strings in DNA:
        inner_loop_two(k, d, DNA, combo, strings)


def inner_loop_two(k, d, DNA, combo, strings):
    global pattern
    c=0
    while c+k < len(DNA[0]):
        print(combo, strings[c:c+k], hammingDistance(combo, strings[c:c+k]))
        if d >= hammingDistance(combo, strings[c:c+k]) and strings == DNA[len(DNA)-1]:
            #if we've reached the last string and the condition is met,
            #that means that the combo is suitable for each string of DNA
            pattern += [combo]
        elif d >= hammingDistance(combo, strings[c:c+k]):
            #condition is met for one string, now move onto next
            break
        elif d < hammingDistance(combo, strings[c:c+k]) and c+k == len(DNA[0])-1:
            #Name error causes this inner loop two to crash, thus causing the first inner loop
            #to pass
            int(a)
        elif d < hammingDistance(combo, strings[c:c+k]):
            #change the window to see if the combo is valid later in the string
            c += 1


pattern = []
DNA=['ATTTGGC', 'TGCCTTA', 'CGGTATC', 'GAAAATT']
print(motif_enumeration(3,1,DNA))
print(pattern)

我认为由于我的第二个内循环崩溃，这将导致我的第一个内循环通过，然后 motif_enumeration 中的另一个组合将被测试，但我的 inner_loop_two 中的第一个条件从不打印任何东西。我还注意到，当内循环崩溃并且 motif_enumeration 继续时，外循环和内循环都会继续。这是我的意思的一个例子...

AAA ATT 2
AAA TTT 3
AAA TTG 3
AAA TGG 3
AAT ATT 1
AAT TGC 3
AAT GCC 3
AAT CCT 2
AAT CTT 2
AAG ATT 2
AAG TTT 3
AAG TTG 2
AAG TGG 2
AAC ATT 2
AAC TTT 3
AAC TTG 3
AAC TGG 3
ATA ATT 1 etc...

我的预期输出是 pattern=[ATA, ATT, GTT, TTT]

Answer 1

逻辑的核心组成部分是，如果组合在 all[=36] 上的 any 位置匹配，我们希望将组合收集到模式集中=] 的目标字符串。我们可以使用 Python 的 all 和 any 函数来做到这一点。这些函数工作效率很高，因为它们会在结果确定后立即停止测试。

import itertools def combination(k): return (''.join(p) for p in itertools.product('ATCG', repeat=k)) def hamming_distance(pattern, seq): return sum(c1 != c2 for c1, c2 in zip(pattern, seq)) def window(s, k): for i in range(1 + len(s) - k): yield s[i:i+k] def motif_enumeration(k, d, DNA): pattern = set() for combo in combination(k): if all(any(hamming_distance(combo, pat) <= d for pat in window(string, k)) for string in DNA): pattern.add(combo) return pattern DNA = ['ATTTGGC', 'TGCCTTA', 'CGGTATC', 'GAAAATT'] print(motif_enumeration(3, 1, DNA))

输出

{'GTT', 'ATA', 'TTT', 'ATT'}

我对您的代码进行了一些其他更改。我们可以通过将生成器传递给 sum 函数来有效地计算汉明距离。我们可以通过使用生成器将组合元组转换为字符串而不是将它们放入列表中来节省时间和 RAM。

motif_enumeration 函数可以进一步压缩成集合理解，但我必须承认它相当密集，甚至比以前的版本更难阅读。不过，它的效率可能会稍微高一些。

def motif_enumeration(k, d, DNA): return {combo for combo in combination(k) if all(any(hamming_distance(combo, pat) <= d for pat in window(string, k)) for string in DNA)}

这里是一个稍微更易读的版本，我在其中提供了 motif_enumeration 辅助函数 in_window 来执行内部测试。

# Return True if combo is within d in any window of string def in_window(combo, string, k, d): return any(hamming_distance(combo, pat) <= d for pat in window(string, k)) def motif_enumeration(k, d, DNA): pattern = set() for combo in combination(k): if all(in_window(combo, string, k, d) for string in DNA): pattern.add(combo) return pattern

Try/Except 和循环的 DNA 基序枚举 - Python3

DNA Motif Enumeration with Try/Except and Loops - Python3

python

loops

try-catch

bioinformatics

hamming-distance