生物医学网络科学Biomedical Network Science
Example01 / 09

软阈值幂次 β 的无标度拟合准则

所属专题
基因调控网络推断与共表达分析
下载
Notebook · .py

对 β = 1..12 构建无符号软阈值邻接矩阵 A = |r|^β，计算连接度 k_i = Σ_j a_ij（不含自身）， 把 k 分 10 个箱，拟合 log10 p(k) ~ log10 k，拟合指数 = −sign(斜率)·R²；选拟合指数首次 ≥ 0.85 的最小 β。 输出解读：β 增大时拟合指数先快速上升、随后在高位波动，平均连接度单调下降。 拟合指数只是选 β 的经验准则，不是对网络无标度性的检验。

Output输出图

软阈值幂次 β 的无标度拟合准则：输出图 1

Output文字输出

beta =  1：拟合指数 = -0.17，平均连接度 =  21.68
beta =  2：拟合指数 =  0.12，平均连接度 =   8.18
beta =  3：拟合指数 =  0.87，平均连接度 =   4.45
beta =  4：拟合指数 =  0.84，平均连接度 =   2.84
beta =  5：拟合指数 =  0.81，平均连接度 =   1.95
beta =  6：拟合指数 =  0.89，平均连接度 =   1.40
beta =  7：拟合指数 =  0.92，平均连接度 =   1.03
beta =  8：拟合指数 =  0.89，平均连接度 =   0.77
beta =  9：拟合指数 =  0.94，平均连接度 =   0.59
beta = 10：拟合指数 =  0.91，平均连接度 =   0.45
beta = 11：拟合指数 =  0.88，平均连接度 =   0.36
beta = 12：拟合指数 =  0.87，平均连接度 =   0.28
选定 beta = 3

Code代码

# 对 β = 1..12 构建无符号软阈值邻接矩阵 A = |r|^β，计算连接度 k_i = Σ_j a_ij（不含自身），
# 把 k 分 10 个箱，拟合 log10 p(k) ~ log10 k，拟合指数 = −sign(斜率)·R²；选拟合指数首次 ≥ 0.85 的最小 β。
# 输出解读：β 增大时拟合指数先快速上升、随后在高位波动，平均连接度单调下降。
# 拟合指数只是选 β 的经验准则，不是对网络无标度性的检验。
import numpy as np
from scipy import stats
import matplotlib.pyplot as plt

def simulate(seed=13):
    """模拟表达矩阵 X（样本 x 基因）、真实模块标签、真实载荷、临床性状。"""
    rng = np.random.default_rng(seed)
    n, sizes = 60, [30, 25, 20, 15]
    trait = rng.normal(size=n)
    factors = [0.8 * trait + 0.6 * rng.normal(size=n)] + [rng.normal(size=n) for _ in sizes[1:]]
    X, labels, load = [], [], []
    for m, (s, f) in enumerate(zip(sizes, factors)):
        for l in rng.uniform(0.5, 0.95, size=s):          # 载荷越大，越像模块“核心”
            X.append(l * f + np.sqrt(1 - l**2) * rng.normal(size=n))
            labels.append(m + 1); load.append(l)
    for _ in range(30):                                     # 噪声基因，标签 0
        X.append(rng.normal(size=n)); labels.append(0); load.append(0.0)
    return np.array(X).T, np.array(labels), np.array(load), trait

def scale_free_fit(A, nbins=10):
    """返回 (无标度拟合指数, 平均连接度)。A 的对角线为 1，计算连接度时减去。"""
    k = A.sum(axis=1) - 1.0
    cnt, edges = np.histogram(k, bins=nbins)
    mid = (edges[:-1] + edges[1:]) / 2
    ok = cnt > 0
    res = stats.linregress(np.log10(mid[ok]), np.log10(cnt[ok] / cnt.sum()))
    return -np.sign(res.slope) * res.rvalue**2, k.mean()

X, true_lab, true_load, trait = simulate()
R = np.corrcoef(X.T)                                        # 基因 x 基因相关矩阵

betas = list(range(1, 13))
res = [scale_free_fit(np.abs(R) ** b) for b in betas]
fits = [r[0] for r in res]; mean_k = [r[1] for r in res]
for b, f, mk in zip(betas, fits, mean_k):
    print(f"beta = {b:2d}：拟合指数 = {f:5.2f}，平均连接度 = {mk:6.2f}")
beta_star = next(b for b, f in zip(betas, fits) if f >= 0.85)
print("选定 beta =", beta_star)

fig, ax = plt.subplots(1, 2, figsize=(9, 3.5))
ax[0].plot(betas, fits, "o-"); ax[0].axhline(0.85, color="r", ls="--")
ax[0].set_xlabel("beta"); ax[0].set_ylabel("scale-free fit index")
ax[1].plot(betas, mean_k, "o-"); ax[1].set_xlabel("beta"); ax[1].set_ylabel("mean connectivity")
plt.tight_layout(); plt.show()

Run在浏览器中运行

首次加载 Python 环境约需 20–40 秒；可以修改参数后重新运行。

软阈值幂次 β 的无标度拟合准则