2017-01-10 83 views
3

我有一個FASTA文件與DNA序列和序列的名稱,我需要做一個矩陣的重疊分數。我在Biopython中找到了模塊pairwise2,這看起來很好。除了我的序列已經對齊,並且當我使用pairwise2時,它再次嘗試對齊需要很長時間的序列,並且顯然對於每個對齊獲得相同的重疊分數。所以我的問題是如何在沒有嘗試重新排列序列的情況下獲得重疊分數? 這是我到目前爲止有:重疊分數矩陣biopython

from Bio.Alphabet import IUPAC 
from Bio import SeqIO 
from Bio import pairwise2 

fasta_file = SeqIO.parse('unambiguous.fasta', 'fasta', alphabet=IUPAC.ambiguous_dna) 

all_seq = [] 
for seq_record in fasta_file: 
    all_seq += [str(seq_record.seq)] 

compare = pairwise2.align.globalms(all_seq[0], all_seq[1], 2, -1, -1, 0) 
print(compare) 

我從FASTA文件只使用第一和第二序列這裏試訓。正如你在腳本中看到的,匹配應該獎勵2分,不匹配和差距-1。當兩個序列在​​同一個位置上有差距時,0應該是獎勵。我知道把0放在第4位不會給我想要的結果,但我還沒有解決這個問題的方法。此時對齊問題似乎更大。 因此,任何人都有一些與pairwise2或其他python/biopython模塊的經驗,可以讓我的重疊分數?

+0

你的意思是'unambiguous.fasta'包含對齊的序列嗎? –

+0

請[編輯]你的問題,包括示例你的問題的輸入。 – MattDMo

回答

0

就我所知,unambiguous.fasta含有比對的基因序列。您可以採用適合您需求的計分函數得分他們:

from itertools import starmap, combinations 


def score(seq1, seq2): 
    def score_(a, b): 
     return (0 if a == b == "-" # both are gaps 
       else -1 if a != b # mismatch or gap 
       else 2)   # match 

    return sum(starmap(score_, zip(seq1, seq2))) 

您可能要修改它忽略不明確的基地位置,因爲人們通常做。這裏是一個整潔的方式來比較所有序列:

sequences = SeqIO.parse('unambiguous.fasta', 'fasta', alphabet=IUPAC.ambiguous_dna) 
scores = starmap(score, combinations(sequences, 2)) 

一旦執行,scores(注意,這是一個懶惰迭代器)將生成分數的成對基質的扁平上三角。 score應該工作得很快,但如果有成千上萬的序列(即數百萬次計算比較),則可能需要重新實現Cython或Numba。

編輯在Python 2.x中,您可能需要用替換izip

+0

謝謝,它似乎是以一種簡單易懂的方式來實現的。它只是不想打印出分數,因爲我得到: 任何想法可能意味着什麼? – JDh

+0

@JDh不客氣。在我寫的答案中:「注意它是一個懶惰的迭代器」。您可以迭代它並逐個執行某些分數(如果您不想將所有分數一次加載到RAM中,您可以這樣做)或將其轉換爲非惰性容器,例如, '名單(分數)'。 Python大量使用懶惰評估。如果你想有效地使用語言,你應該讓自己對這個概念感到滿意。 –