DNA 基序枚举与 try/except 和循环 - Python3.

Question

在过去的几天里，我一直在为一个特定的生物信息学问题而苦苦挣扎，我想知道是否有人可以找到我的逻辑中的错误或我的代码中的错误(或两者兼而有之(。 该函数应该找到所有汉明距离最多为d的k-mers与所有DNA链。 这就是我要做的：

1. 从所有可能的 k-mers 的迭代开始，并将它们与每个字符串进行比较

2. 这意味着我需要另一个穿过DNA链的环路。

3. 我为c+k <= len(DNA[0])-1做了一段时间循环.c+k是我大小k的窗口，我想在每个DNA链中找到至少一个窗口，其中我的组合与该串的汉明距离等于或小于任意d。如果汉明距离满足标准，则 while 循环中断，允许比较下一个字符串。如果没有，则更改窗口，如果c+k==len(DNA[0])-1，并且汉明距离仍然不符合标准，我创建一个名称错误int(a)并且异常会杀死inner_loop。

但是，我的函数除了我不明白的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

逻辑的核心组件是，如果组合在所有目标字符串上的任何位置匹配，我们希望将组合收集到模式集中。我们可以使用 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