#!/usr/bin/env python3
"""敌意复核：把 STRING 的 textmining 通道整个拿掉，只留 experimental + database
两个实验性通道重建网络，重跑主分析。

为什么必须做：STRING 的 combined_score 里 textmining 通道是从文献里挖出来的，
而 GWAS Catalog 的心血管标签也来自同一批文献。不拆开的话，
「痛风邻域富集心血管基因」完全可能只是「这两个词经常出现在同一篇摘要里」。
"""
import gzip, json, os
import numpy as np
import scipy.sparse as sp
import importlib.util

HERE = os.path.dirname(os.path.abspath(__file__))
RES = os.path.join(HERE, "..", "results")
DATA = os.path.join(HERE, "..", "data")
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


def main():
    ensp2sym = {}
    with gzip.open(os.path.join(DATA, "string_info.txt.gz"), "rt", encoding="utf-8") as fh:
        next(fh)
        for line in fh:
            p = line.split("\t")
            ensp2sym[p[0]] = p[1]

    # protein.links.detailed 列: protein1 protein2 neighborhood fusion cooccurence
    #                            coexpression experimental database textmining combined_score
    edges = []
    with gzip.open(os.path.join(DATA, "string_links_detailed.txt.gz"), "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
            # STRING 官方合并公式（去掉先验后概率合并），只用这两个通道
            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))

    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)

    ms = json.load(open(os.path.join(RES, "mapped_sets.json"), encoding="utf-8"))
    gout_idx = np.array(sorted({sym_index[g] for g in ms["gout"] if g in sym_index}), dtype=np.int64)
    cvd_set = {sym_index[g] for g in ms["cvd"] if g in sym_index}

    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, _ = 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.append(a)
    null_all = np.concatenate(nulls)
    emp_p = float((1 + (null_all >= auc).sum()) / (len(null_all) + 1))

    out = {"network": "STRING v12.0 仅 experimental+database 通道，合并分 >=700，textmining 完全剔除",
           "nodes": n, "undirected_edges": int(A.nnz // 2),
           "gout_seeds_in_subnetwork": int(len(gout_idx)),
           "n_eval": int(mask.sum()), "n_positive": int(labels.sum()),
           "auroc_gout_rwr": auc, "auroc_degree_only": auc_deg,
           "null_mean": float(null_all.mean()), "null_sd": float(null_all.std(ddof=1)),
           "null_p95": float(np.quantile(null_all, 0.95)), "empirical_p": emp_p,
           "n_null": int(len(null_all))}
    json.dump(out, open(os.path.join(RES, "sensitivity_channels.json"), "w"), ensure_ascii=False, indent=1)
    print(json.dumps(out, ensure_ascii=False, indent=1))


if __name__ == "__main__":
    main()