代码之家 › 专栏 › 技术社区 › DrJessop

使用Try/Except和循环进行DNA基序计数-Python3

hamming-distance bioinformatics try-catch loops python

DrJessop · 技术社区 · 7 年前

该函数应该从所有DNA字符串中找到汉明距离最多为d的所有k-mer。

从所有可能的k-mer的迭代开始,并将其与每个字符串进行比较
我做了一个while循环 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 从不打印任何内容。我还注意到,当内部循环崩溃时继续,它对外部和内部循环都继续。这是一个我的意思的例子。。。

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]

1 回复 | 直到 7 年前

PM 2Ring 7 年前

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

这个 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\u枚举 辅助函数 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