数学推导与算法 · 生物医学网络科学

本节补全 WGCNA 中软阈值、TOM 和模块特征基因的数学细节，推导数据处理不等式、收缩强度和差异共表达检验，给出 GENIE3 与模块保守性检验的伪代码，并写出用 AUPR 评估推断网络的具体计算。符号： $n$ 个样本、 $p$ 个基因， $x_i\in\mathbb R^n$ 为基因 $i$ 的表达向量， $r_{ij}$ 为相关系数， $a_{ij}\in[0,1]$ 为邻接权重， $k_i=\sum_{u\ne i}a_{iu}$ 为连接度。

软阈值与 TOM 的性质

软阈值保持排序、拉开差距。 $a_{ij}=|r_{ij}|^\beta$ 是 $|r|$ 的单调增函数，所以不改变边的强弱排序；但两条边的比值被放大为 $(|r_1|/|r_2|)^\beta$ 。例如 $r=0.9$ 与 $r=0.5$ 的比值从 1.8 变成 $1.8^6\approx34$ 。这就是“软”阈值：没有一刀切，但弱边的贡献被压到可以忽略。

TOM 为什么在 $[0,1]$ 内。由于 $a_{uj}\le1$ ，

$l_{ij}=\sum_{u\ne i,j}a_{iu}a_{uj}\le\sum_{u\ne i,j}a_{iu}=k_i-a_{ij},$

同理 $l_{ij}\le k_j-a_{ij}$ ，所以 $l_{ij}+a_{ij}\le\min(k_i,k_j)$ ；而分母 $\min(k_i,k_j)+1-a_{ij}\ge\min(k_i,k_j)$ 。故 $0\le\mathrm{TOM}_{ij}\le1$ 。 $\mathrm{TOM}_{ij}=1$ 要求 $a_{ij}=1$ 且度较小的那个基因的所有邻居都是另一个基因的满权邻居； $\mathrm{TOM}_{ij}=0$ 要求二者既不相连也无共同邻居。在无权图上，TOM 就是“直接相连 + 共同邻居数”除以可能的最大值。

模块特征基因。把模块 $q$ 的 $m$ 个基因的表达逐基因标准化，得 $X^{(q)}\in\mathbb R^{n\times m}$ ，做奇异值分解 $X^{(q)}=UDV^\top$ 。特征基因取第一左奇异向量 $\mathrm{ME}_q=u_1$ （即第一主成分得分的方向），其解释的方差比例为 $d_1^2/\sum_k d_k^2$ 。奇异向量的符号不唯一，WGCNA 把 ME 的符号调整为与模块平均表达正相关。模块成员度 $k_{\mathrm{ME}}(i)=\operatorname{cor}(x_i,\mathrm{ME}_q)$ 本质上是基因 $i$ 在第一主成分上的（标准化）载荷。

若 $d_1^2/\sum d_k^2$ 很低（如低于 0.3），说明模块内部并不“一维”，用一个 ME 概括整个模块会丢失大量信息，ME–性状关联也就难以解释。

收缩强度的推导

设要估计的量为 $\sigma_{ij}$ （这里是真实相关 $\rho_{ij}$ ），无偏估计为 $s_{ij}$ （ $E s_{ij}=\sigma_{ij}$ ），目标 $t_{ij}$ 为常数（相关矩阵向单位阵收缩时，非对角 $t_{ij}=0$ ）。收缩估计 $s^*_{ij}=\lambda t_{ij}+(1-\lambda)s_{ij}$ 的总均方误差为

$\mathcal L(\lambda)=\sum_{i\neq j}E\big(s^*_{ij}-\sigma_{ij}\big)^2=\sum_{i\neq j}\Big[(1-\lambda)^2\operatorname{Var}(s_{ij})+\lambda^2(t_{ij}-\sigma_{ij})^2\Big].$

对 $\lambda$ 求导令其为 0： $-(1-\lambda)\sum\operatorname{Var}(s_{ij})+\lambda\sum(t_{ij}-\sigma_{ij})^2=0$ ，得

$\lambda^*=\frac{\sum_{i\neq j}\operatorname{Var}(s_{ij})}{\sum_{i\neq j}\big[\operatorname{Var}(s_{ij})+(t_{ij}-\sigma_{ij})^2\big]}=\frac{\sum_{i\neq j}\operatorname{Var}(s_{ij})}{\sum_{i\neq j}E(t_{ij}-s_{ij})^2}.$

最后一步用了 $E(t-s)^2=\operatorname{Var}(s)+(t-\sigma)^2$ 。把未知的期望换成样本估计：分母用 $\sum r_{ij}^2$ ，分子用样本乘积 $w_{kij}$ 的经验方差，就得到正文中的 Schäfer–Strimmer 公式（Ledoit & Wolf, 2004 给出了以单位阵倍数为目标的同类结果）。公式表明：估计越不确定、或真值离目标越近，就越应该收缩。

数据处理不等式（DPI）

互信息的链式法则把 $I(A;B,C)$ 按两种顺序展开：

$I(A;B,C)=I(A;C)+I(A;B\mid C)=I(A;B)+I(A;C\mid B).$

若 $A\to B\to C$ 构成马尔可夫链，则给定 $B$ 后 $A$ 与 $C$ 独立， $I(A;C\mid B)=0$ 。于是 $I(A;C)=I(A;B)-I(A;B\mid C)\le I(A;B)$ ；同理 $I(A;C)\le I(B;C)$ 。ARACNe 据此对每个三角形删去 MI 最小的边。

