我正在尝试查找多序列比对 (MSA) 之间的互信息 (MI)。
它背后的数学对我来说没问题。不过,我不知道如何用 Python 实现它,至少不知道如何快速实现。
我应该如何计算总体频率P(i;x)
; P(j;y)
; P(ij;xy)
。 Px
和 Py
频率很容易计算,哈希可以处理它,但是 P(ij;xy)
呢? p>
所以我真正的问题是,我应该如何计算给定 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 的 P 和位置 3 处 S 的 P: 在此示例中为 2/3。
因此 P(ij, xy) 表示 i 列中的氨基酸 x 与 j 列中的氨基酸 y 同时出现的概率(或频率)。
Ps:关于MI更简单的解释请引用这个链接mistic.leloir.org.ar/docs/help.html “感谢亚伦”
最佳答案
我不能 100% 确定这是否正确(例如,'-'
应该如何处理)?我假设总和是对数内分子和分母的频率均非零的所有对的总和,此外,我假设它应该是自然对数:
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/