2017-01-10 85 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大量使用懒惰评估。如果你想有效地使用语言,你应该让自己对这个概念感到满意。 –