#!/usr/bin/env python3
"""检查一个会系统性削弱本轮结论的偏倚：洞穴种里已经假基因化/丢失的视觉基因，
根本进不了 1:1 同源集，也就进不了 RELAX。被丢掉的恰恰是退化最彻底的那些基因。"""
import json, sys
from pathlib import Path


def main(S, pkg):
    S, pkg = Path(S), Path(pkg)
    cand = set(open(pkg / "work/code/candidate_symbols.txt").read().split())
    rosa = json.loads((S / "prep/GCF_024868665.1.meta.json").read_text())["protein_to_gene"]
    dala = json.loads((S / "prep/GCF_015846415.1.meta.json").read_text())["protein_to_gene"]
    barb = json.loads((S / "prep/GCF_947034865.1.meta.json").read_text())["protein_to_gene"]
    sets = {"Triplophysa_rosa_cave": set(rosa.values()),
            "Triplophysa_dalaica_surface": set(dala.values()),
            "Barbatula_barbatula_outgroup": set(barb.values())}
    ortho = {r["gene_symbol"] for r in json.loads((pkg / "work/results/orthologs.json").read_text())["orthologs"]}
    out = {"n_candidate_symbols_listed": len(cand)}
    for name, s in sets.items():
        out[name] = {"candidate_symbols_present_after_QC": len(cand & s),
                     "all_genes_after_QC": len(s)}
    out["candidate_in_1to1_all_six"] = len(cand & ortho)
    # 只在地表种/外群里有、洞穴种里没有的候选基因 —— 丢失或注释不到的嫌疑名单
    missing_in_rosa = sorted((cand & sets["Triplophysa_dalaica_surface"] &
                              sets["Barbatula_barbatula_outgroup"]) - sets["Triplophysa_rosa_cave"])
    missing_in_dalaica = sorted((cand & sets["Triplophysa_rosa_cave"] &
                                 sets["Barbatula_barbatula_outgroup"]) - sets["Triplophysa_dalaica_surface"])
    out["candidate_absent_in_cave_but_present_in_both_surface_and_outgroup"] = missing_in_rosa
    out["candidate_absent_in_surface_but_present_in_both_cave_and_outgroup"] = missing_in_dalaica
    out["asymmetry_note"] = ("前者比后者明显多，说明洞穴种确有候选基因丢失/未通过质控；"
                             "这些基因进不了 RELAX，方向上会让本轮低估松弛程度。"
                             "两者相当则说明主要是注释噪声，不构成方向性偏倚。")
    (pkg / "work/results/gene_presence_bias.json").write_text(json.dumps(out, ensure_ascii=False, indent=1))
    print(json.dumps({k: v for k, v in out.items() if k != "asymmetry_note"}, ensure_ascii=False, indent=1))


if __name__ == "__main__":
    main(*sys.argv[1:])