差异共表达——Fisher z 检验与模块水平的置换检验

# 两种条件各 50 个样本、40 个基因：基因 0–9 在条件 A 中共表达（共享一个潜在因子），在条件 B 中失去协同（“重布线”）；
# 基因 10–19 在两种条件中都共表达，但条件 B 中平均表达整体上调（差异表达而非差异共表达）；其余为噪声基因。
# 逐对用 Fisher z 检验 H0: ρ_A = ρ_B，BH 法控制 FDR；再对“模块平均相关之差”做样本标签置换检验。
# 输出解读：差异共表达对集中在基因 0–9 内部；t 检验找到的差异表达基因是 10–19，两类信号互不重叠。
import numpy as np
from scipy import stats

rng = np.random.default_rng(5)
n, p = 50, 40
def condition(rewired, shift):
    X = rng.normal(size=(n, p))
    f = rng.normal(size=n); g = rng.normal(size=n)
    if not rewired:
        X[:, :10] = 0.8 * f[:, None] + 0.6 * rng.normal(size=(n, 10))
    X[:, 10:20] = 0.8 * g[:, None] + 0.6 * rng.normal(size=(n, 10)) + shift
    return X
XA, XB = condition(False, 0.0), condition(True, 1.0)

RA, RB = np.corrcoef(XA.T), np.corrcoef(XB.T)
iu = np.triu_indices(p, 1)
z = (np.arctanh(RA[iu]) - np.arctanh(RB[iu])) / np.sqrt(1 / (n - 3) + 1 / (n - 3))
pv = 2 * stats.norm.sf(np.abs(z))
m = len(pv); o = np.argsort(pv); q = np.empty(m)
q[o] = np.minimum.accumulate((pv[o] * m / np.arange(1, m + 1))[::-1])[::-1]   # BH 校正
sig = q < 0.05
in_mod = (iu[0] < 10) & (iu[1] < 10)
print(f"FDR<0.05 的差异共表达基因对：{sig.sum()} 个，其中位于基因 0–9 内部的 {np.sum(sig & in_mod)} 个（该区域共 {in_mod.sum()} 对）")
t, pt = stats.ttest_ind(XA, XB)
print("差异表达（t 检验 p<0.001）的基因：", np.where(pt < 1e-3)[0].tolist())

def mean_r(X, idx):
    R = np.corrcoef(X[:, idx].T); return R[np.triu_indices(len(idx), 1)].mean()
for name, idx in [("基因 0–9", np.arange(10)), ("基因 10–19", np.arange(10, 20))]:
    obs = mean_r(XA, idx) - mean_r(XB, idx)
    pooled = np.vstack([XA, XB]); null = []
    for _ in range(2000):
        perm = rng.permutation(2 * n)
        null.append(mean_r(pooled[perm[:n]], idx) - mean_r(pooled[perm[n:]], idx))
    p_perm = (1 + np.sum(np.abs(null) >= abs(obs))) / (1 + len(null))
    print(f"{name}：平均相关 A = {mean_r(XA, idx):.2f}，B = {mean_r(XB, idx):.2f}，差 = {obs:.2f}，置换 p = {p_perm:.4f}")

import matplotlib.pyplot as plt
fig, ax = plt.subplots(1, 2, figsize=(10, 4.2), constrained_layout=True)
for a, Rm, title in zip(ax, [RA, RB], ["condition A", "condition B"]):
    im = a.imshow(Rm, cmap="RdBu_r", vmin=-1, vmax=1, interpolation="nearest")
    a.set_title(title); a.set_xlabel("gene"); a.set_ylabel("gene")
    for b in (9.5, 19.5):
        a.axhline(b, color="k", lw=0.6); a.axvline(b, color="k", lw=0.6)
fig.colorbar(im, ax=ax, shrink=0.85, label="Pearson r")
fig.suptitle("Genes 0–9 lose co-expression in B (rewired); genes 10–19 keep it")
plt.show()