如何从几种蛋白质序列中提取多态位置？

Question

我有几个fasta蛋白比对(~5000)，我想要识别多态位置和氨基酸残基之间的变化。我尝试过自己编写代码，但这是非常困难的(我在编程方面是新手)，我看过BioPython，但还没有找到任何东西。我想要这样的东西：

蛋白比对：

> sp1 MQGAAYMQAAAYYMQA
> sp2 MQGAARMQGAAYYMQA
> sp3 MQGAARMQGAAYYMQM
> sp4 MQGAARMQGAAYYMQA
> sp5 MQGAARMQAAAYYMQA
           ^  ^      ^

在上面的示例中，对齐有3个多态位置(标记为^)。第一位位于第六位，第二位位于第九位，第三位位于第十六位。多态位点coult的一个常见表示法是: R6Y，这意味着从R到Y的第6位发生了变化。变化的方向(R->Y或Y->R)是基于该位置上最频繁的字母。因此，在这种情况下，R有最高的频率，人们可以推断方向是R->Y。

正如你所看到的，第6和第16位有单一的变化(不同的字母频率为1)。然而，第9位有两个序列(sp1和sp5)随着这一变化。我想区分这两种类型的多态性。因此，在本例中，我希望输出如下所示：

输出:

# Alignment #1
#   Single polymorphisms:
#     R6Y: sp1
#     A16M: sp3

#   Non-single polymorphisms:
#     G9A: sp1, sp5

我希望这有助于澄清(如果时间太长了，对不起)。

任何建议都是非常感谢的，谢谢！

Answer 1

这里有一个函数，它将发现两个等位基因之间的差异。第一个应该是由每个位置上最常见的字母组成的标准多角形(在您的例子中是MQGAARMQGAAYYMQA)。

def polymorphic_positions(allele1, allele2):
    return [p[1][0] + str(p[0] + 1) + p[1][1] 
            for p in enumerate(zip(allele1, allele2))
            if p[1][0] != p[1][1]]

示例：

>>> polymorphic_positions('MQGAARMQGAAYYMQA', 'MQGAARMQGAAYYMQM')
['A16M']
>>> polymorphic_positions('MQGAARMQGAAYYMQA', 'MQGAARMQAAAYYMQA')
['G9A']

以下是一些参考资料：

• 清单理解
• zip
• 枚举