互信息网络、CLR 背景校正与 ARACNe 的数据处理不等式

# 在一个已知的有向无环调控网络（30 个基因，含线性、饱和与非单调 U 形调控）上生成 400 个稳态样本，
# 分别用 |皮尔逊相关|、互信息（等频分箱的插入估计）、CLR 与 ARACNe（DPI，容差 0）给所有基因对打分，
# 以真实的“父–子”对为金标准（无向），计算 AUROC、AUPR、早期精度比 EPR 与随机基线。
# 输出解读：U 形调控使相关几乎为 0，互信息仍能发现；CLR 与 DPI 去除部分间接边，早期精度通常提高。
# 结果依赖模拟设定、分箱数与样本量；不要把这里的排名推广到真实数据。
import numpy as np

def simulate_dag(p=30, n=400, seed=1):
    rng = np.random.default_rng(seed)
    X = np.zeros((n, p)); edges = set(); kinds = {}
    X[:, :4] = rng.normal(size=(n, 4))                      # 4 个根节点（上游调控因子）
    for j in range(4, p):
        parents = rng.choice(j, size=rng.integers(1, 3), replace=False)
        signal = np.zeros(n)
        for k in parents:
            kind = rng.choice(["linear", "saturating", "u-shape"], p=[0.4, 0.4, 0.2])
            z = (X[:, k] - X[:, k].mean()) / X[:, k].std()
            signal += {"linear": z, "saturating": 2 * np.tanh(1.5 * z), "u-shape": z**2 - 1}[kind]
            edges.add((min(j, k), max(j, k))); kinds[(min(j, k), max(j, k))] = kind
        X[:, j] = signal + 0.7 * signal.std() * rng.normal(size=n)
    return X, edges, kinds

def mutual_information(X, bins=8):
    """等频分箱后的插入（plug-in）互信息估计，单位 nat。"""
    n, p = X.shape
    D = np.argsort(np.argsort(X, axis=0), axis=0) * bins // n   # 每列按秩切成 bins 个等频箱
    M = np.zeros((p, p))
    for i in range(p):
        for j in range(i + 1, p):
            joint = np.zeros((bins, bins)); np.add.at(joint, (D[:, i], D[:, j]), 1); joint /= n
            pi, pj = joint.sum(1), joint.sum(0); nz = joint > 0
            M[i, j] = M[j, i] = np.sum(joint[nz] * np.log(joint[nz] / np.outer(pi, pj)[nz]))
    return M

def clr(M):
    """CLR：每个 MI 相对于两个基因各自 MI 背景分布的 z 分数，负值截为 0，合并为 sqrt(zi^2 + zj^2)。"""
    p = M.shape[0]; off = ~np.eye(p, dtype=bool)
    mu = np.array([M[i, off[i]].mean() for i in range(p)]); sd = np.array([M[i, off[i]].std() for i in range(p)])
    Z = np.maximum(0, (M - mu[:, None]) / sd[:, None])
    S = np.sqrt(Z**2 + Z.T**2); np.fill_diagonal(S, 0); return S

def aracne_dpi(M, eps=0.0):
    """对每个三角形，若某边 MI 小于另两边的最小值（乘 1−eps），标记删除；最后一次性删除。"""
    p = M.shape[0]; keep = M.copy(); drop = np.zeros_like(M, dtype=bool)
    for i in range(p):
        for j in range(i + 1, p):
            for k in range(p):
                if k != i and k != j and M[i, j] < (1 - eps) * min(M[i, k], M[j, k]):
                    drop[i, j] = drop[j, i] = True; break
    keep[drop] = 0; return keep

def evaluate(S, gold, p):
    """S：对称打分矩阵；gold：真边集合 (i<j)。返回 AUROC、AUPR、随机 AUPR、早期精度比 EPR。"""
    iu = np.triu_indices(p, 1); score = S[iu]
    y = np.array([(i, j) in gold for i, j in zip(*iu)])
    order = np.argsort(-score, kind="stable"); ys = y[order]
    tp = np.cumsum(ys); prec = tp / np.arange(1, len(ys) + 1)
    aupr = np.sum(prec[ys]) / ys.sum()                      # 平均精确率：不做线性插值
    ranks = np.argsort(np.argsort(score)) + 1               # 秩和法计算 AUROC（并列忽略）
    auroc = (ranks[y].sum() - y.sum() * (y.sum() + 1) / 2) / (y.sum() * (~y).sum())
    k = y.sum(); base = y.mean()
    return auroc, aupr, base, (tp[k - 1] / k) / base

X, gold, kinds = simulate_dag()
p = X.shape[1]
R = np.abs(np.corrcoef(X.T)); np.fill_diagonal(R, 0)
MI = mutual_information(X)
print(f"真边 {len(gold)} 条 / 候选 {p*(p-1)//2} 条；其中 U 形调控 {sum(v=='u-shape' for v in kinds.values())} 条")
methods = [("|Pearson r|", R), ("互信息 MI", MI), ("CLR", clr(MI)), ("ARACNe(DPI)", aracne_dpi(MI))]
for name, S in methods:
    auroc, aupr, base, epr = evaluate(S, gold, p)
    print(f"{name:12s} AUROC = {auroc:.3f}  AUPR = {aupr:.3f}（随机 {base:.3f}，比值 {aupr/base:4.1f}）  EPR = {epr:4.1f}")
u = [e for e, v in kinds.items() if v == "u-shape"]
print("U 形调控边上的平均 |r| =", np.round(np.mean([R[e] for e in u]), 3), "；平均 MI =", np.round(np.mean([MI[e] for e in u]), 3),
      "；全部候选对的平均 MI =", np.round(MI[np.triu_indices(p, 1)].mean(), 3))

import matplotlib.pyplot as plt
iu = np.triu_indices(p, 1); y = np.array([(i, j) in gold for i, j in zip(*iu)])
fig, ax = plt.subplots(figsize=(6.4, 4.2))
for (name, S), lab, ls in zip(methods, ["|Pearson r|", "MI", "CLR", "ARACNe (DPI)"], ["-", "--", "-.", ":"]):
    ys = y[np.argsort(-S[iu], kind="stable")]; tp = np.cumsum(ys)
    ax.plot(tp / ys.sum(), tp / np.arange(1, len(ys) + 1), ls, lw=2,
            label=f"{lab} (AUPR = {evaluate(S, gold, p)[1]:.2f})")
ax.axhline(y.mean(), color="0.6", lw=0.8, label="random")
ax.set_xlabel("recall"); ax.set_ylabel("precision"); ax.set_ylim(0, 1.02)
ax.set_title("Precision–recall against the true regulatory edges"); ax.legend()
plt.tight_layout(); plt.show()