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

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 

当我测试功能： possible_sequences（ 'Â尝试-C'） 我得到：

'ATAC C' 

但是我应该得到：

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

有人可以帮我吗？我明白，我要回顾一下，当有歧义存在，但我不知道怎么写了第二poss_seq ...

Answer 1

您可以使用itertools.product生成的可能性：

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'] 

然后迭代被取出作为给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')] 

最后，产品的每一个被关到一个字符串join：

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

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

注意possible_sequences返回一个生成器，而不是一次构建所有可能的序列，因此您可以随时停止迭代，而无需等待每个要生成的序列。