#!/usr/bin/env python3 """按 claims.json 预注册的判据跑全部实验臂。不做任何调参。 产出 work/results/scores.npz（原始分数，供复算）与 work/results/analysis.json。 """ import json, os, sys, time import numpy as np import scipy.sparse as sp HERE = os.path.dirname(os.path.abspath(__file__)) RES = os.path.join(HERE, "..", "results") ALPHA = 0.3 # 预注册，全程不调 TOL = 1e-10 MAXIT = 200 N_NULL = 1000 # 每个种子值抽 1000 次 SEEDS = [20260912, 20260913, 20260914] K_FOLD = 5 TEMPORAL_CUTOFF = 2021 def log(msg): print(f"[{time.strftime('%H:%M:%S')}] {msg}", flush=True) def row_normalize(A): d = np.asarray(A.sum(axis=1)).ravel() d[d == 0] = 1.0 return sp.diags(1.0 / d) @ A def rwr_batch(P_T, seed_matrix): """seed_matrix: (n, k) 列为重启分布（每列和为 1）。返回 (n, k) 稳态分数。""" x = seed_matrix.copy() for it in range(MAXIT): x_new = (1 - ALPHA) * (P_T @ x) + ALPHA * seed_matrix delta = np.abs(x_new - x).sum(axis=0).max() x = x_new if delta < TOL: break return x, it + 1 def auroc(scores, labels): """标准 Mann-Whitney U / AUROC，带并列秩平均。scores/labels 同长。""" order = np.argsort(scores, kind="mergesort") s = scores[order] y = labels[order] n = len(s) ranks = np.empty(n, dtype=np.float64) i = 0 while i < n: j = i while j + 1 < n and s[j + 1] == s[i]: j += 1 ranks[i:j + 1] = (i + j) / 2.0 + 1.0 i = j + 1 n1 = y.sum() n0 = n - n1 if n1 == 0 or n0 == 0: return float("nan") return float((ranks[y == 1].sum() - n1 * (n1 + 1) / 2.0) / (n1 * n0)) def auprc(scores, labels): order = np.argsort(-scores, kind="mergesort") y = labels[order] tp = np.cumsum(y) prec = tp / np.arange(1, len(y) + 1) n1 = y.sum() return float((prec * y).sum() / n1) if n1 else float("nan") def degree_bins(deg, n_bins=20): """按度数排名分 n_bins 个等大小的桶，用于度数配对抽样。""" order = np.argsort(deg, kind="mergesort") bin_of = np.empty(len(deg), dtype=np.int32) edges = np.array_split(order, n_bins) members = [] for b, idxs in enumerate(edges): bin_of[idxs] = b members.append(idxs) return bin_of, members def sample_degree_matched(rng, seed_idx, bin_of, members, forbidden): """为每个种子基因，从它所在的度数桶里抽一个不在 forbidden 里的替身。""" out = [] used = set() for g in seed_idx: pool = members[bin_of[g]] for _ in range(200): cand = int(pool[rng.integers(len(pool))]) if cand not in forbidden and cand not in used: out.append(cand) used.add(cand) break else: out.append(int(pool[rng.integers(len(pool))])) return np.array(out, dtype=np.int64) def seed_vectors(n, seed_sets): M = np.zeros((n, len(seed_sets)), dtype=np.float64) for k, s in enumerate(seed_sets): M[np.asarray(s), k] = 1.0 / len(s) return M def run_null(P_T, n, seed_idx, bin_of, members, forbidden, labels, mask, n_null, rng, batch=100): """返回 n_null 个度数配对随机种子集的 AUROC。""" aucs = [] drawn = [] for start in range(0, n_null, batch): k = min(batch, n_null - start) sets = [sample_degree_matched(rng, seed_idx, bin_of, members, forbidden) for _ in range(k)] drawn.extend(sets) X, _ = rwr_batch(P_T, seed_vectors(n, sets)) for c in range(k): aucs.append(auroc(X[mask, c], labels)) return np.array(aucs), drawn def main(): A = sp.load_npz(os.path.join(RES, "string_adj.npz")).tocsr() nodes = json.load(open(os.path.join(RES, "network_nodes.json"), encoding="utf-8")) symbols = nodes["symbols"] ms = json.load(open(os.path.join(RES, "mapped_sets.json"), encoding="utf-8")) gsets = json.load(open(os.path.join(RES, "gene_sets.json"), encoding="utf-8")) n = A.shape[0] deg_w = np.asarray(A.sum(axis=1)).ravel() gout_idx = np.array(sorted(ms["gout"].values()), dtype=np.int64) cvd_idx = np.array(sorted(ms["cvd"].values()), dtype=np.int64) gout_set, cvd_set = set(gout_idx.tolist()), set(cvd_idx.tolist()) P_T = row_normalize(A).T.tocsr() bin_of, members = degree_bins(deg_w) # ---------- 评估宇宙：剔除全部痛风种子（泄漏控制） ---------- mask = np.ones(n, dtype=bool) mask[gout_idx] = False labels = np.array([1 if i in cvd_set else 0 for i in np.where(mask)[0]], dtype=np.int64) log(f"评估宇宙 {mask.sum()} 个基因，其中心血管阳性 {labels.sum()}") # ---------- 主实验臂 ---------- X, iters = rwr_batch(P_T, seed_vectors(n, [gout_idx])) rwr_full = X[:, 0] auc_rwr = auroc(rwr_full[mask], labels) ap_rwr = auprc(rwr_full[mask], labels) auc_deg = auroc(deg_w[mask], labels) ap_deg = auprc(deg_w[mask], labels) log(f"RWR 收敛 {iters} 次迭代; AUROC gout_rwr={auc_rwr:.4f} degree_only={auc_deg:.4f}") # ---------- 零模型 1：度数配对随机种子 ---------- null_main = {} for s in SEEDS: rng = np.random.default_rng(s) aucs, _ = run_null(P_T, n, gout_idx, bin_of, members, gout_set, labels, mask, N_NULL, rng) null_main[s] = aucs log(f"零模型(度数配对) seed={s}: mean={aucs.mean():.4f} p95={np.quantile(aucs,0.95):.4f} max={aucs.max():.4f}") null_all = np.concatenate([null_main[s] for s in SEEDS]) emp_p_main = float((1 + (null_all >= auc_rwr).sum()) / (len(null_all) + 1)) # ---------- 零模型 2：标签置换 ---------- perm_aucs = [] for s in SEEDS: rng = np.random.default_rng(s + 1) sc = rwr_full[mask] for _ in range(N_NULL): perm_aucs.append(auroc(sc, rng.permutation(labels))) perm_aucs = np.array(perm_aucs) log(f"标签置换零模型: mean={perm_aucs.mean():.4f} ci95=[{np.quantile(perm_aucs,0.025):.4f},{np.quantile(perm_aucs,0.975):.4f}]") # ---------- 配对五折（H2）+ 正对照（H0） ---------- folds_rows = [] pos_rows = [] for s in SEEDS: rng = np.random.default_rng(s + 2) perm = rng.permutation(gout_idx) folds = np.array_split(perm, K_FOLD) for f, held in enumerate(folds): train = np.array([g for g in gout_idx if g not in set(held.tolist())], dtype=np.int64) # H2：用 4/5 痛风种子预测心血管（评估宇宙固定为全体非痛风种子） Xf, _ = rwr_batch(P_T, seed_vectors(n, [train])) a_real = auroc(Xf[mask, 0], labels) rs = sample_degree_matched(rng, train, bin_of, members, gout_set) Xr, _ = rwr_batch(P_T, seed_vectors(n, [rs])) a_null = auroc(Xr[mask, 0], labels) a_deg = auc_deg folds_rows.append({"seed": s, "fold": f, "n_train": int(len(train)), "auroc_gout_rwr": a_real, "auroc_matched_null": a_null, "auroc_degree_only": a_deg}) # H0 正对照：用 4/5 痛风种子预测留出的 1/5 痛风基因 pmask = np.ones(n, dtype=bool) pmask[train] = False plabels = np.array([1 if i in set(held.tolist()) else 0 for i in np.where(pmask)[0]], dtype=np.int64) pa_real = auroc(Xf[pmask, 0], plabels) pa_null = auroc(Xr[pmask, 0], plabels) pa_deg = auroc(deg_w[pmask], plabels) pos_rows.append({"seed": s, "fold": f, "n_held": int(len(held)), "auroc_gout_rwr": pa_real, "auroc_matched_null": pa_null, "auroc_degree_only": pa_deg}) log(f"五折完成 seed={s}") # ---------- 时间留出（H3） ---------- gy = {g: gsets["gout_genes"][g]["first_year"] for g in ms["gout"]} cy = {g: gsets["cvd_genes"][g]["first_year"] for g in ms["cvd"]} t_seed_idx = np.array(sorted(ms["gout"][g] for g in ms["gout"] if gy[g] <= TEMPORAL_CUTOFF), dtype=np.int64) old_cvd = {ms["cvd"][g] for g in ms["cvd"] if cy[g] <= TEMPORAL_CUTOFF} new_cvd = {ms["cvd"][g] for g in ms["cvd"] if cy[g] > TEMPORAL_CUTOFF} tmask = np.ones(n, dtype=bool) tmask[gout_idx] = False for i in old_cvd: tmask[i] = False tlabels = np.array([1 if i in new_cvd else 0 for i in np.where(tmask)[0]], dtype=np.int64) log(f"时间留出：种子 {len(t_seed_idx)}，宇宙 {tmask.sum()}，新增心血管阳性 {tlabels.sum()}") Xt, _ = rwr_batch(P_T, seed_vectors(n, [t_seed_idx])) auc_t = auroc(Xt[tmask, 0], tlabels) auc_t_deg = auroc(deg_w[tmask], tlabels) null_t = {} for s in SEEDS: rng = np.random.default_rng(s + 3) aucs, _ = run_null(P_T, n, t_seed_idx, bin_of, members, gout_set, tlabels, tmask, N_NULL, rng) null_t[s] = aucs log(f"时间留出零模型 seed={s}: mean={aucs.mean():.4f} p95={np.quantile(aucs,0.95):.4f}") null_t_all = np.concatenate([null_t[s] for s in SEEDS]) emp_p_t = float((1 + (null_t_all >= auc_t).sum()) / (len(null_t_all) + 1)) log(f"时间留出 AUROC={auc_t:.4f} degree={auc_t_deg:.4f} emp_p={emp_p_t:.4g}") np.savez_compressed(os.path.join(RES, "scores.npz"), rwr_full=rwr_full, deg_w=deg_w, mask=mask, labels=labels, rwr_temporal=Xt[:, 0], tmask=tmask, tlabels=tlabels, gout_idx=gout_idx, cvd_idx=cvd_idx, t_seed_idx=t_seed_idx, null_main=null_all, null_temporal=null_t_all, perm_aucs=perm_aucs) out = { "params": {"alpha": ALPHA, "tol": TOL, "score_min": 700, "seeds": SEEDS, "n_null_per_seed": N_NULL, "k_fold": K_FOLD, "temporal_cutoff": TEMPORAL_CUTOFF}, "universe": {"n_eval": int(mask.sum()), "n_positive": int(labels.sum()), "n_gout_seeds": int(len(gout_idx)), "prevalence": float(labels.mean())}, "main": {"auroc_gout_rwr": auc_rwr, "auprc_gout_rwr": ap_rwr, "auroc_degree_only": auc_deg, "auprc_degree_only": ap_deg, "null_mean": float(null_all.mean()), "null_sd": float(null_all.std(ddof=1)), "null_p95": float(np.quantile(null_all, 0.95)), "null_max": float(null_all.max()), "n_null": int(len(null_all)), "empirical_p": emp_p_main, "null_by_seed": {str(s): {"mean": float(v.mean()), "p95": float(np.quantile(v, 0.95)), "max": float(v.max()), "emp_p": float((1 + (v >= auc_rwr).sum()) / (len(v) + 1))} for s, v in null_main.items()}}, "label_permutation": {"mean": float(perm_aucs.mean()), "ci95": [float(np.quantile(perm_aucs, 0.025)), float(np.quantile(perm_aucs, 0.975))], "n": int(len(perm_aucs))}, "paired_folds_h2": folds_rows, "positive_control_h0": pos_rows, "temporal_h3": {"n_seeds": int(len(t_seed_idx)), "n_eval": int(tmask.sum()), "n_positive": int(tlabels.sum()), "auroc_gout_rwr": auc_t, "auroc_degree_only": auc_t_deg, "null_mean": float(null_t_all.mean()), "null_sd": float(null_t_all.std(ddof=1)), "null_p95": float(np.quantile(null_t_all, 0.95)), "empirical_p": emp_p_t, "n_null": int(len(null_t_all)), "null_by_seed": {str(s): {"mean": float(v.mean()), "emp_p": float((1 + (v >= auc_t).sum()) / (len(v) + 1))} for s, v in null_t.items()}}, } json.dump(out, open(os.path.join(RES, "analysis.json"), "w"), ensure_ascii=False, indent=1) log("写出 analysis.json") if __name__ == "__main__": main()