随附材料 · 实验代码
kg_experiment.py
"""
OA (Disease::DOID:8398) gene-candidate ranking robustness on Hetionet v1.0.
Only stdlib + numpy/scipy/sklearn. Deterministic given seed.
"""
import os, sys, gzip, json, time, argparse
os.environ.setdefault("OMP_NUM_THREADS", "2")
import numpy as np
import scipy.sparse as sp
from scipy.stats import kendalltau
RNG_GLOBAL = None
def log(msg):
print(f"[{time.strftime('%H:%M:%S')}] {msg}", flush=True)
def load_graph(nodes_path, edges_path):
ids, kinds = [], []
with open(nodes_path, "r", encoding="utf-8") as f:
header = f.readline()
for line in f:
parts = line.rstrip("\n").split("\t")
if len(parts) < 3:
continue
ids.append(parts[0])
kinds.append(parts[2])
id2idx = {nid: i for i, nid in enumerate(ids)}
N = len(ids)
kinds = np.array(kinds)
edges_by_meta = {}
with gzip.open(edges_path, "rt", encoding="utf-8") as f:
header = f.readline()
for line in f:
s, me, t = line.rstrip("\n").split("\t")
si, ti = id2idx[s], id2idx[t]
edges_by_meta.setdefault(me, []).append((si, ti))
return ids, kinds, id2idx, N, edges_by_meta
def build_meta_matrix(pairs, N):
if not pairs:
return sp.csr_matrix((N, N))
rows = np.array([p[0] for p in pairs])
cols = np.array([p[1] for p in pairs])
data = np.ones(len(pairs), dtype=np.float64)
M = sp.coo_matrix((data, (rows, cols)), shape=(N, N)).tocsr()
M.sum_duplicates()
M.data[:] = 1.0 # collapse multi-edges to single presence
return M
def build_full_homogeneous(edges_by_meta, N):
rows, cols = [], []
for me, pairs in edges_by_meta.items():
for s, t in pairs:
rows.append(s); cols.append(t)
rows.append(t); cols.append(s)
data = np.ones(len(rows), dtype=np.float64)
A = sp.coo_matrix((data, (rows, cols)), shape=(N, N)).tocsr()
A.sum_duplicates()
A.data[:] = 1.0
return A
def degree_weight(A, w):
"""Return D^-w * A * D^-w (elementwise edge reweighting by endpoint degree)."""
deg = np.asarray(A.sum(axis=1)).ravel()
deg_safe = np.where(deg > 0, deg, 1.0)
dinv = deg_safe ** (-w)
Dinv = sp.diags(dinv)
return (Dinv @ A @ Dinv).tocsr(), deg
EDGE_WEIGHT_EXP = {
"unweighted": 0.0,
"degree_normalized_w0.4": 0.4,
"inverse_degree": 1.0,
}
def ppr_scores(W_row_normalized, restart_idx, alpha, n_iter=30):
N = W_row_normalized.shape[0]
p = np.zeros(N)
p[restart_idx] = 1.0
restart = p.copy()
for _ in range(n_iter):
p = alpha * restart + (1 - alpha) * (W_row_normalized.T @ p)
return p
def row_normalize(A):
deg = np.asarray(A.sum(axis=1)).ravel()
deg_safe = np.where(deg > 0, deg, 1.0)
Dinv = sp.diags(1.0 / deg_safe)
return (Dinv @ A).tocsr()
def rank_of_positive(scores, pos_idx, exclude_idx_set):
"""1-indexed rank of pos_idx among all candidates (excluding exclude_idx_set, pos itself stays)."""
pos_score = scores[pos_idx]
mask = np.ones(len(scores), dtype=bool)
if exclude_idx_set:
mask[list(exclude_idx_set)] = False
mask[pos_idx] = True
cand_scores = scores[mask]
# rank = 1 + number of candidates strictly greater, ties get half-credit avg
greater = np.sum(cand_scores > pos_score)
equal = np.sum(cand_scores == pos_score) - 1 # subtract self
rank = 1 + greater + equal / 2.0
n_neg = mask.sum() - 1
return rank, n_neg
def auc_from_rank(rank, n_neg):
if n_neg <= 0:
return 1.0
return 1.0 - (rank - 1) / n_neg
def main():
ap = argparse.ArgumentParser()
ap.add_argument("--data-dir", default="data")
ap.add_argument("--out", default="results.json")
ap.add_argument("--seeds", nargs="+", type=int, default=[0, 1, 2])
ap.add_argument("--quick", action="store_true", help="tiny smoke-test run")
args = ap.parse_args()
t0 = time.time()
log("loading graph...")
ids, kinds, id2idx, N, edges_by_meta = load_graph(
os.path.join(args.data_dir, "nodes.tsv"),
os.path.join(args.data_dir, "edges.sif.gz"),
)
log(f"N={N} nodes, {sum(len(v) for v in edges_by_meta.values())} directed edge rows, "
f"{len(edges_by_meta)} metaedge types")
OA = "Disease::DOID:8398"
oa_idx = id2idx[OA]
gene_mask = (kinds == "Gene")
gene_indices = np.where(gene_mask)[0]
log(f"gene candidate pool size = {len(gene_indices)}")
DaG = build_meta_matrix(edges_by_meta.get("DaG", []), N)
DlA = build_meta_matrix(edges_by_meta.get("DlA", []), N)
AeG = build_meta_matrix(edges_by_meta.get("AeG", []), N)
GiG = build_meta_matrix(edges_by_meta.get("GiG", []), N)
GiG_sym = (GiG + GiG.T).tocsr()
GiG_sym.data[:] = 1.0
A_full = build_full_homogeneous(edges_by_meta, N)
log(f"A_full nnz={A_full.nnz}")
oa_gene_row = DaG.getrow(oa_idx)
positive_genes = oa_gene_row.indices.copy()
log(f"OA positive DaG genes = {len(positive_genes)} (expect 62)")
if args.quick:
positive_genes = positive_genes[:6]
log(f"[quick] truncated to {len(positive_genes)} positives")
other_disease_gene_idx = set(np.where(DaG.sum(axis=0) > 0)[1] if False else [])
# genes associated (DaG) with any OTHER disease -> hard-negative pool
DaG_csc = DaG.tocsc()
gene_disease_deg = np.asarray(DaG_csc.sum(axis=0)).ravel()
other_disease_mask = np.zeros(N, dtype=bool)
other_disease_mask[gene_indices] = gene_disease_deg[gene_indices] > 0
# remove OA's own positives from that "other disease" pool definition per-fold (handled later)
full_deg = np.asarray(A_full.sum(axis=1)).ravel()
# ---- Precompute weighted matrices per edge_weight scheme (reused across folds) ----
weighted_full = {}
weighted_dla = {}
weighted_aeg = {}
weighted_gig = {}
for scheme, w in EDGE_WEIGHT_EXP.items():
Wf, _ = degree_weight(A_full, w)
weighted_full[scheme] = Wf
Wa, _ = degree_weight(DlA, w)
weighted_dla[scheme] = Wa
Wb, _ = degree_weight(AeG, w)
weighted_aeg[scheme] = Wb
Wc, _ = degree_weight(GiG_sym, w)
weighted_gig[scheme] = Wc
log(f"weighted matrices built in {time.time()-t0:.1f}s")
# row-normalized transition matrices for PPR, per scheme
ppr_transition = {scheme: row_normalize(weighted_full[scheme]) for scheme in EDGE_WEIGHT_EXP}
# DWPC score vectors per scheme (fold-invariant: no DaG involved)
dwpc_scores = {}
for scheme in EDGE_WEIGHT_EXP:
MP1 = (weighted_dla[scheme] @ weighted_aeg[scheme]).tocsr()
MP2 = (weighted_dla[scheme] @ weighted_aeg[scheme] @ weighted_gig[scheme]).tocsr()
row = np.asarray(MP1.getrow(oa_idx).todense()).ravel() + np.asarray(MP2.getrow(oa_idx).todense()).ravel()
dwpc_scores[scheme] = row
log("dwpc score vectors ready")
# AA/CN score vectors per scheme, fold dependent only via zeroing oa-gene entries of weighted_full
def aa_scores(scheme, zero_pairs=()):
Wf = weighted_full[scheme]
row = Wf.getrow(oa_idx).copy().tolil()
for g in zero_pairs:
row[0, g] = 0.0
row = row.tocsr()
score = (row @ Wf)
return np.asarray(score.todense()).ravel()
def degree_null_scores(zero_gene=None):
d = full_deg.copy()
if zero_gene is not None:
d = d.copy()
d[zero_gene] -= 1
d[oa_idx] -= 1
return d
ALPHAS = [0.15, 0.30, 0.50]
SCHEMES = list(EDGE_WEIGHT_EXP.keys())
TOPKS = [20, 50, 100, 200]
NEG_SCHEMES = ["uniform", "degree_matched", "other_disease_genes"]
N_NEG_PER_FOLD = 20 if not args.quick else 5
scorer_names = ["scorer-ppr", "scorer-dwpc", "scorer-adamic-adar",
"baseline-degree-null"]
# ---- canonical config used for headline arm metrics (auroc/auprc/recall/mrr) ----
CANON_SCHEME = "degree_normalized_w0.4"
CANON_ALPHA = 0.15
def rank_gene_pool_scores(scores):
return scores[gene_indices]
def score_vector_ppr(scheme, alpha, zero_pair=None):
W = ppr_transition[scheme]
restart_idx = oa_idx
if zero_pair is None:
p = ppr_scores(W, restart_idx, alpha)
return p
gene = zero_pair
Wm = W.copy()
# zero row entries corresponding to edge oa<->gene in transition (approx: zero column gene in row oa's original out-mass already baked into row-normalization; we zero the raw weighted edge before row-normalizing for exactness)
Wf = weighted_full[scheme].tolil()
old_og = Wf[oa_idx, gene]
old_go = Wf[gene, oa_idx]
Wf[oa_idx, gene] = 0.0
Wf[gene, oa_idx] = 0.0
Wf = Wf.tocsr()
Wn = row_normalize(Wf)
p = ppr_scores(Wn, restart_idx, alpha)
return p
# -------------------- Part A: LOO evaluation for arms (canonical config, 3 seeds) --------------------
log("Part A: leave-one-out evaluation of arms (auroc/auprc/recall/mrr)")
arms_results = {name: {"per_seed": {}} for name in scorer_names + ["baseline-random-permutation"]}
fold_genes = positive_genes
per_fold_records = [] # jsonl-able raw records
g2l = {int(g): i for i, g in enumerate(gene_indices)} # global node idx -> gene-pool local idx
for seed in args.seeds:
rng = np.random.default_rng(seed)
fold_auc = {name: [] for name in scorer_names}
fold_recall = {name: {k: [] for k in TOPKS} for name in scorer_names}
fold_mrr = {name: [] for name in scorer_names}
fold_auprc_scores = {name: {ns: {"y": [], "s": []} for ns in NEG_SCHEMES} for name in scorer_names}
dm_negs = {} # gene -> degree-matched negatives, reused by the random-permutation arm
for gene in fold_genes:
exclude = set(fold_genes.tolist()) - {gene} # filtered ranking: exclude other true positives
# scorer: ppr canonical
p_ppr = score_vector_ppr(CANON_SCHEME, CANON_ALPHA, zero_pair=gene)
s_ppr = rank_gene_pool_scores(p_ppr)
rank_ppr, nneg_ppr = rank_of_positive(s_ppr, np.where(gene_indices == gene)[0][0],
{np.where(gene_indices == g2)[0][0] for g2 in exclude if g2 in set(gene_indices.tolist())})
auc_ppr = auc_from_rank(rank_ppr, nneg_ppr)
# scorer: dwpc canonical (fold invariant, but exclude-set differs -> recompute rank only)
s_dwpc_full = dwpc_scores[CANON_SCHEME]
s_dwpc = rank_gene_pool_scores(s_dwpc_full)
rank_dwpc, nneg_dwpc = rank_of_positive(s_dwpc, np.where(gene_indices == gene)[0][0],
{np.where(gene_indices == g2)[0][0] for g2 in exclude})
auc_dwpc = auc_from_rank(rank_dwpc, nneg_dwpc)
# scorer: adamic-adar canonical
s_aa_full = aa_scores(CANON_SCHEME, zero_pairs=[gene])
s_aa = rank_gene_pool_scores(s_aa_full)
rank_aa, nneg_aa = rank_of_positive(s_aa, np.where(gene_indices == gene)[0][0],
{np.where(gene_indices == g2)[0][0] for g2 in exclude})
auc_aa = auc_from_rank(rank_aa, nneg_aa)
# baseline: degree null
s_dn_full = degree_null_scores(zero_gene=gene)
s_dn = rank_gene_pool_scores(s_dn_full)
rank_dn, nneg_dn = rank_of_positive(s_dn, np.where(gene_indices == gene)[0][0],
{np.where(gene_indices == g2)[0][0] for g2 in exclude})
auc_dn = auc_from_rank(rank_dn, nneg_dn)
for name, rank, auc in [
("scorer-ppr", rank_ppr, auc_ppr),
("scorer-dwpc", rank_dwpc, auc_dwpc),
("scorer-adamic-adar", rank_aa, auc_aa),
("baseline-degree-null", rank_dn, auc_dn),
]:
fold_auc[name].append(auc)
fold_mrr[name].append(1.0 / rank)
for k in TOPKS:
fold_recall[name][k].append(1.0 if rank <= k else 0.0)
per_fold_records.append({
"seed": seed, "gene_idx": int(gene),
"rank_ppr": rank_ppr, "rank_dwpc": rank_dwpc, "rank_aa": rank_aa, "rank_degnull": rank_dn,
})
# negative sampling for AUPRC (per scorer, per neg scheme)
candidate_pool = np.array([g for g in gene_indices if g != gene and g not in exclude])
for ns in NEG_SCHEMES:
if ns == "uniform":
negs = rng.choice(candidate_pool, size=min(N_NEG_PER_FOLD, len(candidate_pool)), replace=False)
elif ns == "degree_matched":
pos_deg = full_deg[gene]
diffs = np.abs(full_deg[candidate_pool] - pos_deg)
order = np.argsort(diffs)
negs = candidate_pool[order[:min(N_NEG_PER_FOLD, len(candidate_pool))]]
dm_negs[int(gene)] = negs
else: # other_disease_genes
hard_pool = candidate_pool[other_disease_mask[candidate_pool]]
if len(hard_pool) < N_NEG_PER_FOLD:
extra = candidate_pool[~other_disease_mask[candidate_pool]]
hard_pool = np.concatenate([hard_pool, rng.choice(extra, size=N_NEG_PER_FOLD - len(hard_pool), replace=False)]) if len(extra) > 0 else hard_pool
negs = rng.choice(hard_pool, size=min(N_NEG_PER_FOLD, len(hard_pool)), replace=False)
for name, full_scores in [
("scorer-ppr", p_ppr), ("scorer-dwpc", s_dwpc_full),
("scorer-adamic-adar", s_aa_full), ("baseline-degree-null", s_dn_full),
]:
fold_auprc_scores[name][ns]["y"].append(1)
fold_auprc_scores[name][ns]["s"].append(full_scores[gene])
for ng in negs:
fold_auprc_scores[name][ns]["y"].append(0)
fold_auprc_scores[name][ns]["s"].append(full_scores[ng])
# random permutation baseline
rp_auc, rp_mrr, rp_recall = [], [], {k: [] for k in TOPKS}
rp_auprc_y, rp_auprc_s = [], []
for gene in fold_genes:
perm_scores = rng.permutation(len(gene_indices)).astype(float)
exclude = set(fold_genes.tolist()) - {gene}
excl_local = {np.where(gene_indices == g2)[0][0] for g2 in exclude}
gidx_local = np.where(gene_indices == gene)[0][0]
rank, nneg = rank_of_positive(perm_scores, gidx_local, excl_local)
rp_auc.append(auc_from_rank(rank, nneg))
rp_mrr.append(1.0 / rank)
for k in TOPKS:
rp_recall[k].append(1.0 if rank <= k else 0.0)
# same degree-matched negatives as the other arms, so AUPRC is comparable
rp_auprc_y.append(1)
rp_auprc_s.append(perm_scores[gidx_local])
for ng in dm_negs[int(gene)]:
rp_auprc_y.append(0)
rp_auprc_s.append(perm_scores[g2l[int(ng)]])
from sklearn.metrics import average_precision_score
arms_results.setdefault("baseline-random-permutation", {"per_seed": {}})
arms_results["baseline-random-permutation"]["per_seed"][seed] = {
"auroc": float(np.mean(rp_auc)), "mrr": float(np.mean(rp_mrr)),
**{f"recall_at_{k}": float(np.mean(rp_recall[k])) for k in TOPKS},
"auprc_degmatched": float(average_precision_score(rp_auprc_y, rp_auprc_s)),
}
for name in scorer_names:
y = fold_auprc_scores[name]["degree_matched"]["y"]
s = fold_auprc_scores[name]["degree_matched"]["s"]
auprc = average_precision_score(y, s) if len(set(y)) > 1 else float("nan")
arms_results[name]["per_seed"][seed] = {
"auroc": float(np.mean(fold_auc[name])),
"mrr": float(np.mean(fold_mrr[name])),
**{f"recall_at_{k}": float(np.mean(fold_recall[name][k])) for k in TOPKS},
"auprc_degmatched": float(auprc),
}
log(f"seed {seed} done: " + ", ".join(f"{n}:auroc={arms_results[n]['per_seed'][seed]['auroc']:.3f}" for n in scorer_names))
# -------------------- Part B: 36-setting grid for Jaccard/Kendall robustness (P1) + overlap w/ degree-null (P2) --------------------
log("Part B: 36-setting grid (full graph, no fold removal) for Jaccard/Kendall robustness")
GRID_SCORERS = ["ppr", "dwpc", "adamic_adar", "degree_null"]
grid_top50 = {}
grid_full_rank = {}
for scorer in GRID_SCORERS:
for scheme in SCHEMES:
for alpha in ALPHAS:
key = f"{scorer}|{scheme}|{alpha}"
if scorer == "ppr":
p = score_vector_ppr(scheme, alpha, zero_pair=None)
s = rank_gene_pool_scores(p)
elif scorer == "dwpc":
s = rank_gene_pool_scores(dwpc_scores[scheme])
elif scorer == "adamic_adar":
s = rank_gene_pool_scores(aa_scores(scheme))
else:
s = rank_gene_pool_scores(full_deg)
order = np.argsort(-s)
top50 = set(gene_indices[order[:50]].tolist())
grid_top50[key] = top50
grid_full_rank[key] = s
log(f"grid computed: {len(grid_top50)} settings")
keys = list(grid_top50.keys())
jacc_vals = []
kendall_vals = []
rng_sub = np.random.default_rng(0)
sub_idx = rng_sub.choice(len(gene_indices), size=min(1500, len(gene_indices)), replace=False)
for i in range(len(keys)):
for j in range(i + 1, len(keys)):
a, b = grid_top50[keys[i]], grid_top50[keys[j]]
inter = len(a & b); union = len(a | b)
jacc_vals.append(inter / union if union > 0 else 0.0)
ra = grid_full_rank[keys[i]][sub_idx]
rb = grid_full_rank[keys[j]][sub_idx]
tau, _ = kendalltau(ra, rb)
kendall_vals.append(float(tau) if tau == tau else 0.0)
jacc_vals = np.array(jacc_vals)
kendall_vals = np.array(kendall_vals)
log(f"pairwise settings compared = {len(jacc_vals)} (expect C(36,2)=630); "
f"jaccard median={np.median(jacc_vals):.3f}, kendall median={np.median(kendall_vals):.3f}")
# P2: degree-null canonical AUROC + top50 overlap vs canonical non-baseline arms
degnull_auroc_mean = float(np.mean([arms_results["baseline-degree-null"]["per_seed"][s]["auroc"] for s in args.seeds]))
degnull_top50_canon = grid_top50[f"degree_null|{CANON_SCHEME}|{CANON_ALPHA}"]
overlaps = {}
for scorer, gname in [("ppr", "scorer-ppr"), ("dwpc", "scorer-dwpc"), ("adamic_adar", "scorer-adamic-adar")]:
key = f"{scorer}|{CANON_SCHEME}|{CANON_ALPHA}"
ov = len(grid_top50[key] & degnull_top50_canon)
overlaps[gname] = ov
max_overlap_frac = max(overlaps.values()) / 50.0
log(f"degree-null canonical AUROC={degnull_auroc_mean:.3f}; top50 overlap with pipelines={overlaps} "
f"(max frac={max_overlap_frac:.2f})")
# -------------------- verdicts --------------------
p1_median_jaccard = float(np.median(jacc_vals))
p1_outcome = "confirmed" if p1_median_jaccard < 0.5 else "refuted"
p2_outcome = "confirmed" if (degnull_auroc_mean >= 0.70 and max_overlap_frac >= 0.30) else "refuted"
runtime = time.time() - t0
log(f"TOTAL runtime = {runtime:.1f}s")
# -------------------- assemble results.json --------------------
def arm_metrics(name):
seeds_data = arms_results[name]["per_seed"]
metric_keys = ["auroc", "mrr", "auprc_degmatched"] + [f"recall_at_{k}" for k in TOPKS]
out = {}
for mk in metric_keys:
vals = [seeds_data[s][mk] for s in args.seeds]
vals = [v for v in vals if v == v] # drop nan
out[mk] = {
"per_seed": [seeds_data[s][mk] for s in args.seeds],
"mean": float(np.mean(vals)) if vals else float("nan"),
"std": float(np.std(vals)) if vals else float("nan"),
}
return out
arms = [
{"name": "baseline-degree-null", "is_baseline": True,
"what": "忽略OA节点,仅用Gene在全图(homogeneous graph)中的度排序,走同样的62折留一评测",
"metrics": arm_metrics("baseline-degree-null")},
{"name": "baseline-random-permutation", "is_baseline": True,
"what": "每fold对候选基因随机打分,给出评测流水线下限",
"metrics": arm_metrics("baseline-random-permutation")},
{"name": "scorer-ppr", "is_baseline": False,
"what": f"以OA为重启源的个性化PageRank,行归一化转移矩阵,30轮幂迭代,canonical config: edge_weight={CANON_SCHEME}, alpha={CANON_ALPHA}",
"metrics": arm_metrics("scorer-ppr")},
{"name": "scorer-dwpc", "is_baseline": False,
"what": f"度加权路径计数,元路径 D-DlA-A-AeG-G (2跳) 与 D-DlA-A-AeG-G-GiG-G (3跳) 之和,canonical edge_weight={CANON_SCHEME}",
"metrics": arm_metrics("scorer-dwpc")},
{"name": "scorer-adamic-adar", "is_baseline": False,
"what": f"OA与Gene在同质化全图上的广义Adamic-Adar/共邻指标(度幂次由edge_weight方案决定),canonical edge_weight={CANON_SCHEME}",
"metrics": arm_metrics("scorer-adamic-adar")},
]
headline = {
"metric": "auroc",
"baseline_mean": arm_metrics("baseline-degree-null")["auroc"]["mean"],
# 必须对应 arms 里某一条,取 canonical PPR;三臂均值在 arms 里查不到出处
"proposed_mean": arm_metrics("scorer-ppr")["auroc"]["mean"],
"delta": None,
"claim": (f"36套设定两两Top-50 Jaccard中位数={p1_median_jaccard:.3f}"
f"({'<0.5,P1成立' if p1_outcome=='confirmed' else '>=0.5,P1被推翻,候选名单对设定选择不敏感'})；"
f"纯度数零模型AUROC={degnull_auroc_mean:.3f},与三条推理臂canonical Top-50最大重合率={max_overlap_frac:.2f}"
f"({'均达标,P2成立' if p2_outcome=='confirmed' else '未同时达标,P2被推翻'})。"),
}
headline["delta"] = headline["proposed_mean"] - headline["baseline_mean"]
negative_result = (p1_outcome == "refuted") or (p2_outcome == "refuted")
results = {
"task_slug": "oa-kg-gene-ranking-robustness",
"question": "在Hetionet上做OA候选基因排序时,候选名单对边权/打分函数/重启率/负样本采样等设定的敏感程度,以及纯度数零模型的水平",
"falsifiable_prediction": ("P1: 36套设定两两Top-50 Jaccard中位数<0.5。"
"P2: 纯度数零模型AUROC>=0.70且其Top-50与完整流程Top-50重合>=30%。"),
# schema 只收 confirmed|refuted|inconclusive;两条独立预测一成立一推翻,合并只能是混合
"prediction_outcome": (p1_outcome if p1_outcome == p2_outcome else "inconclusive"),
"prediction_outcome_detail": f"P1={p1_outcome}, P2={p2_outcome}",
"negative_result": negative_result,
"dataset": f"Hetionet v1.0 (CC0), {N} nodes, {A_full.nnz//2} undirected homogeneous edges, "
f"OA={OA} 有 {len(positive_genes)} 条DaG边(留一交叉验证), Gene候选池={len(gene_indices)}",
"env": {"python": "3.12.3", "torch": "n/a (no torch used)",
"extra_packages": []},
"seeds": list(args.seeds),
"arms": arms,
"headline": headline,
"sweep_robustness": {
"n_settings": len(keys),
"n_pairs": len(jacc_vals),
"jaccard_top50_median": p1_median_jaccard,
"jaccard_top50_p25": float(np.percentile(jacc_vals, 25)),
"jaccard_top50_p75": float(np.percentile(jacc_vals, 75)),
"kendall_tau_median": float(np.median(kendall_vals)),
"kendall_tau_p25": float(np.percentile(kendall_vals, 25)),
"kendall_tau_p75": float(np.percentile(kendall_vals, 75)),
"kendall_subsample_n_genes": int(len(sub_idx)),
},
"degree_null_vs_pipelines": {
"degnull_auroc_mean": degnull_auroc_mean,
"top50_overlap_with_pipelines": overlaps,
"max_overlap_fraction": max_overlap_frac,
},
"deviations": [
"DWPC只实现了2条代表性元路径 (D-DlA-A-AeG-G 与 D-DlA-A-AeG-G-GiG-G),未覆盖原计划提到的DaG/GiG/GpBP全部元路径组合,因40分钟预算不够枚举更多;两条元路径均不使用DaG边,天然避免与留一评测标签泄漏。",
"PPR/AA的edge_weight加权应用于全图同质邻接矩阵的度幂次归一化,不是逐metaedge类型分别加权;DWPC的edge_weight则是在各自metaedge的二部图度上归一化——三者共用同一个w取值(0/0.4/1.0)但数学形式不同,这是同一个'边权强度旋钮'在不同打分函数里的具体化,不是同一个矩阵。",
"PPR的edge_weight归一化度数用的是全图整体度(不含留一fold的单边移除),留一只体现在:1)转移矩阵里把该OA-gene边置零后重新行归一化,2)degree-null分数里从全图度数减1。度数归一化系数D^-w本身不随fold重算,是预算内的简化,已在NOTES.md说明。",
"Kendall tau为提速在1500个随机抽样基因上计算(不是全部20945个候选),抽样种子固定为0,已记录在sweep_robustness字段。",
"AUPRC(auprc_degmatched)每fold只采样20个负样本(--quick模式5个)做pooled average precision,不是对全量负样本做PR曲线。",
],
"runtime_sec": runtime,
}
with open(args.out, "w", encoding="utf-8") as f:
json.dump(results, f, ensure_ascii=False, indent=2)
log(f"results written to {args.out}")
with open("logs/per_fold_raw.jsonl", "w", encoding="utf-8") as f:
for rec in per_fold_records:
f.write(json.dumps(rec) + "\n")
log(f"per-fold raw ranks written to logs/per_fold_raw.jsonl ({len(per_fold_records)} records)")
with open("logs/grid_top50.json", "w", encoding="utf-8") as f:
json.dump({k: sorted(list(v)) for k, v in grid_top50.items()}, f)
log("grid top50 lists written to logs/grid_top50.json")
return results
if __name__ == "__main__":
main()
make_figs.py
import json
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
with open("results.json", "r", encoding="utf-8") as f:
R = json.load(f)
fig, axes = plt.subplots(1, 3, figsize=(15, 4.5))
# Panel 1: AUROC per arm across seeds
ax = axes[0]
names = [a["name"] for a in R["arms"]]
means = [a["metrics"]["auroc"]["mean"] for a in R["arms"]]
stds = [a["metrics"]["auroc"]["std"] for a in R["arms"]]
colors = ["#888" if a["is_baseline"] else "#2b6cb0" for a in R["arms"]]
ax.bar(range(len(names)), means, yerr=stds, color=colors, capsize=4)
ax.set_xticks(range(len(names)))
ax.set_xticklabels(names, rotation=30, ha="right", fontsize=8)
ax.set_ylabel("AUROC (62-fold LOO, mean over seeds)")
ax.axhline(0.5, color="k", linestyle="--", linewidth=0.8)
ax.set_title("Arm comparison: AUROC")
# Panel 2: Jaccard(top50) distribution across 630 setting-pairs
ax = axes[1]
sw = R["sweep_robustness"]
# reconstruct approx distribution shape using summary stats only (median/p25/p75) as a box-ish view;
# also load raw grid top50 to compute the real histogram
with open("logs/grid_top50.json", "r", encoding="utf-8") as f:
grid = json.load(f)
keys = list(grid.keys())
sets = {k: set(v) for k, v in grid.items()}
jvals = []
for i in range(len(keys)):
for j in range(i + 1, len(keys)):
a, b = sets[keys[i]], sets[keys[j]]
inter = len(a & b); union = len(a | b)
jvals.append(inter / union if union else 0.0)
ax.hist(jvals, bins=30, color="#2b6cb0", alpha=0.8)
ax.axvline(0.5, color="red", linestyle="--", label="P1 threshold 0.5")
ax.axvline(np.median(jvals), color="k", linestyle="-", label=f"median={np.median(jvals):.2f}")
ax.set_xlabel("Top-50 Jaccard similarity (630 setting pairs)")
ax.set_ylabel("count")
ax.set_title("Robustness across 36 settings")
ax.legend(fontsize=8)
# Panel 3: degree-null vs pipeline overlap
ax = axes[2]
ov = R["degree_null_vs_pipelines"]["top50_overlap_with_pipelines"]
names2 = list(ov.keys())
vals2 = [ov[n] for n in names2]
ax.bar(names2, vals2, color="#c05621")
ax.axhline(15, color="red", linestyle="--", label="P2 threshold (30% of 50)")
ax.set_ylabel("# genes shared with degree-null Top-50 (out of 50)")
ax.set_xticklabels(names2, rotation=20, ha="right", fontsize=8)
ax.set_title(f"Degree-null AUROC={R['degree_null_vs_pipelines']['degnull_auroc_mean']:.2f}")
ax.legend(fontsize=8)
plt.tight_layout()
plt.savefig("figs/summary.png", dpi=150)
print("wrote figs/summary.png")
run.sh
#!/usr/bin/env bash
set -euo pipefail
cd "$(dirname "$0")"
export OMP_NUM_THREADS=2
PY=~/joulebeat-evolve/outreach/env/bin/python
mkdir -p data figs logs
if [ ! -f data/edges.sif.gz ]; then
curl -sL -o data/edges.sif.gz "https://github.com/hetio/hetionet/raw/main/hetnet/tsv/hetionet-v1.0-edges.sif.gz"
fi
if [ ! -f data/nodes.tsv ]; then
curl -sL -o data/nodes.tsv "https://raw.githubusercontent.com/hetio/hetionet/main/hetnet/tsv/hetionet-v1.0-nodes.tsv"
fi
# 整轮 stdout/stderr 落盘,保证 logs/run.log 永远是当前这批数字的日志
{
"$PY" kg_experiment.py --seeds 0 1 2 --out results.json
"$PY" make_figs.py
} 2>&1 | tee logs/run.log