#!/usr/bin/env python3 """预注册判据的正式统计 + 效应量 + 配对补充分析。 08_analyze.py 已按判据出过一遍数；这一份补三件 08 没给全的东西： 1. claims.json 的 H2 原文写的是「Mann-Whitney U 双侧 p<0.05」，08 跑的是单侧(less)。 以 claims.json 为准，双侧作为判定用的 p，单侧一并列出。 2. 效应量：odds ratio 的 95% 置信区间（Woolf 对数法，含 Haldane 修正）、 Mann-Whitney 的 rank-biserial 相关、K 中位数差。 3. 配对补充：候选与背景是按比对长度 1:1 配对的，配对内的 K 差可以用符号秩检验， 比不配对的 MWU 更贴近本轮设计。非预注册，只作补充，不参与判定。 """ import json, math, statistics, sys from pathlib import Path from scipy.stats import fisher_exact, mannwhitneyu, wilcoxon def or_ci(a, b, c, d, z=1.959963984540054): """[[a,b],[c,d]] 的 odds ratio 与 95% CI（Woolf 对数法，全格 +0.5 Haldane 修正）。""" A, B, C, D = a + 0.5, b + 0.5, c + 0.5, d + 0.5 lor = math.log((A * D) / (B * C)) se = math.sqrt(1 / A + 1 / B + 1 / C + 1 / D) return math.exp(lor), (math.exp(lor - z * se), math.exp(lor + z * se)) def main(pkg): pkg = Path(pkg); res = pkg / "work" / "results" rows = [json.loads(l) for l in open(res / "relax_raw.jsonl") if l.strip()] rows = [r for r in rows if r.get("ok") and r.get("K") is not None] idx = {(r["gene"], r["set"], r["branch"]): r for r in rows} pairs = json.loads((res / "pairs.json").read_text())["pairs"] complete = [p for p in pairs if all( k in idx for k in ((p["candidate"], "candidate", "cave"), (p["candidate"], "candidate", "surface"), (p["background"], "background", "cave"), (p["background"], "background", "surface")))] out = {"n_pairs_planned": len(pairs), "n_pairs_complete": len(complete), "n_relax_runs_ok": len(rows), "note": "只用四格全部跑完的成对基因，保证 2x2 严格配对。"} def K(setname, branch): key = "candidate" if setname == "candidate" else "background" return [idx[(p[key], setname, branch)]["K"] for p in complete] def relaxed(setname, branch): key = "candidate" if setname == "candidate" else "background" rs = [idx[(p[key], setname, branch)] for p in complete] return sum(1 for r in rs if r["K"] < 1 and r["p_value"] < 0.05) def intensified(setname, branch): key = "candidate" if setname == "candidate" else "background" rs = [idx[(p[key], setname, branch)] for p in complete] return sum(1 for r in rs if r["K"] > 1 and r["p_value"] < 0.05) n = len(complete) for branch in ("cave", "surface"): a, c = relaxed("candidate", branch), relaxed("background", branch) b, d = n - a, n - c or_pt, (lo, hi) = or_ci(a, b, c, d) or_raw, p_two = fisher_exact([[a, b], [c, d]], alternative="two-sided") _, p_greater = fisher_exact([[a, b], [c, d]], alternative="greater") kc, kb = K("candidate", branch), K("background", branch) u_two, p_mw_two = mannwhitneyu(kc, kb, alternative="two-sided") u_less, p_mw_less = mannwhitneyu(kc, kb, alternative="less") # rank-biserial：U / (n1*n2) 换算，正值 = 候选 K 更小 rb = 1 - 2 * u_two / (len(kc) * len(kb)) out[branch] = { "counts": {"n_candidate": n, "n_background": n, "n_relaxed_candidate": a, "n_relaxed_background": c, "n_intensified_candidate": intensified("candidate", branch), "n_intensified_background": intensified("background", branch), "frac_relaxed_candidate": round(a / n, 4), "frac_relaxed_background": round(c / n, 4), "table_[[cand_rel,cand_not],[bg_rel,bg_not]]": [[a, b], [c, d]]}, "fisher": {"odds_ratio_sample": round(float(or_raw), 4), "odds_ratio_haldane": round(or_pt, 4), "odds_ratio_95CI": [round(lo, 4), round(hi, 4)], "p_two_sided": round(float(p_two), 6), "p_one_sided_greater": round(float(p_greater), 6)}, "K_distribution": { "median_candidate": round(statistics.median(kc), 4), "median_background": round(statistics.median(kb), 4), "median_difference_cand_minus_bg": round(statistics.median(kc) - statistics.median(kb), 4), "mean_candidate": round(statistics.mean(kc), 4), "mean_background": round(statistics.mean(kb), 4), "mannwhitney_U": float(u_two), "p_two_sided": round(float(p_mw_two), 6), "p_one_sided_candidate_lower": round(float(p_mw_less), 6), "rank_biserial_positive_means_candidate_lower": round(rb, 4)}, } # 预注册判定（严格照 claims.json 的 prediction / decision_rule） cav, sur = out["cave"], out["surface"] h1_ok = cav["fisher"]["p_one_sided_greater"] < 0.05 and cav["fisher"]["odds_ratio_sample"] > 1 h2_ok = (cav["K_distribution"]["p_two_sided"] < 0.05 and cav["K_distribution"]["median_candidate"] < cav["K_distribution"]["median_background"]) h3_refuted = (sur["fisher"]["p_one_sided_greater"] < 0.05 and sur["fisher"]["odds_ratio_sample"] > 1) out["prereg_verdict"] = { "H1": "成立" if h1_ok else "被推翻", "H1_basis": f"Fisher 单侧 p={cav['fisher']['p_one_sided_greater']}（判据 <0.05），" f"OR={cav['fisher']['odds_ratio_sample']}（判据 >1）；双侧 p={cav['fisher']['p_two_sided']}", "H2": "成立" if h2_ok else "被推翻", "H2_basis": f"Mann-Whitney 双侧 p={cav['K_distribution']['p_two_sided']}（claims.json 原文判据），" f"候选 K 中位数 {cav['K_distribution']['median_candidate']} < 背景 " f"{cav['K_distribution']['median_background']}；单侧 p=" f"{cav['K_distribution']['p_one_sided_candidate_lower']}", "H3_negative_control": "被推翻" if h3_refuted else "成立", "H3_basis": f"地表分支 Fisher 单侧 p={sur['fisher']['p_one_sided_greater']}，" f"OR={sur['fisher']['odds_ratio_sample']}（<1，方向与洞穴分支相反）", } # 配对补充（非预注册） supp = {} for branch in ("cave", "surface"): diffs = [x - y for x, y in zip(K("candidate", branch), K("background", branch))] stat, p = wilcoxon(diffs, alternative="two-sided") supp[f"paired_wilcoxon_{branch}"] = { "median_paired_diff_cand_minus_bg": round(statistics.median(diffs), 4), "n_pairs_candidate_lower": sum(1 for d in diffs if d < 0), "W": float(stat), "p_two_sided": round(float(p), 6)} # 同一批基因在两条分支上的 K 差（尺子的分支特异性；两个基因集分别看） for s in ("candidate", "background"): diffs = [x - y for x, y in zip(K(s, "cave"), K(s, "surface"))] stat, p = wilcoxon(diffs, alternative="two-sided") supp[f"branch_paired_wilcoxon_{s}"] = { "median_K_cave": round(statistics.median(K(s, "cave")), 4), "median_K_surface": round(statistics.median(K(s, "surface")), 4), "median_paired_diff_cave_minus_surface": round(statistics.median(diffs), 4), "W": float(stat), "p_two_sided": round(float(p), 6)} supp["_note"] = ("以上均非预注册判据，不参与 H1/H2/H3 判定，只作补充。" "分支间比较受该分支的分支长度、组装/注释质量等谱系特异因素影响，" "候选与背景两个集合上同向同量级出现，说明它不是候选基因集特有的信号。") out["supplementary_not_prereg"] = supp (res / "final_stats.json").write_text(json.dumps(out, ensure_ascii=False, indent=1)) print(json.dumps(out, ensure_ascii=False, indent=1)) if __name__ == "__main__": main(*sys.argv[1:])