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

Question

我正在尝试计算邻接矩阵的加权拓扑重叠，但我不知道如何使用 numpy 正确地计算它。正确实现的 R 函数来自 WGCNA（https://www.rdocumentation.org/packages/WGCNA/versions/1.67/topics/TOMsimilarity). The formula for computing this (I THINK) is detailed in equation 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 版本的源代码是 here 但它使用 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) but it is not for weighted adjacencies. The author of GTOM also took a look at this problem（sophisticated/efficient NumPy 实现了很多，但它仍然没有产生预期的结果）。

有谁知道如何重现 WGCNA 实现？

编辑：2019.06.20 我改编了来自@scleronomic 的一些代码和来自 v2016.06 的 @benmaier with credits in the doc string. The function is available in soothsayer 等。希望这将使人们能够更轻松地在 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

Answer 1

首先让我们看一下二元邻接矩阵的方程部分 a_ij:

a_ij：表示节点i是否连接到节点j
k_i：节点i的邻居数（连通性）
l_ij：节点i和节点j

so w_ij 测量具有较低连通性的节点的邻居中有多少也是另一个节点的邻居（即 w_ij 测量“它们的相对 inter-connectedness”） .

我的猜测是他们将 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 的对角线应该为零而不是 1 的直觉：

 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 的给定描述说明：

Note that w_ij = 1 if the node with fewer connections satisfies two conditions:

(a) all of its neighbors are also neighbors of the other node and

(b) it is connected to the other node.

In contrast, w_ij = 0 if i and j are un-connected and the two nodes do not share any neighbors.

所以 A-D 之间的连接应该满足这个标准并且是 w_14=1.

案例零对角线：
案例一对角线：

应用公式时仍然缺少的是对角线值不匹配。我默认将它们设置为一个。无论如何，节点自身的 inter-connectedness 是什么？不同于一（或零，取决于定义）的值对我来说没有意义。 Case Zero 和 Case One 都不会导致简单示例中的 w_ii=1。在案例零中，k_i+1 == l_ii是必要的，而在案例一中，k_i == l_ii+1是必要的我觉得两者都不对。

总而言之，我会将邻接矩阵的对角线设置为 zero，使用给定的方程并将结果的对角线默认设置为 one。

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)

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

How to compute the Topological Overlap Measure [TOM] for a weighted adjacency matrix in Python?

python

arrays

r

statistics

adjacency-matrix