#!/usr/bin/env python3 """交付后补充分析（非预注册，记入 AMENDMENTS.md 第 5 条）：时间留出臂（H3）换冻结网络快照重跑。 缺陷：原轮 H3 的标签按 2021-12-31 切了，但特征网络用的是 STRING v12.0（2023 年快照）， 含 2022 年后积累的互作知识（含从 2022 年后文献挖出的 textmining 证据），不是严格前瞻。 修复：把网络换成 STRING v11.5（文件时间戳 2021-06 / 2021-10，早于 2022）。**唯一改动 = 网络快照。** 不改的东西（全部直接复用 03_score.py / 02 / 04 的函数与常量，不复制参数）： combined_score >= 700、权重 score/1000、无向取 max、去自环；符号只做精确匹配 preferred_name； alpha=0.3、TOL=1e-10、MAXIT=200 幂迭代；种子 = 首次报道年份 <= 2021 的已映射痛风基因； 阳性 = 首次报道年份 > 2021 的已映射心血管基因；评估宇宙剔除全部已映射痛风基因 + 2021 年底前已知心血管基因； 零模型 = 20 个等大小加权度数桶内逐个配对抽替身、forbidden = 全部痛风基因、每种子 1000 次、 3 个种子 20260912/13/14，时间臂 rng = default_rng(seed + 3)，主臂 rng = default_rng(seed)（与 03/04 完全相同）； 经验 p = (1 + #{null >= 实测}) / (N + 1)；z = (实测 - null 均值) / null 标准差(ddof=1)（与 06 相同）。 用法（需 numpy + scipy；本机 /usr/bin/python3）： python3 10_temporal_frozen_v115.py --snapshot v12.0 --data-dir <含 string_info/links/links_detailed 的目录> 从原始文件重建 v12.0 网络，逐位核对原轮结果（证明本脚本与原流程等价），再补跑 v12.0 去 textmining 的时间臂 python3 10_temporal_frozen_v115.py --snapshot v11.5 [--data-dir ../data/string_v11.5] python3 10_temporal_frozen_v115.py --assemble python3 10_temporal_frozen_v115.py --seed-matched 种子数对照：v12.0 网络只用 v11.5 时间臂的同一组种子（需先跑完两个 --snapshot） python3 10_temporal_frozen_v115.py --assemble 汇总 -> work/results/temporal_frozen_v115.json（含纯标准库独立复算） """ import argparse, gzip, hashlib, importlib.util, json, os, sys, time import numpy as np import scipy.sparse as sp HERE = os.path.dirname(os.path.abspath(__file__)) PKG = os.path.abspath(os.path.join(HERE, "..", "..")) RES = os.path.join(HERE, "..", "results") OUT = os.path.join(RES, "temporal_frozen") spec = importlib.util.spec_from_file_location("score", os.path.join(HERE, "03_score.py")) S = importlib.util.module_from_spec(spec); spec.loader.exec_module(S) SCORE_MIN = 700 # 与 02/04 相同 SNAP = { "v12.0": { "default_dir": os.path.join(HERE, "..", "data"), "info": "string_info.txt.gz", "links": "string_links.txt.gz", "detailed": "string_links_detailed.txt.gz", # RAW_DATA_PROVENANCE.md 记录的原轮 SHA256，必须对上才算复现原轮 "expected_sha256": { "string_info.txt.gz": "144de4b0d98c6a7dfde6ddc2591cf88657f27b989eadff4f501450c3ed1f0f1c", "string_links.txt.gz": "3e22f32572211aa341d5b4bd08d30c32e693e294603202120936872f87719d4f", "string_links_detailed.txt.gz": "28a2ea3b22c964e691c0f2f934c5e3e397a4b0a12c833d8787fa02b7847b91b6", }, }, "v11.5": { "default_dir": os.path.join(HERE, "..", "data", "string_v11.5"), "info": "9606.protein.info.v11.5.txt.gz", "links": "9606.protein.links.v11.5.txt.gz", "detailed": "9606.protein.links.detailed.v11.5.txt.gz", "urls": { "9606.protein.info.v11.5.txt.gz": "https://stringdb-static.org/download/protein.info.v11.5/9606.protein.info.v11.5.txt.gz", "9606.protein.links.v11.5.txt.gz": "https://stringdb-static.org/download/protein.links.v11.5/9606.protein.links.v11.5.txt.gz", "9606.protein.links.detailed.v11.5.txt.gz": "https://stringdb-downloads.org/download/protein.links.detailed.v11.5/9606.protein.links.detailed.v11.5.txt.gz", }, "expected_sha256": None, }, } # 结论判定规则 —— 在看到任何 v11.5 结果之前写死 DECISION_RULE = { "primary_arm": "v11.5 全通道（combined_score）时间留出臂；这是对审稿意见（特征未按时间切）的直接修复", "H3_rule_from_claims_json": "经验 p < 0.05 且 AUROC 高于同宇宙 degree_only（claims.json decision_rule.H3_refuted_if 原文的反面）", "不成立": "v11.5 全通道时间臂不满足 H3 规则", "减弱": "v11.5 满足 H3 规则，但 AUROC 相对零模型均值的增益 < v12.0 原轮增益的一半", "仍然成立": "v11.5 满足 H3 规则，且增益 >= v12.0 原轮增益的一半", "secondary_textmining": "去 textmining 的 v11.5 时间臂单独按 H3 规则判 passed/REFUTED 并单列，不并入主结论；" "理由：v11.5 的 textmining 语料本身早于 2022，它不构成时间泄漏，而是『与标签共享文献池』的另一类混杂", } # 种子数混杂的归因规则 —— 在看到 v12.0/34 种子结果之前写死（v11.5/34 已知不显著） ATTRIBUTION_RULE = { "why": "v11.5 上能映射的时间种子 37→34、新阳性 334→339；『不显著』可能是 ① 网络里混进 2022 年后知识（泄漏）或 ② 少 3 个种子功效下降（假象）", "control": "在 v12.0 上只用与 v11.5 相同的 34 个种子重跑时间留出臂；评估宇宙 / 阳性 / 零模型口径与原轮完全相同", "leakage": "v12.0/34 仍满足 H3 规则（p<0.05 且胜过度数基线），而 v11.5/34 不满足 → 归因为网络快照泄漏", "unstable_underpowered": "v12.0/34 也不满足 H3 规则 → 原 p=0.0077 依赖那 3 个种子，时间留出结论不稳、功效不足", } def log(msg): print(f"[{time.strftime('%H:%M:%S')}] {msg}", flush=True) def sha256(path): h = hashlib.sha256() with open(path, "rb") as fh: for chunk in iter(lambda: fh.read(1 << 20), b""): h.update(chunk) return h.hexdigest() def read_info(path): ensp2sym = {} with gzip.open(path, "rt", encoding="utf-8") as fh: next(fh) for line in fh: p = line.split("\t") ensp2sym[p[0]] = p[1] return ensp2sym def to_adjacency(edges, ensp2sym): """与 02_build_network.py / 04_sensitivity_channels.py 逐行相同的建图逻辑。""" nodes = sorted({a for a, b, s in edges} | {b for a, b, s in edges}) idx = {p: i for i, p in enumerate(nodes)} n = len(nodes) rows = np.fromiter((idx[a] for a, b, s in edges), dtype=np.int32, count=len(edges)) cols = np.fromiter((idx[b] for a, b, s in edges), dtype=np.int32, count=len(edges)) vals = np.fromiter((s / 1000.0 for a, b, s in edges), dtype=np.float64, count=len(edges)) A = sp.coo_matrix((vals, (rows, cols)), shape=(n, n)).tocsr() A = A.maximum(A.T) A.setdiag(0); A.eliminate_zeros() symbols = [ensp2sym[p] for p in nodes] sym_index = {} for i, s in enumerate(symbols): sym_index.setdefault(s, i) return A, nodes, symbols, sym_index def build_combined(links_path, ensp2sym): """02_build_network.py：combined_score >= 700。""" edges, kept, total = [], 0, 0 with gzip.open(links_path, "rt", encoding="utf-8") as fh: next(fh) for line in fh: a, b, s = line.split() total += 1 s = int(s) if s < SCORE_MIN: continue kept += 1 edges.append((a, b, s)) A, nodes, symbols, sym_index = to_adjacency(edges, ensp2sym) return A, nodes, symbols, sym_index, {"edges_total_in_file": total, "edges_kept_directed_rows": kept} def build_no_textmining(detailed_path, ensp2sym): """04_sensitivity_channels.py：只留 experimental + database，按 STRING 合并公式合成，>= 700。""" edges = [] with gzip.open(detailed_path, "rt", encoding="utf-8") as fh: header = next(fh).split() ie, idb = header.index("experimental"), header.index("database") for line in fh: f = line.split() e, d = int(f[ie]), int(f[idb]) if e == 0 and d == 0: continue comb = 1.0 - (1.0 - e / 1000.0) * (1.0 - d / 1000.0) s = int(round(comb * 1000)) if s < SCORE_MIN: continue edges.append((f[0], f[1], s)) A, nodes, symbols, sym_index = to_adjacency(edges, ensp2sym) return A, nodes, symbols, sym_index def map_sets(sym_index, gene_list): """02 的 map_set：精确匹配，不做别名。""" return {g: sym_index[g] for g in gene_list if g in sym_index} def null_summary(null_by_seed, observed): allv = np.concatenate([null_by_seed[s] for s in S.SEEDS]) mean, sd = float(allv.mean()), float(allv.std(ddof=1)) return allv, { "null_mean": mean, "null_sd": sd, "null_p95": float(np.quantile(allv, 0.95)), "null_max": float(allv.max()), "n_null": int(len(allv)), "n_null_ge_observed": int((allv >= observed).sum()), "empirical_p": float((1 + (allv >= observed).sum()) / (len(allv) + 1)), "z_vs_null": float((observed - mean) / sd), "delta_vs_null_mean": float(observed - mean), "null_by_seed": {str(s): {"mean": float(v.mean()), "sd": float(v.std(ddof=1)), "emp_p": float((1 + (v >= observed).sum()) / (len(v) + 1))} for s, v in null_by_seed.items()}, } def verdict_h3(auc, auc_deg, emp_p): return "passed" if (emp_p < 0.05 and auc > auc_deg) else "REFUTED" def temporal_arm(A, symbols, gout_map, cvd_map, gsets, tag, seed_whitelist=None): """03_score.py 的 H3 段，逐行对齐。seed_whitelist 仅用于种子数对照：只把白名单内的早期痛风基因当种子， 评估宇宙、forbidden、阳性口径一律不变（仍剔除全部已映射痛风基因）。""" n = A.shape[0] deg_w = np.asarray(A.sum(axis=1)).ravel() gout_idx = np.array(sorted(gout_map.values()), dtype=np.int64) gout_set = set(gout_idx.tolist()) P_T = S.row_normalize(A).T.tocsr() bin_of, members = S.degree_bins(deg_w) gy = {g: gsets["gout_genes"][g]["first_year"] for g in gout_map} cy = {g: gsets["cvd_genes"][g]["first_year"] for g in cvd_map} t_seed_idx = np.array(sorted(gout_map[g] for g in gout_map if gy[g] <= S.TEMPORAL_CUTOFF and (seed_whitelist is None or g in seed_whitelist)), dtype=np.int64) old_cvd = {cvd_map[g] for g in cvd_map if cy[g] <= S.TEMPORAL_CUTOFF} new_cvd = {cvd_map[g] for g in cvd_map if cy[g] > S.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"[{tag}] 时间留出：种子 {len(t_seed_idx)}，宇宙 {tmask.sum()}，新增心血管阳性 {tlabels.sum()}") Xt, iters = S.rwr_batch(P_T, S.seed_vectors(n, [t_seed_idx])) auc_t = S.auroc(Xt[tmask, 0], tlabels) auc_t_deg = S.auroc(deg_w[tmask], tlabels) null_t = {} for s in S.SEEDS: rng = np.random.default_rng(s + 3) aucs, _ = S.run_null(P_T, n, t_seed_idx, bin_of, members, gout_set, tlabels, tmask, S.N_NULL, rng) null_t[s] = aucs log(f"[{tag}] 时间留出零模型 seed={s}: mean={aucs.mean():.4f} p95={np.quantile(aucs, 0.95):.4f}") null_all, ns = null_summary(null_t, auc_t) res = { "n_seeds": int(len(t_seed_idx)), "n_eval": int(tmask.sum()), "n_positive": int(tlabels.sum()), "n_gout_mapped_excluded": int(len(gout_idx)), "n_old_cvd_excluded": int(len(old_cvd)), "n_new_cvd_mapped": int(len(new_cvd)), "rwr_iterations": int(iters), "auroc_gout_rwr": auc_t, "auroc_degree_only": auc_t_deg, "auprc_gout_rwr": S.auprc(Xt[tmask, 0], tlabels), "auprc_degree_only": S.auprc(deg_w[tmask], tlabels), "prevalence": float(tlabels.mean()), **ns, "delta_vs_degree_only": float(auc_t - auc_t_deg), "beats_degree_only": bool(auc_t > auc_t_deg), "verdict_H3_rule": verdict_h3(auc_t, auc_t_deg, ns["empirical_p"]), "seed_symbols": sorted(symbols[i] for i in t_seed_idx), "positive_symbols_in_universe": sorted(symbols[i] for i in np.where(tmask)[0][tlabels == 1]), } arrays = {"score": Xt[:, 0], "deg_w": deg_w, "tmask": tmask, "tlabels": tlabels, "null": null_all} return res, arrays def main_arm(A, gout_idx, cvd_set, tag): """03_score.py 主臂段 / 04_sensitivity_channels.py，逐行对齐（rng = default_rng(seed)）。""" n = A.shape[0] deg_w = np.asarray(A.sum(axis=1)).ravel() P_T = S.row_normalize(A).T.tocsr() bin_of, members = S.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) X, iters = S.rwr_batch(P_T, S.seed_vectors(n, [gout_idx])) auc = S.auroc(X[mask, 0], labels) auc_deg = S.auroc(deg_w[mask], labels) nulls = {} for s in S.SEEDS: rng = np.random.default_rng(s) a, _ = S.run_null(P_T, n, gout_idx, bin_of, members, set(gout_idx.tolist()), labels, mask, S.N_NULL, rng) nulls[s] = a log(f"[{tag}] 主臂零模型 seed={s}: mean={a.mean():.4f}") null_all, ns = null_summary(nulls, auc) log(f"[{tag}] 主臂 AUROC={auc:.4f} degree={auc_deg:.4f} emp_p={ns['empirical_p']:.4g}") return {"n_gout_seeds": int(len(gout_idx)), "n_eval": int(mask.sum()), "n_positive": int(labels.sum()), "rwr_iterations": int(iters), "auroc_gout_rwr": auc, "auroc_degree_only": auc_deg, "auprc_gout_rwr": S.auprc(X[mask, 0], labels), "auprc_degree_only": S.auprc(deg_w[mask], labels), **ns, "beats_degree_only": bool(auc > auc_deg), "verdict_H1_rule": "passed" if (ns["empirical_p"] < 0.05 and auc > auc_deg) else "REFUTED"}, null_all def gene_status(genes, info_names, sym_index): out = {"mapped": [], "symbol_absent_from_protein_info": [], "in_protein_info_but_no_edge_ge700": []} for g in sorted(genes): if g in sym_index: out["mapped"].append(g) elif g in info_names: out["in_protein_info_but_no_edge_ge700"].append(g) else: out["symbol_absent_from_protein_info"].append(g) return out def run_snapshot(snap, data_dir): cfg = SNAP[snap] data_dir = os.path.abspath(data_dir or cfg["default_dir"]) os.makedirs(OUT, exist_ok=True) gsets = json.load(open(os.path.join(RES, "gene_sets.json"), encoding="utf-8")) cells = {"snapshot": snap, "inputs": {}} # 不记录本机目录名（可能含需求方标识） for key in ("info", "links", "detailed"): fn = cfg[key] path = os.path.join(data_dir, fn) rec = {"file": fn, "bytes": os.path.getsize(path), "sha256": sha256(path)} if cfg.get("urls"): rec["url"] = cfg["urls"][fn] if cfg["expected_sha256"]: rec["matches_original_provenance_sha256"] = rec["sha256"] == cfg["expected_sha256"][fn] if not rec["matches_original_provenance_sha256"]: sys.exit(f"{fn} 的 SHA256 与 RAW_DATA_PROVENANCE.md 不符，不能当作原轮 v12.0 输入") cells["inputs"][key] = rec log(f"[{snap}] 输入文件哈希完成") ensp2sym = read_info(os.path.join(data_dir, cfg["info"])) info_names = set(ensp2sym.values()) # ---------- 全通道网络（02） ---------- A, nodes, symbols, sym_index, filemeta = build_combined(os.path.join(data_dir, cfg["links"]), ensp2sym) gout_map = map_sets(sym_index, gsets["gout_genes"]) cvd_map = map_sets(sym_index, gsets["cvd_genes"]) deg_c = np.diff(A.indptr) net = {"proteins_in_info": len(ensp2sym), **filemeta, "nodes": int(A.shape[0]), "undirected_edges": int(A.nnz // 2), "median_weighted_degree": float(np.median(np.asarray(A.sum(axis=1)).ravel())), "median_degree": int(np.median(deg_c)), "gout_mapped": len(gout_map), "gout_unmapped": len(gsets["gout_genes"]) - len(gout_map), "cvd_mapped": len(cvd_map), "cvd_unmapped": len(gsets["cvd_genes"]) - len(cvd_map), "overlap_mapped": len(set(gout_map) & set(cvd_map))} log(f"[{snap}] 全通道网络 {json.dumps(net, ensure_ascii=False)}") cells["network_combined"] = net cut = S.TEMPORAL_CUTOFF gout_early = [g for g, v in gsets["gout_genes"].items() if v["first_year"] <= cut] cvd_new = [g for g, v in gsets["cvd_genes"].items() if v["first_year"] > cut] cells["mapping_status"] = { "temporal_seed_candidates_raw_gwas": len(gout_early), "temporal_seed_status": gene_status(gout_early, info_names, sym_index), "new_cvd_candidates_raw_gwas": len(cvd_new), "new_cvd_status_counts": {k: len(v) for k, v in gene_status(cvd_new, info_names, sym_index).items()}, "new_cvd_status": gene_status(cvd_new, info_names, sym_index), } checks = {} if snap == "v12.0": A0 = sp.load_npz(os.path.join(RES, "string_adj.npz")).tocsr() nn0 = json.load(open(os.path.join(RES, "network_nodes.json"), encoding="utf-8")) ms0 = json.load(open(os.path.join(RES, "mapped_sets.json"), encoding="utf-8")) diff = (A - A0) if A.shape == A0.shape else None checks["adjacency_identical_to_string_adj_npz"] = bool(diff is not None and diff.nnz == 0) checks["node_symbols_identical"] = symbols == nn0["symbols"] checks["gout_mapping_identical"] = gout_map == ms0["gout"] checks["cvd_mapping_identical"] = cvd_map == ms0["cvd"] log(f"[{snap}] 建图复现核对 {checks}") else: sp.save_npz(os.path.join(OUT, f"adj_{snap}.npz"), A) json.dump({"symbols": symbols}, open(os.path.join(OUT, f"nodes_{snap}.json"), "w"), ensure_ascii=False) # ---------- 时间留出臂（主结论） ---------- t_res, t_arr = temporal_arm(A, symbols, gout_map, cvd_map, gsets, f"{snap}/全通道") cells["temporal_combined"] = t_res log(f"[{snap}] 时间留出 AUROC={t_res['auroc_gout_rwr']:.4f} degree={t_res['auroc_degree_only']:.4f} " f"null={t_res['null_mean']:.4f}±{t_res['null_sd']:.4f} emp_p={t_res['empirical_p']:.4g} z={t_res['z_vs_null']:.2f}") # ---------- 主臂（顺手：记录是否变化） ---------- gout_idx = np.array(sorted(gout_map.values()), dtype=np.int64) m_res, m_null = main_arm(A, gout_idx, set(cvd_map.values()), f"{snap}/全通道") cells["main_combined"] = m_res if snap == "v12.0": an = json.load(open(os.path.join(RES, "analysis.json"), encoding="utf-8")) z0 = np.load(os.path.join(RES, "scores.npz")) t0, mn0 = an["temporal_h3"], an["main"] checks["temporal_auroc_identical"] = t_res["auroc_gout_rwr"] == t0["auroc_gout_rwr"] checks["temporal_degree_auroc_identical"] = t_res["auroc_degree_only"] == t0["auroc_degree_only"] checks["temporal_null_draws_identical_3000"] = bool(np.array_equal(t_arr["null"], z0["null_temporal"])) checks["temporal_empirical_p_identical"] = t_res["empirical_p"] == t0["empirical_p"] checks["temporal_scores_identical"] = bool(np.array_equal(t_arr["score"], z0["rwr_temporal"])) checks["temporal_counts_identical"] = (t_res["n_seeds"], t_res["n_eval"], t_res["n_positive"]) == \ (t0["n_seeds"], t0["n_eval"], t0["n_positive"]) checks["main_auroc_identical"] = m_res["auroc_gout_rwr"] == mn0["auroc_gout_rwr"] checks["main_null_draws_identical_3000"] = bool(np.array_equal(m_null, z0["null_main"])) checks["main_empirical_p_identical"] = m_res["empirical_p"] == mn0["empirical_p"] # ---------- 去 textmining（04 的建图） ---------- A2, nodes2, symbols2, sym_index2 = build_no_textmining(os.path.join(data_dir, cfg["detailed"]), ensp2sym) # 与 04 相同：从全通道网络已映射的基因里，再查子网络 gout_map2 = {g: sym_index2[g] for g in gout_map if g in sym_index2} cvd_map2 = {g: sym_index2[g] for g in cvd_map if g in sym_index2} net2 = {"nodes": int(A2.shape[0]), "undirected_edges": int(A2.nnz // 2), "gout_seeds_in_subnetwork": len(gout_map2), "cvd_in_subnetwork": len(cvd_map2)} cells["network_no_textmining"] = net2 log(f"[{snap}] 去 textmining 网络 {net2}") t2_res, t2_arr = temporal_arm(A2, symbols2, gout_map2, cvd_map2, gsets, f"{snap}/去textmining") cells["temporal_no_textmining"] = t2_res log(f"[{snap}] 去textmining 时间留出 AUROC={t2_res['auroc_gout_rwr']:.4f} degree={t2_res['auroc_degree_only']:.4f} " f"null={t2_res['null_mean']:.4f}±{t2_res['null_sd']:.4f} emp_p={t2_res['empirical_p']:.4g}") gout_idx2 = np.array(sorted(gout_map2.values()), dtype=np.int64) m2_res, m2_null = main_arm(A2, gout_idx2, set(cvd_map2.values()), f"{snap}/去textmining") cells["main_no_textmining"] = m2_res if snap == "v12.0": sc0 = json.load(open(os.path.join(RES, "sensitivity_channels.json"), encoding="utf-8")) checks["no_textmining_network_size_identical"] = (net2["nodes"], net2["undirected_edges"]) == \ (sc0["nodes"], sc0["undirected_edges"]) checks["no_textmining_main_auroc_identical"] = m2_res["auroc_gout_rwr"] == sc0["auroc_gout_rwr"] checks["no_textmining_main_empirical_p_identical"] = m2_res["empirical_p"] == sc0["empirical_p"] checks["no_textmining_main_null_mean_identical"] = m2_res["null_mean"] == sc0["null_mean"] checks["all_pass"] = all(v for k, v in checks.items()) cells["reproduction_checks_vs_original_round"] = checks log(f"[{snap}] 复现核对 all_pass={checks['all_pass']} {checks}") np.savez_compressed(os.path.join(OUT, f"arrays_{snap}.npz"), temporal_score=t_arr["score"], temporal_deg_w=t_arr["deg_w"], temporal_tmask=t_arr["tmask"], temporal_tlabels=t_arr["tlabels"], temporal_null=t_arr["null"], temporal_notm_score=t2_arr["score"], temporal_notm_deg_w=t2_arr["deg_w"], temporal_notm_tmask=t2_arr["tmask"], temporal_notm_tlabels=t2_arr["tlabels"], temporal_notm_null=t2_arr["null"], main_null=m_null, main_notm_null=m2_null) if snap == "v11.5": # 纯文本导出，供纯标准库独立复算（与 08_dump_for_audit.py 同格式） tidx = np.where(t_arr["tmask"])[0] with open(os.path.join(OUT, "eval_scores_temporal_v115.tsv"), "w", encoding="utf-8") as fh: fh.write("gene\trwr_temporal_score\tweighted_degree\tnew_cvd_label\n") for pos, i in enumerate(tidx): fh.write("%s\t%r\t%r\t%d\n" % (symbols[i], float(t_arr["score"][i]), float(t_arr["deg_w"][i]), int(t_arr["tlabels"][pos]))) np.savetxt(os.path.join(OUT, "null_auroc_temporal_v115.txt"), t_arr["null"], fmt="%.17g") cells["generated_at_utc"] = time.strftime("%Y-%m-%dT%H:%M:%SZ", time.gmtime()) cells["env"] = {"python": sys.version.split()[0], "numpy": np.__version__, "scipy": sp.__name__ and __import__("scipy").__version__} json.dump(cells, open(os.path.join(OUT, f"cells_{snap}.json"), "w"), ensure_ascii=False, indent=1) log(f"[{snap}] 写出 cells_{snap}.json") def run_seed_matched(): """v12.0 网络（原轮 string_adj.npz，已在 --snapshot v12.0 中核对与原始文件重建逐位一致）， 种子限定为 v11.5 时间臂实际用到的种子符号。""" c115 = json.load(open(os.path.join(OUT, "cells_v11.5.json"), encoding="utf-8")) c12 = json.load(open(os.path.join(OUT, "cells_v12.0.json"), encoding="utf-8")) assert c12["reproduction_checks_vs_original_round"]["adjacency_identical_to_string_adj_npz"] whitelist = set(c115["temporal_combined"]["seed_symbols"]) gsets = json.load(open(os.path.join(RES, "gene_sets.json"), encoding="utf-8")) A = sp.load_npz(os.path.join(RES, "string_adj.npz")).tocsr() symbols = json.load(open(os.path.join(RES, "network_nodes.json"), encoding="utf-8"))["symbols"] ms = json.load(open(os.path.join(RES, "mapped_sets.json"), encoding="utf-8")) missing = sorted(whitelist - set(c12["temporal_combined"]["seed_symbols"])) res, arr = temporal_arm(A, symbols, ms["gout"], ms["cvd"], gsets, "v12.0/v11.5种子", seed_whitelist=whitelist) log(f"[v12.0/v11.5种子] 时间留出 AUROC={res['auroc_gout_rwr']:.4f} degree={res['auroc_degree_only']:.4f} " f"null={res['null_mean']:.4f}±{res['null_sd']:.4f} emp_p={res['empirical_p']:.4g} z={res['z_vs_null']:.2f}") np.savez_compressed(os.path.join(OUT, "arrays_v12.0_seedmatched.npz"), temporal_score=arr["score"], temporal_tmask=arr["tmask"], temporal_tlabels=arr["tlabels"], temporal_null=arr["null"]) out = {"what": "v12.0 网络 + v11.5 时间臂的同一组种子（种子数混杂对照）", "attribution_rule_written_before_this_run": ATTRIBUTION_RULE, "v11.5_seed_symbols_not_in_v12.0_seeds": missing, "v12.0_seeds_dropped": sorted(set(c12["temporal_combined"]["seed_symbols"]) - whitelist), "temporal": res, "generated_at_utc": time.strftime("%Y-%m-%dT%H:%M:%SZ", time.gmtime())} json.dump(out, open(os.path.join(OUT, "cells_v12.0_seedmatched.json"), "w"), ensure_ascii=False, indent=1) log("写出 cells_v12.0_seedmatched.json") # ------------------------------------------------------------------ 汇总 ------------------------------------------------------------------ def edge_pairs(A, symbols, restrict=None): U = sp.triu(A, k=1).tocoo() out = set() for r, c in zip(U.row, U.col): a, b = symbols[r], symbols[c] if restrict is not None and a not in restrict and b not in restrict: continue out.add((a, b) if a < b else (b, a)) return out def assemble(): from scipy import stats c12 = json.load(open(os.path.join(OUT, "cells_v12.0.json"), encoding="utf-8")) c115 = json.load(open(os.path.join(OUT, "cells_v11.5.json"), encoding="utf-8")) if not c12.get("reproduction_checks_vs_original_round", {}).get("all_pass"): sys.exit("v12.0 复现核对未全部通过，不能把差异归因于网络快照") csm = json.load(open(os.path.join(OUT, "cells_v12.0_seedmatched.json"), encoding="utf-8")) t12, t115, t12s = c12["temporal_combined"], c115["temporal_combined"], csm["temporal"] keys = ["n_seeds", "n_eval", "n_positive", "auroc_gout_rwr", "auroc_degree_only", "null_mean", "null_sd", "null_p95", "empirical_p", "z_vs_null", "delta_vs_null_mean", "delta_vs_degree_only", "auprc_gout_rwr", "auprc_degree_only", "prevalence", "verdict_H3_rule"] def table(a, b): rows = [] for k in keys: r = {"metric": k, "v12.0": a[k], "v11.5": b[k]} if isinstance(a[k], (int, float)) and not isinstance(a[k], bool): r["v11.5_minus_v12.0"] = b[k] - a[k] rows.append(r) return rows def table3(a, b, c): rows = [] for k in keys: r = {"metric": k, "v12.0_37seeds_original": a[k], "v12.0_34seeds_seedmatched": b[k], "v11.5_34seeds_frozen": c[k]} rows.append(r) return rows # ---- 种子数混杂归因（按跑前写死的 ATTRIBUTION_RULE） ---- if t115["verdict_H3_rule"] == "passed": attribution = "not_applicable_v11.5_passed" elif t12s["verdict_H3_rule"] == "passed": attribution = "leakage" else: attribution = "unstable_underpowered" # ---- 结论（按跑前写死的 DECISION_RULE） ---- if t115["verdict_H3_rule"] != "passed": verdict = "不成立" elif t115["delta_vs_null_mean"] < 0.5 * t12["delta_vs_null_mean"]: verdict = "减弱" else: verdict = "仍然成立" # ---- 网络差异诊断 ---- A12 = sp.load_npz(os.path.join(RES, "string_adj.npz")).tocsr() s12 = json.load(open(os.path.join(RES, "network_nodes.json"), encoding="utf-8"))["symbols"] A115 = sp.load_npz(os.path.join(OUT, "adj_v11.5.npz")).tocsr() s115 = json.load(open(os.path.join(OUT, "nodes_v11.5.json"), encoding="utf-8"))["symbols"] e12, e115 = edge_pairs(A12, s12), edge_pairs(A115, s115) seeds12, seeds115 = set(t12["seed_symbols"]), set(t115["seed_symbols"]) seedunion = seeds12 | seeds115 se12, se115 = edge_pairs(A12, s12, seedunion), edge_pairs(A115, s115, seedunion) pos12, pos115 = set(t12["positive_symbols_in_universe"]), set(t115["positive_symbols_in_universe"]) a12 = np.load(os.path.join(OUT, "arrays_v12.0.npz")) a115 = np.load(os.path.join(OUT, "arrays_v11.5.npz")) sc12 = {s12[i]: float(a12["temporal_score"][i]) for i in np.where(a12["temporal_tmask"])[0]} sc115 = {s115[i]: float(a115["temporal_score"][i]) for i in np.where(a115["temporal_tmask"])[0]} common = sorted(set(sc12) & set(sc115)) rho = stats.spearmanr([sc12[g] for g in common], [sc115[g] for g in common]).statistic # 公共宇宙 + 公共阳性上的对照（排除『只是阳性集合变了』的解释）：只用于诊断，不做判定 lab_common = np.array([1 if (g in pos12 and g in pos115) else 0 for g in common]) both_pos = int(lab_common.sum()) auc_common_12 = S.auroc(np.array([sc12[g] for g in common]), lab_common) auc_common_115 = S.auroc(np.array([sc115[g] for g in common]), lab_common) asm = np.load(os.path.join(OUT, "arrays_v12.0_seedmatched.npz")) sc12s = {s12[i]: float(asm["temporal_score"][i]) for i in np.where(asm["temporal_tmask"])[0]} auc_common_12s = S.auroc(np.array([sc12s[g] for g in common]), lab_common) network_comparison = { "combined_score_ge_700": { "v12.0": {k: c12["network_combined"][k] for k in ("proteins_in_info", "edges_total_in_file", "nodes", "undirected_edges", "median_degree", "gout_mapped", "cvd_mapped", "overlap_mapped")}, "v11.5": {k: c115["network_combined"][k] for k in ("proteins_in_info", "edges_total_in_file", "nodes", "undirected_edges", "median_degree", "gout_mapped", "cvd_mapped", "overlap_mapped")}, }, "no_textmining_subnetwork": {"v12.0": c12["network_no_textmining"], "v11.5": c115["network_no_textmining"]}, "edge_overlap_by_gene_symbol_pairs": { "v12.0_edges": len(e12), "v11.5_edges": len(e115), "shared": len(e12 & e115), "only_in_v12.0": len(e12 - e115), "only_in_v11.5": len(e115 - e12), "share_of_v12.0_edges_absent_in_v11.5": len(e12 - e115) / len(e12), }, "edges_touching_temporal_seeds": { "seed_union_size": len(seedunion), "v12.0": len(se12), "v11.5": len(se115), "shared": len(se12 & se115), "only_in_v12.0": len(se12 - se115), "only_in_v11.5": len(se115 - se12), }, "temporal_seeds": {"v12.0": len(seeds12), "v11.5": len(seeds115), "lost_in_v11.5": sorted(seeds12 - seeds115), "gained_in_v11.5": sorted(seeds115 - seeds12), "v11.5_mapping_status": c115["mapping_status"]["temporal_seed_status"]}, "new_cvd_positives_in_universe": {"v12.0": len(pos12), "v11.5": len(pos115), "shared": len(pos12 & pos115), "lost_in_v11.5": sorted(pos12 - pos115), "gained_in_v11.5": sorted(pos115 - pos12), "v11.5_new_cvd_mapping_status_counts": c115["mapping_status"]["new_cvd_status_counts"], "v11.5_new_cvd_symbol_absent_from_protein_info": c115["mapping_status"]["new_cvd_status"]["symbol_absent_from_protein_info"]}, "temporal_rwr_score_spearman_on_shared_universe": {"n_genes": len(common), "spearman_rho": float(rho)}, "diagnostic_auroc_on_shared_universe_and_shared_positives": { "note": "诊断用，非判定：只在两个版本都在宇宙里的基因上、只把两个版本都算阳性的基因当阳性，比较两份分数。用来区分『网络变了』与『种子/阳性集合变了』", "n_genes": len(common), "n_positives": both_pos, "auroc_v12.0_37seeds_scores": auc_common_12, "auroc_v12.0_34seeds_scores": auc_common_12s, "auroc_v11.5_34seeds_scores": auc_common_115}, } # ---- 纯标准库独立复算（复用 audit_delivery.py 的实现） ---- aspec = importlib.util.spec_from_file_location("audit", os.path.join(PKG, "audit_delivery.py")) AU = importlib.util.module_from_spec(aspec); aspec.loader.exec_module(AU) _, trows = AU.read_tsv(os.path.join(OUT, "eval_scores_temporal_v115.tsv")) pairs = [(float(r[1]), int(r[3])) for r in trows] deg_pairs = [(float(r[2]), int(r[3])) for r in trows] with open(os.path.join(OUT, "null_auroc_temporal_v115.txt"), encoding="utf-8") as fh: nulls = [float(x) for x in fh if x.strip()] a_std = AU.auroc(pairs) p_std = (1 + sum(1 for v in nulls if v >= a_std)) / (len(nulls) + 1.0) recheck = {"by": "audit_delivery.py 的纯标准库 auroc（不经 numpy/scipy），读 eval_scores_temporal_v115.tsv + null_auroc_temporal_v115.txt", "auroc_gout_rwr": a_std, "auroc_degree_only": AU.auroc(deg_pairs), "empirical_p": p_std, "n_rows": len(pairs), "n_null": len(nulls)} recheck["matches_numpy_path"] = (abs(a_std - t115["auroc_gout_rwr"]) < 1e-12 and abs(recheck["auroc_degree_only"] - t115["auroc_degree_only"]) < 1e-12 and abs(p_std - t115["empirical_p"]) < 1e-15) notm12, notm115 = c12["temporal_no_textmining"], c115["temporal_no_textmining"] strip = lambda d: {k: v for k, v in d.items() if k not in ("seed_symbols", "positive_symbols_in_universe")} out = { "_what": "交付后补充分析：原轮时间留出臂（H3）把 STRING v12.0 换成 2021 年发布的 STRING v11.5 重跑。唯一改动 = 网络快照；其余参数、种子、零模型、判据全部沿用 claims.json / 03_score.py。非预注册，见 AMENDMENTS.md 第 5 条。", "generated_at_utc": time.strftime("%Y-%m-%dT%H:%M:%SZ", time.gmtime()), "code": "work/code/10_temporal_frozen_v115.py（--snapshot v12.0 / --snapshot v11.5 / --assemble）", "env": c115["env"], "why": "审稿意见成立：原 H3 标签按 2021-12-31 切，特征网络 STRING v12.0 却含 2022 年后积累的互作与文献挖掘证据，不是严格前瞻。", "only_change": {"from": "STRING v12.0 (9606)", "to": "STRING v11.5 (9606)；三个下载文件 Last-Modified 为 2021-06-22（protein.info）与 2021-10-30（links / links.detailed），均早于 2022-01-01"}, "held_fixed": ["combined_score >= 700，权重 score/1000，无向取 max，去自环（02）", "符号只做精确匹配 preferred_name，不做别名（02）", "RWR alpha=0.3，幂迭代 L1 残差 < 1e-10 或 200 次（03）", "种子 = GWAS 首次报道年份 <= 2021 的已映射痛风基因；阳性 = 首次报道年份 >= 2022 的已映射心血管基因（03，AMENDMENTS 第 2 条）", "评估宇宙剔除全部已映射痛风基因 + 2021 年底前已知心血管基因（03）", "度数配对随机种子零模型：20 个等大小加权度数桶、逐个配对、forbidden=全部痛风基因、1000 次 × 3 种子 20260912/13/14、rng=default_rng(seed+3)（03）", "经验 p=(1+#{null>=实测})/(N+1)，z=(实测-null均值)/null标准差 ddof=1（03/06）", "基因集 work/results/gene_sets.json 原样复用（未重新抽取 GWAS）"], "inputs": {"v11.5": c115["inputs"], "v12.0": c12["inputs"]}, "reproduction_check_v12.0": {"what": "用同一脚本从原始 v12.0 文件重建网络并重跑，逐位对上原轮 analysis.json / scores.npz / sensitivity_channels.json 才允许把差异归因于快照", **c12["reproduction_checks_vs_original_round"]}, "decision_rule_written_before_v11.5_results": DECISION_RULE, "temporal_frozen_v11.5": strip(t115), "temporal_original_v12.0_reproduced": strip(t12), "temporal_v12.0_seedmatched_34": {"what": csm["what"], "v12.0_seeds_dropped": csm["v12.0_seeds_dropped"], "v11.5_seed_symbols_not_in_v12.0_seeds": csm["v11.5_seed_symbols_not_in_v12.0_seeds"], **strip(t12s)}, "seed_count_confound_control": {"attribution_rule_written_before_run": csm["attribution_rule_written_before_this_run"], "three_way": table3(t12, t12s, t115), "attribution": attribution}, "comparison_temporal_combined": table(t12, t115), "sensitivity_no_textmining_temporal": { "note": "原轮只在主臂上做过去 textmining（v12.0 主臂 0.5755→0.5431），时间臂从未做过；这里两个版本都补上", "v12.0": strip(notm12), "v11.5": strip(notm115), "comparison": table(notm12, notm115)}, "main_arm_not_temporal": { "note": "主臂不声称时间性，不在本次修复范围；顺手在 v11.5 上跑一遍，只记录不改结论", "combined": {"v12.0": c12["main_combined"], "v11.5": c115["main_combined"]}, "no_textmining": {"v12.0": c12["main_no_textmining"], "v11.5": c115["main_no_textmining"]}}, "network_comparison": network_comparison, "stdlib_recheck_v11.5_temporal": recheck, "conclusion": {"verdict": verdict, "seed_count_attribution": attribution, "v11.5_temporal_H3_rule": t115["verdict_H3_rule"], "v11.5_temporal_no_textmining_H3_rule": notm115["verdict_H3_rule"], "delta_vs_null_mean": {"v12.0": t12["delta_vs_null_mean"], "v11.5": t115["delta_vs_null_mean"]}}, } json.dump(out, open(os.path.join(RES, "temporal_frozen_v115.json"), "w"), ensure_ascii=False, indent=1) print(json.dumps({"verdict": verdict, "attribution": attribution, "v12.0_34": {k: t12s[k] for k in keys}, "v11.5": {k: t115[k] for k in keys}, "v12.0": {k: t12[k] for k in keys}, "notm_v12.0": {k: notm12[k] for k in keys}, "notm_v11.5": {k: notm115[k] for k in keys}, "main_v11.5": {k: c115["main_combined"][k] for k in ("auroc_gout_rwr", "auroc_degree_only", "null_mean", "empirical_p", "z_vs_null", "verdict_H1_rule")}, "main_notm_v11.5": {k: c115["main_no_textmining"][k] for k in ("auroc_gout_rwr", "auroc_degree_only", "null_mean", "empirical_p", "z_vs_null", "verdict_H1_rule")}, "network": {k: v for k, v in network_comparison.items() if k not in ("temporal_seeds", "new_cvd_positives_in_universe")}, "seeds": {k: v for k, v in network_comparison["temporal_seeds"].items() if k != "v11.5_mapping_status"}, "pos": {k: (v if not isinstance(v, list) or len(v) < 40 else f"{len(v)} genes") for k, v in network_comparison["new_cvd_positives_in_universe"].items()}, "recheck": recheck}, ensure_ascii=False, indent=1)) if __name__ == "__main__": ap = argparse.ArgumentParser() ap.add_argument("--snapshot", choices=list(SNAP)) ap.add_argument("--data-dir") ap.add_argument("--assemble", action="store_true") ap.add_argument("--seed-matched", action="store_true") a = ap.parse_args() if a.assemble: assemble() elif a.seed_matched: run_seed_matched() elif a.snapshot: run_snapshot(a.snapshot, a.data_dir) else: ap.error("需要 --snapshot 或 --assemble")