算法 1  ARACNe（简化版）
输入：表达矩阵 X（n×p）；MI 显著性阈值 α；DPI 容差 ε（0 表示严格 DPI，常用 0–0.15）
输出：无向网络
1. 对每对基因估计 I(i;j)（核密度或自适应分箱）              // O(p²·n)
2. 用置换得到零分布，删去 I(i;j) 不显著的边
3. 对每个三角形 (i,j,k)，若 I(i,j) < (1−ε)·min{I(i,k), I(j,k)}：标记 (i,j) 待删
4. 一次性删除所有标记边（先标记后删除，结果不依赖遍历顺序）
复杂度：第 3 步最坏 O(p³)；实际只需遍历第 2 步后剩余边构成的三角形

GENIE3

GENIE3 的边权是随机森林（或极端随机树）的方差减少重要性。在树的一个内部节点 $\mathcal N$ ，样本集 $S$ 被特征 $g$ 划分为 $S_l$ 与 $S_r$ ，该节点的贡献为

$\Delta(\mathcal N)=|S|\operatorname{Var}(S)-|S_l|\operatorname{Var}(S_l)-|S_r|\operatorname{Var}(S_r),$

把特征 $g$ 在所有树、所有节点上的 $\Delta$ 求和并取平均，即为 $g$ 对该靶基因的重要性。

算法 2  GENIE3
输入：表达矩阵 X（n×p，每个基因标准化为单位方差）；候选调控因子集合 R；
      树的数目 T（如 1000）；每次分裂考察的候选特征数 K（随机森林常取 sqrt(|R|)）
输出：按权重排序的有向边列表 (g → j, w_gj)
1. 对每个靶基因 j = 1..p：
2.     输入特征 ← R \ {j}；输出 ← x_j
3.     训练 T 棵树（随机森林：每棵树用 bootstrap 样本，每次分裂在随机选出的 K 个特征中
       找最优划分；极端随机树：用全部样本，分裂阈值也随机抽取）
4.     对每个 g ∈ R \ {j}：w_gj ← 所有树中以 g 分裂的节点的 Δ 之和 / T
5. 把全部 (g, j, w_gj) 汇总，按 w_gj 从大到小排序
复杂度：约 O(p · T · K · n log n)；对 p 个靶基因完全可并行

第 1 步中先把每个靶基因标准化为单位方差，是为了让不同靶基因的重要性可以放在一起排序。有向边的“方向”完全来自集合 $R$ 的先验设定。

差异共表达检验与模块保守性

Fisher $z$ 检验的依据：若 $(x_i,x_j)$ 近似二元正态， $\operatorname{artanh}r$ 近似服从 $\mathcal N\big(\operatorname{artanh}\rho,\,1/(n-3)\big)$ 。两组样本独立，所以差的方差为两组方差之和，标准化后在 $H_0:\rho^A=\rho^B$ 下近似标准正态。对 $p(p-1)/2$ 个基因对同时检验时，检验统计量之间高度相关（共享基因），BH 程序在正相关依赖下仍控制 FDR，但功效有限，因此常把检验聚合到基因或模块层面。

算法 3  模块保守性的置换 Z 统计量（简化版）
输入：参考数据相关矩阵 R_ref；检验数据 X_test；模块 M（大小 m）；置换次数 B
1. 计算观测统计量：密度类 D = (meanCor, meanAdj, propVarExpl, meanKME)；
   连接模式类 C = (cor.kIM, cor.kME, cor.cor)，均在检验数据中、以参考数据为对照计算
2. 对 b = 1..B：在检验数据中随机抽取 m 个基因作为“伪模块”（置换模块标签），重算 D_b、C_b
3. 每个统计量的 Z = (观测值 − 置换均值) / 置换标准差
4. Zdensity = median(Z_D)；Zconnectivity = median(Z_C)；Zsummary = (Zdensity + Zconnectivity) / 2
复杂度：O(B · m² · n)；B 取 100–200 通常足以估计均值与标准差

AUPR 的计算与基线

把推断结果排成边列表 $e_1,e_2,\dots$ ，金标准真边集合为 $G$ ，评估的候选边全集为 $\mathcal C$ 。前 $k$ 条预测中真阳性数为 $\mathrm{TP}_k$ ：

$P_k=\frac{\mathrm{TP}_k}{k},\qquad R_k=\frac{\mathrm{TP}_k}{|G|},\qquad \mathrm{AUPR}\approx\sum_k (R_k-R_{k-1})P_k .$

随机排序的期望精确率处处等于真边比例 $\pi=|G|/|\mathcal C|$ ，因此随机 AUPR $\approx\pi$ 。DREAM5 中 $\pi$ 常远低于 1%，所以“AUPR=0.05”可能已是随机水平的几倍，而 AUROC 在这种极端不平衡下会显得很高、区分度低。实践要点： $\mathcal C$ 只应包含金标准覆盖的调控因子与靶基因之间的边；PR 曲线不能在相邻点之间线性插值（Davis & Goadrich, 2006）；报告 AUPR 时应同时报告 $\pi$ 或 AUPR/ $\pi$ 。

不完整金标准的影响：设真实真边集合为 $G^*$ ，金标准只覆盖其中比例 $c$ （ $G\subset G^*$ ， $|G|=c|G^*|$ ），且漏标与预测得分无关。则前 $k$ 条预测中被记为真阳性的期望数为 $c\cdot\mathrm{TP}^*_k$ ，表观精确率 $P_k\approx c\,P^*_k$ ，随机基线也同比变为 $c\,\pi^*$ 。所以比值 AUPR/ $\pi$ 对随机漏标近似不变，而 AUPR 本身被按比例低估；若漏标与得分相关（例如数据库偏好研究充分的高表达基因，而方法也偏好它们），比值也会有偏。