重叠分数矩阵biopython

overlap score matrix biopython

我有一个包含 DNA 序列和序列名称的 FASTA 文件,我需要制作一个重叠分数矩阵。我在 Biopython 中找到了 pairwise2 模块,它似乎做得很好。除了我的序列已经对齐并且当我使用 pairwise2 它再次尝试对齐需要很长时间的序列并且显然每次对齐都获得相同的重叠分数。所以我的问题是如何在不尝试再次对齐序列的情况下获得重叠分数? 这是我目前所拥有的:

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 模块有一定经验的人都可以得到重叠分数吗?

据我了解,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 上,您可能希望将 zip 替换为 izip