python - 蛋白质互信息

Question

我正在尝试查找多序列比对 (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 出现的概率(或频率)。

Ps:关于MI更简单的解释请引用这个链接mistic.leloir.org.ar/docs/help.html “感谢亚伦”

Answer 1

我不是 100% 确定这是否正确(例如，应该如何处理 '-')？我假设求和是针对 log 中分子和分母中的频率均为非零的所有对，此外，我假设它应该是 natural 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)