定义一个返回不同可能性的循环

Question

你好，我对python很陌生。我有以下问题： 我想写一个脚本，给定一个有歧义的（dna）序列，写下所有可能的序列（如果少于100个，如果超过100个可能的序列，则打印一条适当的错误消息） 对于DNA核苷酸的模糊性：http://www.bioinformatics.org/sms/iupac.html

示例：对于序列“AYGH”，脚本的输出将是“ACGA”, “ACGC”, “ACGT”, “ATGA”, “ATGC”和“ATGT”。A、 C，G和T是默认的核苷酸。所有其他的可以有不同的值（参见链接）。你知道吗

所以我写了这个：

def possible_sequences (seq):
    poss_seq = ''
    for i in seq:
        if i=='A'or i=='C'or i=='G'or i=='T': 
            poss_seq += i 
        else: 
            if i== 'R':  
                poss_seq += 'A' # OR 'G', how should i implement this? 
            elif i == 'Y': 
                poss_seq += 'C' # OR T 
            elif i == 'S': 
                poss_seq += 'G' # OR C
            elif i == 'W': 
                poss_seq += 'A' # OR T 
            elif i == 'K': 
                poss_seq += 'G' # OR T
            elif i == 'M': 
                poss_seq += 'A' # OR C
            elif i == 'B': 
                poss_seq += 'C' # OR G OR T 
            elif i == 'D': 
                poss_seq += 'A' # OR G OR T 
            elif i == 'H': 
                poss_seq += 'A' # OR C OR T 
            elif i == 'V': 
                poss_seq += 'A' # OR C OR G 
            elif i == 'N': 
                poss_seq += 'A' # OR C OR G OR T 
            elif i == '-' or i == '.': 
                poss_seq += ' '
    return poss_seq

当我测试我的功能时： 可能的\u序列（'ATRY-C'） 我得到了：

'ATAC C'

但我应该得到：

'ATAC C'
'ATAT C' 
'ATGC C'
'ATGT C'

有人能帮帮我吗？我明白，当出现歧义时，我必须重述并写第二个poss\ seq，但我不知道如何。。。你知道吗

Answer 1

您可以使用^{}生成可能性：

from itertools import product

# List possible nucleotides for each possible item in sequence
MAP = {
    'A': 'A',
    'C': 'C',
    'G': 'G',
    'T': 'T',
    'R': 'AG',
    'Y': 'CT',
    'S': 'GC',
    'W': 'AT',
    'K': 'GT',
    'M': 'AC',
    'B': 'CGT',
    'D': 'AGT',
    'H': 'ACT',
    'V': 'ACG',
    'N': 'ACGT',
    '-': ' ',
    '.': ' '
}

def possible_sequences(seq):
    return (''.join(c) for c in product(*(MAP[c] for c in seq)))

print(list(possible_sequences('AYGH')))
print(list(possible_sequences('ATRY-C')))

输出：

['ACGA', 'ACGC', 'ACGT', 'ATGA', 'ATGC', 'ATGT']
['ATAC C', 'ATAT C', 'ATGC C', 'ATGT C']

在上面，我们首先迭代给定序列中的项目，得到每个项目可能的核苷酸列表：

possibilities = [MAP[c] for c in 'ATRY-C']
print(possibilities)

# ['A', 'T', 'AG', 'CT', ' ', 'C']

然后将iterable解压为给product的参数，该参数将返回笛卡尔积：

products = list(product(*['A', 'T', 'AG', 'CT', ' ', 'C']))
print(products)

# [('A', 'T', 'A', 'C', ' ', 'C'), ('A', 'T', 'A', 'T', ' ', 'C'), 
#  ('A', 'T', 'G', 'C', ' ', 'C'), ('A', 'T', 'G', 'T', ' ', 'C')]

最后，每个产品都变成一个带有^{}的字符串：

print(list(''.join(p) for p in products))

# ['ATAC C', 'ATAT C', 'ATGC C', 'ATGT C']

请注意，possible_sequences返回一个生成器，而不是一次构造所有可能的序列，因此您可以随时轻松停止迭代，而不必等待生成每个序列。你知道吗