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

Question

我正在尝试计算邻接矩阵的加权拓扑重叠，但我不知道如何使用 numpy 正确地计算它。执行正确实现的 R 函数来自 WGCNA ( https://www.rdocumentation.org/packages/WGCNA/versions/1.67/topics/TOMsimilarity )。计算这个(我认为)的公式详见equation 4我相信在下面正确地转载了它。

有谁知道如何正确实现以反射(reflect) 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 的给定描述说明:

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 .

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

应用公式时仍然缺少的是对角线值不匹配。我默认将它们设置为一个。无论如何，节点与自身的相互联系是什么？不同于一(或零，取决于定义)的值对我来说没有意义。案例零 和案例一 都不会导致 w_ii=1在简单的例子中。在零情况中，k_i+1 == l_ii 是必要的，在案例一中，k_i == l_ii+1 是必要的，这对我来说都是错误的。

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