python - 蛋白质相互信息

我正在尝试查找多序列比对（MSA）之间的相互信息（MI）。

对于我来说，它背后的数学是可以的。虽然，但我不知道如何至少以一种快速的方式在Python中实现它。

我应该如何计算总频率P(i;x); P(j;y); P(ij;xy)。 Px和Py频率很容易计算哈希可以处理，但是P(ij;xy)呢？

因此，我真正的问题是，如何计算给定的i和j列中Pxy的概率？

请注意，MI可以定义为：

MI(i,j) = Sum(x->n)Sum(y->m) P(ij,xy) * log(P(ij,xy)/P(i,x)*P(j,y))

其中i和j是列中的氨基酸位置，x和y是在给定的i或j列中找到的不同氨基酸。

谢谢，

编辑

我的输入数据看起来像一个df：

A = [
['M','T','S','K','L','G','-'.'-','S','L','K','P'],
['M','A','A','S','L','A','-','A','S','L','P','E'],
...,
['M','T','S','K','L','G','A','A','S','L','P','E'],
]

因此，确实很容易计算给定位置的任何氨基酸频率，
 例如：

P(M) at position 1: 1
P(T) at position 2: 2/3
P(A) at position 2: 1/3
P(S) at position 3: 2/3
P(A) at position 3: 1/3

例如，我应该如何继续获取位置2处的T和位置3处的S的P：
在此示例中为2/3。

因此，P（ij，xy）表示第i列中的氨基酸x与第j列中的氨基酸y同时出现的概率（或频率）。

附言：有关MI的更简单说明，请参考此链接mistic.leloir.org.ar/docs/help.html“感谢亚伦”

最佳答案

我不确定100％是否正确（例如，应该如何处理'-'）？我假设总和在log内的分子和分母中的频率均为非零的所有对中，此外，我假定它应该是自然对数：

from math import log
from collections import Counter

def MI(sequences,i,j):
    Pi = Counter(sequence[i] for sequence in sequences)
    Pj = Counter(sequence[j] for sequence in sequences)
    Pij = Counter((sequence[i],sequence[j]) for sequence in sequences)

    return sum(Pij[(x,y)]*log(Pij[(x,y)]/(Pi[x]*Pj[y])) for x,y in Pij)

该代码通过使用3个Counter对象获取相关计数，然后返回总和（该总和是公式的直接转换）来工作。

如果这是不正确的，那么如果您编辑问题以使它具有一些预期的输出可进行测试，则将很有帮助。

在编辑上。这是一个版本，它不将'-'视为另一个氨基酸，而是过滤掉出现在两列中任一列中的序列，将这些序列解释为无法获得必需信息的序列：

def MI(sequences,i,j):
    sequences = [s for s in sequences if not '-' in [s[i],s[j]]]
    Pi = Counter(s[i] for s in sequences)
    Pj = Counter(s[j] for s in sequences)
    Pij = Counter((s[i],s[j]) for s in sequences)

    return sum(Pij[(x,y)]*log(Pij[(x,y)]/(Pi[x]*Pj[y])) for x,y in Pij)

关于python - 蛋白质相互信息，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/42772986/