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

Question

我正在尝试计算邻接矩阵的加权拓扑重叠，但我无法弄清楚如何使用numpy 正确计算。正确实现的R 函数来自WGCNA (https://www.rdocumentation.org/packages/WGCNA/versions/1.67/topics/TOMsimilarity)。 equation 4 中详细说明了计算此 (I THINK) 的公式，我相信它在下面正确复制。

有谁知道如何正确实现它以反映 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) 的包，但它不适用于加权邻接。 The author of GTOM also took a look at this problem（这是一个更加复杂/高效的 NumPy 实现，但仍然没有产生预期的结果）。

有人知道如何重现 WGCNA 实现吗？

编辑：2019.06.20 我已经改编了来自@scleronomic 和@benmaier 的一些代码，并在文档字符串中添加了学分。该功能在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

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 的对角线应该为零而不是 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的给定描述解释：

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

• (a) 它的所有邻居也是另一个节点的邻居，并且
• (b) 它连接到另一个节点。

相比之下，w_ij = 0 如果 i 和 j 未连接且两个节点不共享任何邻居。

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

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

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