NetMoss2微生物组网络分析
Xiao, L., Zhang, F. & Zhao, F. Large-scale microbiome data integration enables robust biomarker identification. Nat Comput Sci 2, 307–316 (2022). https://doi.org/10.1038/s43588-022-00247-8
10.6.1 NetMoss2论文解析
Network Module Structure Shift (NetMoss): 网络模块结构转换算法,该算法侧重于网络模块的转换,以评估细菌在不同状态之间的重要性。
因为大型研究构建的网络表现出微弱的微生物相互作用,将不同大小的数据集直接整合到一个网络中可能会掩盖大型数据集的真实微生物特征。为了解决这个问题,我们在分析中引入了**一种单变量加权方法(a univariate weighting method)**, 即为较大的数据集分配更大的权重,以增加其在最终集成网络中的贡献。我们首先在成对排列测试中验证了该方法,其中任意两个网络被集成到一个网络中,如图2b所示,在两个研究组的所有组合中,大型研究对整合网络的贡献更大,贡献随着纳入的样本量增加而增加。在所有七个网络中的整合中观察到类似的结果: 社区规模越大,其对最终整合网络的贡献越大,表明这种单变量权重方法可以有效地突出大型研究在最终网络中的优势,并减少集成过程中的偏差。
问题:
- **一种单变量加权方法(a univariate weighting method)**的具体含义是什么?
- 图2g中的网络距离是如何度量的?Combat和Limma是如何去除批次效应的?
- 预测过渡的过程中,如何生成模拟网络?为什么NetMoss算法能够测量不同状态之间网络结构的变化呢?
- Neighbor Shift (NESH) 分数和 Jaccard Edge Index (JEI)是如何定义的?其如何衡量网络中节点的变化?
- NetMoss如何判断其识别的生物标记物是准确的?并且效果比Wilcox.test效果好?
- 从 gutMDisorder 数据库中检索了 66 种 CRC 相关细菌,发现在疾病和对照组中,标记细菌的连接强度显著高于非标记细菌。
- 66个已知的CRC相关细菌中,有55个出现在整合的CRC数据集中,其中传统的统计检验方法只能识别32%的标记细菌,而NetMoss识别了其中68%的细菌标记物。
- 图4f-i中的AUC是如何计算的?是数据集内部的交叉验证?还是迁移的结果?
10.6.2 Network analysis
微生物的共生网络是用 SPIEC-EASI 38构建的。使用R包“igraph” 计算网络的拓扑系数。对于每个病例或对照组,使用单变量加权方法去除批次效应, 不同规模的研究被整合到一个网络中。考虑到由大型研究构建的网络表现出弱的共丰度模式,向大型网络添加了很大的权重以增加它们对最终集成网络的贡献,从而减少集成中的偏差。该过程如下实施。
首先,基于每项研究的丰度矩阵构建了一个共现网络。接下来,比对任意两个分类群之间的每一对共丰度模式。缺失的共丰度模式用值 0 填充。实施单变量加权方法以根据每个研究的大小为不同的共同丰度模式对添加不同的权重。在此过程中,Hedges 和 Olkin 方法40用于评估每项研究中相关系数的条件偏差。对于某项研究 n i n_i ni,相关系数 r i r_i ri的条件偏差 v i v_i vi计算为:
补充:
Hedges 和 Olkin 方法是做什么的?有何意义?
回答: 估计的线性加权组合。当k项独立研究共享相同的效应量 σ \sigma σ时,很自然地通过汇集每项研究的估计值来估计 σ \sigma σ。如果研究的样本量不同,则较大的研究的估计值比较小研究的估计值更加精确。在这种情况下,在合并时为了更精确的估计赋予较大的研究更多权重时合理的。这导致加权估计的形式为:
d w = w 1 d 1 + . . . + w k d k d_w = w_1d_1 + ... + w_kd_k dw=w1d1+...+wkdk
其中 w 1 , . . . , w k w_1,...,w_k w1,...,wk时非负权重且和为1。估计加权:使 d w d_w dw方差最小化的权重给出的权重与每个研究中的方差成反比【即方差越大的研究,在最终线性组合中所占权重越小】。这在直觉上是清楚的,因为更小的方差即意味着更精确,应该占据更大的权重。因而,
w i = 1 σ 2 d i / ∑ j = 1 k 1 σ 2 ( d j ) w_i = \frac{1}{\sigma^2{d_i}}/\sum_{j=1}^k \frac{1}{\sigma^2(d_j)} wi=σ2di1/j=1∑kσ2(dj)1
其中 σ 2 ( d i ) \sigma^2(d_i) σ2(di)是 d i d_i di的方差。因为我们使用了大样本理论,权重是
w i = 1 σ 2 d i / ∑ j = 1 k 1 σ ∞ 2 ( d j ) w_i = \frac{1}{\sigma^2{d_i}}/\sum_{j=1}^k \frac{1}{\sigma^2_{\infin}(d_j)} wi=σ2di1/j=1∑kσ∞2(dj)1
其中 σ ∞ 2 ( d i ) \sigma_{\infin}^2(d_i) σ∞2(di)是大样本方差。计算最精确加权估计的实际问题是第i个权重取决于 d i d_i di的方差,而 d i d_i di又取决于未知参数 σ \sigma σ,并且必须对其进行估计。
为了解决这个困难,可以使用不依赖于 S S S的方差近似值的权重。此过程会产生无偏的合并估计量,并且通常不如使用最佳权重时精确。例如,近似权重由下式给出:
w i ≅ n i ~ ∑ j = 1 k n j ~ w_i \cong \frac{\tilde{n_i}}{\sum_{j=1}^k \tilde{n_j}} wi≅∑j=1knj~ni~
其中 n ~ i = n i E n i C n i E + n i C \tilde n_i = \frac{n_i^En_i^C}{n_i^E + n_i^C} n~i=niE+niCniEniC。当 σ \sigma σ接近于0或者 n i ~ \tilde{n_i} ni~很大时,由此得出的权重接近最优。
v i = 1 − r i 2 n i − 1 v_i = \frac{1-r_i^2}{n_i -1 } vi=ni−11−ri2
每对共丰度模式的权重定义为:
ρ
=
∑
i
=
1
k
w
i
r
i
∑
i
=
1
k
w
i
\rho = \frac{\sum_{i=1}^{k}w_ir_i}{\sum_{i=1}^{k}w_i}
ρ=∑i=1kwi∑i=1kwiri
其中
w
i
w_i
wi是
v
i
v_i
vi的倒数,k代表研究的数量。
为了展示真实的生态过程,使用 WGCNA(Weighted correlation network analysis, 加权共表达网络分析, http://blog.genesino.com/2018/04/wgcna/)进行了模块划分,微生物在一个模块中相互合作,同时在任意两个模块之间竞争性交互。在加权网络中,节点
i
i
i的连接强度定义为该节点与网络中所有其他节点之间的连接强度之和,如:
k
i
=
∑
j
=
1
n
a
i
j
k_i = \sum_{j=1}^n a_{ij}
ki=j=1∑naij
其中
a
i
j
a_{ij}
aij表示节点
i
i
i和节点
j
j
j之间的相关系数。
为了突出节点在网络模块结构中的重要性,我们将节点
i
i
i的连接强度重新定义为:
k
i
=
∑
j
=
1
k
a
i
j
−
∑
n
1
+
1
n
a
i
j
k_i = \sum_{j=1}^ka_{ij} - \sum_{n_1+1}^na_{ij}
ki=j=1∑kaij−n1+1∑naij
其中n代表网络中所有节点的数量,
n
i
n_i
ni代表特定模块内的节点数量。
10.6.3 NetMoss算法
NetMoss score 用于衡量每个节点在网络结构转变中的驱动力。 这被定义如下。
首先,健康状态的相关矩阵为
A
=
[
c
i
j
]
A=[c_{ij}]
A=[cij],疾病状态的相关矩阵为
A
′
=
[
c
′
i
j
]
A^′=[c^′{ij}]
A′=[c′ij],其中
c
i
j
c_{ij}
cij为节点
i
i
i和节点
j
j
j的相关系数:
c
i
j
=
c
o
r
(
i
,
j
)
c_{ij}=cor(i,j)
cij=cor(i,j)
为了获得优化的模块结构,实施了线性变换以将
c
i
j
c_{ij}
cij转换为
s
i
j
s_{ij}
sij
s
i
j
=
1
+
c
i
j
2
s_{ij} = \frac{1+c_{ij}}{2}
sij=21+cij
因此,变换后健康和疾病状态的相关矩阵分别为
B
=
[
s
i
j
]
B =[s_{ij}]
B=[sij]和
B
′
=
[
s
i
j
′
]
B'=[s'_{ij}]
B′=[sij′]。
基于WGCNA算法计算软阈值
β
\beta
β,加权网络
a
i
j
a_{ij}
aij为:
a
i
j
=
∣
s
i
j
∣
β
a_{ij}=|s_{ij}|^β
aij=∣sij∣β
因此,健康和疾病状态的加权相关矩阵分别为 $C =[a_{ij}] $和
C
′
=
[
a
i
j
′
]
C'=[a'_{ij}]
C′=[aij′]。
计算节点 i i i 和节点$ j 的拓扑重叠矩阵 的拓扑重叠矩阵 的拓扑重叠矩阵 ω_{ij}$ 为
ω
i
j
=
l
i
j
+
a
i
j
m
i
n
{
k
i
,
k
j
}
+
1
−
a
i
j
ω_{ij}=\frac{l_{ij}+a_{ij}}{min\{k_i,k_j\}+1−a_{ij}}
ωij=min{ki,kj}+1−aijlij+aij
其中
l
i
j
=
Σ
u
a
i
u
a
u
j
l_{ij}=Σ_ua_{iu}a_{uj}
lij=Σuaiuauj,
k
i
=
Σ
a
i
u
ki=Σa_{iu}
ki=Σaiu;
u
u
u 表示网络中除节点
i
i
i 和节点
j
j
j 之外的其他节点。【如何理解呢?
l
i
j
为其他节点与
i
和
j
之间的权重的乘积之和。
l_{ij}为其他节点与i和j之间的权重的乘积之和。
lij为其他节点与i和j之间的权重的乘积之和。】
拓扑重叠矩阵(https://zhuanlan.zhihu.com/p/441952423):
- Topological overlap(TO,拓扑重叠):拓扑重叠是通过比较两个节点和网络中其他节点的加权相关性来定量描述节点之间相似性的方法。
- Topological overlap matrix(TOM,拓扑重叠矩阵):把邻接矩阵转换为拓扑重叠矩阵,以降低噪音和假相关,获得的新距离矩阵,这个信息可拿来构建网络或绘制TOM图。
- Adjacency matrix(邻接矩阵):基因和基因之间的加权相关值构成的矩阵。邻接矩阵由分布在0到1之间的数值组成,用来描述节点间相关性强度,由相关性矩阵根据power值进行次方转化而来。
节点$ i 和节点 和节点 和节点 j$ 之间的距离定义为
d
i
j
=
1
−
ω
i
j
d_{ij}=1−ω_{ij}
dij=1−ωij
然后,健康和疾病状态的距离矩阵分别为
D
=
[
d
i
j
]
D =[d_{ij}]
D=[dij]和
D
′
=
[
d
i
j
′
]
D'=[d'_{ij}]
D′=[dij′]。
模块划分在矩阵
D
D
D和
D
′
D'
D′中进行,矩阵
D
D
D包含
n
n
n个模块,矩阵
D
′
D'
D′包含
m
m
m个模块。
D
D
D和
D
′
D'
D′的交集模块数为
K
(
K
≤
m
n
)
K(K≤ mn)
K(K≤mn)。 计算交叉模块各节点到健康距离矩阵D和疾病距离矩阵D′的平均距离,分别得到N×K阶矩阵
D
m
o
d
D_{mod}
Dmod和
D
m
o
d
′
D^′_{mod}
Dmod′。 差分模块距离矩阵定义为:
Δ
D
=
D
m
o
d
′
−
D
m
o
d
ΔD=D^′_{mod}−D_{mod}
ΔD=Dmod′−Dmod
任意交叉模块中节点
i
i
i 的 NetMoss 得分为
N
M
S
S
(
i
)
A
→
B
=
∑
j
N
e
i
g
h
b
o
r
s
A
Δ
D
i
j
−
∑
l
N
e
i
g
h
b
o
r
s
A
Δ
D
i
l
NMSS(i)_{A\rightarrow B} = \sum_j^{Neighbors A}\Delta D_{ij} - \sum_l^{Neighbors A}\Delta D_{il}
NMSS(i)A→B=j∑NeighborsAΔDij−l∑NeighborsAΔDil
其中 A 和 B 分别代表健康和疾病网络; NeighborsA 表示健康网络中的所有相邻模块,NeighborsB 表示疾病网络中的所有相邻模块。
交叉模块代表从健康到疾病转变期间的稳定元素,其中转变的模块导致网络结构的改变。 NetMoss算法通过衡量网络结构转变的驱动力来评估每个节点的重要性。
10.6.4 模块转换模拟和随机噪声产生
为了验证NetMoss算法在识别网络结构方面的效果,生成了模拟网络。考虑到微生物群网络的稀疏性,我们开发了一种算法来生成具有一定模块结构的模拟相关矩阵,并在矩阵中添加随机噪声以模拟自然干扰。
从向量空间
R
M
k
R^{Mk}
RMk中选取
g
k
g_k
gk个单位向量,构造一个
M
k
−
b
y
−
g
k
M_k-by-g_k
Mk−by−gk矩阵
U
k
U_k
Uk。 第
k
k
k 个模块由
g
k
−
b
y
−
g
k
gk-by-gk
gk−by−gk 矩阵
Σ
k
Σk
Σk 表示:
U
k
=
(
u
1
∣
u
2
∣
.
.
.
∣
u
k
)
U_k = (u_1|u_2|...|u_k)
Uk=(u1∣u2∣...∣uk)
∑ k = ρ k ( U k T U k ) + ( 1 − ρ k ) I \sum_k = \rho_k(U_k^TU_k) + (1-\rho_k)I k∑=ρk(UkTUk)+(1−ρk)I
其中 k k k代表着模拟矩阵中的模块数目, g k g_k gk表示第 k k k个模块的大小, M k M_k Mk表示第 k k k个模块内部相关系数的变化范围, ρ k ρk ρk表示第 k k k个模块内部的最大相关系数, I I I表示单位向量。
据此构造了一个 N × N N×N N×N的矩阵,模块 Σ 1 , Σ 2 , . . . , Σ k Σ_1,Σ_2,...,Σ_k Σ1,Σ2,...,Σk依次排列在对角线上,模块以外的区域用 0 0 0填充。
为了产生随机噪声,从向量空间
R
m
R^m
Rm 中选择
c
c
c 个单位向量来构造一个
m
×
c
m×c
m×c 矩阵
U
k
U_k
Uk。 随机噪声矩阵
S
S
S 选自
c
×
c
c×c
c×c 矩阵
S
k
S_k
Sk:
U
k
=
(
u
1
∣
u
2
∣
…
∣
u
c
)
S
k
=
ε
k
(
U
k
T
U
k
)
ε
k
=
ε
1
+
k
−
1
r
−
1
(
ε
r
−
ε
1
)
U_k=(u_1|u_2|…|u_c) \\ S_k=ε_k(U^T_kU_k)\\ ε_k=ε_1+\frac{k−1}{r−1}(ε_r−ε_1)
Uk=(u1∣u2∣…∣uc)Sk=εk(UkTUk)εk=ε1+r−1k−1(εr−ε1)
其中
r
r
r表示噪声矩阵的行,
c
c
c表示噪声矩阵的列,
m
m
m表示噪声的变化范围,
ε
k
ε_k
εk表示
k
k
k阶噪声,其中
k
=
1
k=1
k=1表示噪声的最小值,
k
=
r
k=r
k=r 表示噪声的最大值。 最后,在相应的模块矩阵中加入随机噪声矩阵
S
S
S来模拟自然干扰。
ROC曲线的绘制:https://blog.csdn.net/duyibo123/article/details/110090088