数学推导与算法

本节把邻近度、分离度写成可执行的算法，说明度匹配抽样的细节与复杂度，并给出 DIAMOnD 的伪代码和基于扩散的替代方法——带重启的随机游走（RWR）与热核。符号：相互作用组 $G=(V,E)$，$|V|=N$，$|E|=M$；$d(u,v)$ 为最短路径长度；$S$ 为疾病基因集合，$T$ 为药物靶点集合；$k_v$ 为节点度。

最近距离 $d_c$ 与度匹配 $z$ 分数

$d_c(S,T)=\frac1{|T|}\sum_{t\in T}\min_{s\in S}d(s,t)$ 有一个计算上的便利：$\min_{s\in S}d(s,t)$ 对所有 $t$ 可以用一次多源广度优先搜索（把 $S$ 中所有节点同时作为起点、距离设为 0）得到，代价 $O(N+M)$，不需要计算全体节点对的距离。注意 $d_c$ 不对称：它对靶点求平均，衡量“每个靶点离疾病模块多远”。

算法 1  网络邻近度 z 分数（Guney 等, 2016 思路）
输入：G；疾病基因 S；靶点 T；随机次数 B（如 1000）；分箱最小容量 b（例如约 100 个节点）
输出：d_c、z 分数、经验 p 值
1. S, T ← 只保留位于 G 最大连通分支中的节点（报告丢弃了多少）
2. 度分箱：把节点按度升序排列，依次装箱，每箱至少 b 个节点（高度端的箱会跨越较宽的度范围）
3. dist ← 以 S 为源的多源 BFS；d_obs ← mean_{t∈T} dist[t]
4. 对 i = 1..B：
5.     S' ← 对 S 中每个节点，从其所在度箱中无放回随机抽一个替身
6.     T' ← 对 T 同样处理
7.     dist' ← 以 S' 为源的多源 BFS；d_i ← mean_{t∈T'} dist'[t]
8. μ, σ ← d_1..d_B 的均值与标准差
9. z ← (d_obs − μ)/σ；p ← (1 + #{d_i ≤ d_obs}) / (B + 1)
复杂度：O(B·(N+M))；对 1.3 万节点、14 万边的互作组，B=1000 时在普通电脑上可接受

为什么必须分箱抽样？最短路径长度强烈依赖度：高度节点到任意节点的平均距离更短。如果 $T$ 中多为枢纽而随机替身均匀抽取，替身的度偏低、距离偏长，$z$ 就会系统性地偏负，产生“药物特异地靠近疾病”的假象。分箱是一种折中：箱太细则替身选择太少、随机性不足；箱太粗则度匹配不准。

分离度 $s_{AB}$ 的计算与含义

按 Menche 等（2015）的定义：

$$\langle d_{AA}\rangle=\frac1{|A|}\sum_{a\in A}\min_{a'\in A,\,a'\ne a}d(a,a'),\qquad \langle d_{AB}\rangle=\frac{1}{|A|+|B|}\Big[\sum_{a\in A}\min_{b\in B}d(a,b)+\sum_{b\in B}\min_{a\in A}d(a,b)\Big].$$

$\langle d_{AB}\rangle$ 对 $A$、$B$ 对称，共有基因的距离为 0。可以这样理解 $s_{AB}=\langle d_{AB}\rangle-\frac12(\langle d_{AA}\rangle+\langle d_{BB}\rangle)$ 的符号：若 $A$、$B$ 是同一个模块被随机分成两半，那么从 $A$ 的基因找最近的 $B$ 基因，与找最近的另一个 $A$ 基因，距离分布相近，$s_{AB}\approx0$；若两模块有共享基因或彼此穿插，跨模块的最近距离比模块内更短，$s_{AB}<0$；若两模块在网络中各居一处，$s_{AB}>0$。$\langle d_{AB}\rangle$ 可用两次多源 BFS 得到；$\langle d_{AA}\rangle$ 需要排除自身，可从每个 $a$ 出发做 BFS，代价 $O(|A|(N+M))$。

DIAMOnD 的伪代码

算法 3  DIAMOnD 疾病模块扩展（Ghiassian 等, 2015）
输入：G（N 个节点）；种子集合 S0；迭代次数 n；种子权重 α（默认 1）
输出：按加入顺序排列的节点列表
1. S ← S0
2. 重复 n 次：
3.     C ← S 中节点的全部邻居 − S
4.     对每个 v ∈ C：
5.         k ← 节点 v 的度；ks ← v 在 S 中的邻居数（种子计 α 份时，相应调整 N、|S|、k、ks）
6.         p(v) ← Σ_{j=ks}^{k} C(|S|, j)·C(N−|S|, k−j) / C(N, k)      # 超几何上尾
7.     v* ← argmin p(v)；S ← S ∪ {v*}；记录 v*
8. 返回记录的顺序；模块边界由留出种子的召回率或独立注释的富集决定
复杂度：每次迭代需更新候选的 ks（只有 v* 的邻居会变），用对数阶乘计算 p 值；总计约 O(n·(|C| + k_max))

超几何 $p$ 值在 $k$ 和 $k_s$ 都较大时数值很小，应在对数空间计算（如 scipy.stats.hypergeom.logsf）以免下溢。

扩散方法：带重启的随机游走

最短路径只看“最近的一条路”，对缺失一条边非常敏感。扩散方法综合了所有路径。设 $A$ 为邻接矩阵，$D=\operatorname{diag}(k_v)$，列随机矩阵 $W=AD^{-1}$。游走者每一步以概率 $1-r$ 走到随机邻居，以概率 $r$ 跳回种子集合，种子分布 $p_0$（如 $S$ 上均匀，其余为 0）：

$$p_{t+1}=(1-r)\,W p_t+r\,p_0 .$$

因为 $W$ 的谱半径为 1，$(1-r)W$ 的谱半径小于 1，迭代收敛到唯一的不动点

$$p^\ast=r\big(I-(1-r)W\big)^{-1}p_0=r\sum_{\ell=0}^{\infty}(1-r)^\ell W^\ell p_0 ,$$

即所有长度 $\ell$ 的游走按几何权重 $(1-r)^\ell$ 加总：$r$ 越大越“局部”，$r$ 越小越“全局”。

算法 2  RWR 候选基因 / 药物打分
输入：邻接矩阵 A（稀疏）；种子集合 S；重启概率 r；收敛阈值 tol（如 1e-6）
输出：每个节点的稳态分数 p*
1. W ← A D^{-1}（孤立节点需先移除）
2. p0 ← 1/|S| 于 S 上，其余为 0；p ← p0
3. 重复：p_new ← (1−r)·W·p + r·p0
4.     若 ||p_new − p||_1 < tol：停止；否则 p ← p_new
5. 返回 p*（排序时排除种子本身）
6. （去度偏倚）以度匹配随机种子集重复 1–5 共 B 次，
   把每个节点的分数转为 z 分数或经验百分位
复杂度：每次迭代 O(M)，迭代次数约 log(tol)/log(1−r)；总计 O(B·M·log(1/tol)/r)

热核与 RWR 的统一

热扩散 $\dot p=-Lp$（$L=I-W$）的解是 $p(t)=e^{-tL}p_0$。由于 $I$ 与 $W$ 可交换， $$e^{-tL}=e^{-t}e^{tW}=\sum_{\ell\ge0}\frac{e^{-t}t^\ell}{\ell!}W^\ell .$$ 与 RWR 的 $r\sum_\ell(1-r)^\ell W^\ell$ 对比可见，两者都是 $\sum_\ell c_\ell W^\ell p_0$ 形式的游走长度加权和，权重分别为泊松分布与几何分布。几何分布的尾部更重，RWR 会给更长的游走留下更多权重；泊松分布集中在 $t$ 附近。对称归一化下的同类核还包括正则化拉普拉斯核 $(I+\beta L_{\mathrm{sym}})^{-1}$ 与扩散核 $e^{-\beta L_{\mathrm{sym}}}$（$L_{\mathrm{sym}}=I-D^{-1/2}AD^{-1/2}$）。大网络中无须显式计算矩阵指数，可用 Krylov 子空间方法直接求 $e^{-tL}p_0$（如 scipy.sparse.linalg.expm_multiply），RWR 则解一个稀疏线性方程组。