如何在Python中计算加权邻接矩阵的拓扑重叠度测量[TOM]?

Question

我正在尝试计算带权拓扑重叠（weighted topological overlap）的邻接矩阵，但是我无法使用 numpy 正确完成它。正确实现的 R 函数来自于 WGCNA (https://www.rdocumentation.org/packages/WGCNA/versions/1.67/topics/TOMsimilarity)。计算这个的公式（我认为）在第4个方程式中详细说明，我相信以下是正确的复制。

有谁知道如何正确实现它以反映 WGCNA 版本吗？

是的，我知道 rpy2，但是如果可能的话，我想要使用轻量级的方法。

首先，我的对角线不是 1，并且值与原始值没有一致的误差（例如，不是所有的值都偏离了 x）。

当我在 R 中计算这个时，我使用以下内容：

> library(WGCNA, quiet=TRUE)
> df_adj = read.csv("https://pastebin.com/raw/sbAZQsE6", row.names=1, header=TRUE, check.names=FALSE, sep="\t")
> df_tom = TOMsimilarity(as.matrix(df_adj), TOMType="unsigned", TOMDenom="min")
# ..connectivity..
# ..matrix multiplication (system BLAS)..
# ..normalization..
# ..done.
# I've uploaded it to this url: https://pastebin.com/raw/HT2gBaZC

我不确定我的代码哪里有误。 R 版本的源代码在这里，但它使用了 C 后端脚本，这对我来说非常难以理解。

以下是我在 Python 中的实现：

import pandas as pd
import numpy as np
from sklearn.datasets import load_iris

def get_iris_data():
    iris = load_iris()
    # Iris dataset
    X = pd.DataFrame(iris.data,
                     index = [*map(lambda x:f"iris_{x}", range(150))],
                     columns = [*map(lambda x: x.split(" (cm)")[0].replace(" ","_"), iris.feature_names)])

    y = pd.Series(iris.target,
                           index = X.index,
                           name = "Species")
    return X, y

# Get data
X, y = get_iris_data()

# Create an adjacency network
# df_adj = np.abs(X.T.corr()) # I've uploaded this part to this url: https://pastebin.com/raw/sbAZQsE6
df_adj = pd.read_csv("https://pastebin.com/raw/sbAZQsE6", sep="\t", index_col=0)
A_adj = df_adj.values

# Correct TOM from WGCNA for the A_adj
# See above for code
# https://www.rdocumentation.org/packages/WGCNA/versions/1.67/topics/TOMsimilarity
df_tom__wgcna = pd.read_csv("https://pastebin.com/raw/HT2gBaZC", sep="\t", index_col=0)

# My attempt
A = A_adj.copy()
dimensions = A.shape
assert dimensions[0] == dimensions[1]
d = dimensions[0]

# np.fill_diagonal(A, 0)

# Equation (4) from http://dibernardo.tigem.it/files/papers/2008/zhangbin-statappsgeneticsmolbio.pdf
A_tom = np.zeros_like(A)
for i in range(d):
    a_iu = A[i]
    k_i = a_iu.sum()
    for j in range(i+1, d):
        a_ju = A[:,j]
        k_j = a_ju.sum()
        l_ij = np.dot(a_iu, a_ju)
        a_ij = A[i,j]
        numerator = l_ij + a_ij
        denominator = min(k_i, k_j) + 1 - a_ij
        w_ij = numerator/denominator
        A_tom[i,j] = w_ij
A_tom = (A_tom + A_tom.T)

有一个名为GTOM的包（https://github.com/benmaier/gtom），但它不适用于加权邻接矩阵。 GTOM的作者也研究了这个问题（使用了更为复杂/高效的NumPy实现，但仍未产生预期结果）。 有人知道如何复制WGCNA实现吗？ 编辑：2019年06月20日 我从@scleronomic和@benmaier那里借鉴了一些代码，并在doc string中列出了他们的贡献。该函数可在soothsayer中使用，版本为v2016.06及以上。希望这将使Python中的拓扑重叠比更容易地使用，而不仅限于R。 https://github.com/jolespin/soothsayer/blob/master/soothsayer/networks/networks.py

import numpy as np
import soothsayer as sy
df_adj = sy.io.read_dataframe("https://pastebin.com/raw/sbAZQsE6")
df_tom = sy.networks.topological_overlap_measure(df_adj)
df_tom__wgcna = sy.io.read_dataframe("https://pastebin.com/raw/HT2gBaZC")
np.allclose(df_tom, df_tom__wgcna)
# True

- O.rka

Answer 1

首先，让我们看一下针对二进制邻接矩阵a_ij的方程式的各个部分：

• a_ij: 表示节点i是否与节点j相连
• k_i: 节点i的邻居数量（连通性）
• l_ij: 节点i和节点j的公共邻居数量

因此，w_ij衡量的是较低连接性节点的多少邻居也是另一个节点的邻居（即w_ij衡量"它们的相对互连性"）。

我猜他们将A的对角线定义为零而不是一。 在这种假设下，我可以重现WGCNA的值。

A[range(d), range(d)] = 0  # Assumption
L = A @ A  # Could be done smarter by using the symmetry
K = A.sum(axis=1)

A_tom = np.zeros_like(A)
for i in range(d):
    for j in range(i+1, d):  
        numerator = L[i, j] + A[i, j]
        denominator = min(K[i], K[j]) + 1 - A[i, j]
        A_tom[i, j] = numerator / denominator
    
A_tom += A_tom.T
A_tom[range(d), range(d)] = 1  # Set diagonal to 1 by default

A_tom__wgcna = np.array(pd.read_csv("https://pastebin.com/raw/HT2gBaZC", 
                        sep="\t", index_col=0))
print(np.allclose(A_tom, A_tom__wgcna))

A的对角线为零而不是一的直觉可以通过对具有二进制A的简单示例进行观察得出：

 Graph      Case Zero    Case One
   B          A B C D      A B C D  
 /   \      A 0 1 1 1    A 1 1 1 1  
A-----D     B 1 0 0 1    B 1 1 0 1  
 \   /      C 1 0 0 1    C 1 0 1 1  
   C        D 1 1 1 0    D 1 1 1 1  

给定方程4的描述解释如下：

注意，如果连接较少的节点满足两个条件，则w_ij = 1：

（a）它的所有邻居也是另一个节点的邻居；

（b）它与另一个节点相连。

相反，如果和j未连接且两个节点没有共享任何邻居，则w_ij = 0。

因此，A-D之间的连接应满足这个标准，即w_14=1。

情况零对角线： 情况一对角线： 当应用公式时还缺失了什么呢？对角线值不匹配。我将默认设置为1。节点自身的互联性是什么意思？与1（或0，根据定义）不同的值对我来说没有意义。在简单的例子中，无论是零情况还是一情况，都没有得到w_ii=1 的结果。在零情况中，必须有k_i+1 == l_ii，在一情况中，必须有k_i == l_ii+1，这两种情况都似乎是错的。

所以总结一下，我会将邻接矩阵的对角线设置为0，使用给定的方程，并将结果的对角线默认设置为1。

- scleronomic

---

感谢您帮助我解决这个问题！我之前对L矩阵的处理有误，哈哈。 - O.rka

你对如何在这个函数中使用广播来提高速度有什么建议吗？https://dev59.com/8rroa4cB1Zd3GeqPgTw2 - O.rka

Answer 2

给定邻接矩阵A，可以在不使用for循环的情况下计算TOM矩阵W，这将极大地加快进程。

L = np.dot(A,A)
k = np.sum(A,axis=0); d = len(k); tile = np.tile(k,(d,1))
K = np.min(np.stack((tile,tile.T),axis=2),axis=2)
W = (L + A)/(K + 1 - A); np.fill_diagonal(W, 1)

- Alon