2016-11-30 47 views
0

你好,我很新的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 ...

+0

你只是循環一次。您需要遍歷序列,然後在遇到要更改的字母(即不是A,C,G或T)時,循環遍歷可能的替換。你將不得不嵌套這些循環來獲得所有的排列。 – samb8s

+0

只是一個提示;你可以做'如果我在'RWMDHVN':poss_seq + ='A''而不是寫作多個if-elif的語句 –

+0

是的,那就是'我想說的最後一句話。我明白,但我不知道如何實施它? –

回答

7

您可以使用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返回一個生成器,而不是一次構建所有可能的序列,因此您可以隨時停止迭代,而無需等待每個要生成的序列。

+0

我還沒有玩過itertools,但它看起來非常強大。這似乎是做到這一點的最佳方式 – samb8s

+0

你打敗了我。 :)你想讓我編輯我的完整版「MAP」到你的答案嗎? FWIW,我可能只是做'MAP [c]'而不是'MAP.get(c,'')',因爲如果我們收到錯誤的數據,我們可能想要捕捉關鍵錯誤。 –

+0

@ PM2Ring謝謝你的提議,但我可以做到。使用'get'的唯一原因是第二個例子中有'-',但結果是缺失。如果需要捕捉錯誤,那麼索引操作符就是要走的路。 – niemmi