#!/usr/bin/env python # build_composite_score.py — composite snipe-vs-GT benchmark score across assay # combinations x k-configs. import os, sys from itertools import combinations import numpy as np import pandas as pd HERE = os.path.dirname(os.path.abspath(__file__)) SRC = os.path.join(HERE, "..", "comparison_subset.tsv") sys.path.insert(0, os.path.join(HERE, "..", "manuscript_main_figures")) from _readlevel_mapping import attach_readlevel_mapping # adds 'mapping_gt' # (key, snipe_col, gt_col) METRICS = [ ("seq_change", "basepair_change_rate", "bamqc_total_event_rate"), ("depth", "mean_depth_of_reference_coverage", "qualimap_mean_coverage_data"), ("mapping", "reference_mapping_rate", "mapping_gt"), ("breadth", "fraction_of_reference_covered_by_sample_(%)", "qualimap_covbreadth_1X"), ] METRIC_KEYS = [m[0] for m in METRICS] ASSAYS = ["WGS_SIM", "INDELS_SIM", "WXS_SIM", "METAGENOMIC_SIM", "RNA_SIM"] SHORT = {"WGS_SIM":"WGS", "INDELS_SIM":"INDELS", "WXS_SIM":"WXS", "METAGENOMIC_SIM":"META", "RNA_SIM":"RNA"} WEIGHTS = {"seq_change": 2.0, "depth": 1.0, "mapping": 1.0, "breadth": 1.0} K1S = [21, 31, 51] DEPTH_GUARD = 250.0 CV_THRESH = 0.05 MIN_N = 3 EPS = 1e-3 BREADTH_UNIT_THRESH = 1.5 # p95 below this => snipe breadth is a 0-1 fraction, rescale to % def geomean(vals, eps=EPS): """Geometric mean over finite values, eps-floored to avoid hard zeros. NaN if none finite.""" v = np.asarray([x for x in vals if np.isfinite(x)], dtype=float) if v.size == 0: return np.nan v = np.clip(v, eps, None) return float(np.exp(np.mean(np.log(v)))) def metric_accuracy(gt, snipe, cv_thresh=CV_THRESH, min_n=MIN_N): """Hybrid per-metric accuracy in [0,1] for one assay at one config. Returns dict(mape, r, informative, n, acc). acc=NaN if not computable.""" x = np.asarray(gt, dtype=float); y = np.asarray(snipe, dtype=float) m = np.isfinite(x) & np.isfinite(y) x, y = x[m], y[m] n = int(x.size) out = dict(mape=np.nan, r=np.nan, informative=False, n=n, acc=np.nan) if n < min_n: return out # too few valid samples: metric dropped (acc stays NaN), per spec nz = x != 0 if nz.any(): out["mape"] = float(np.median(np.abs(y[nz] - x[nz]) / np.abs(x[nz]))) acc_mape = (1.0 - min(out["mape"], 1.0)) if np.isfinite(out["mape"]) else np.nan informative = (n >= min_n and x.std() > 0 and y.std() > 0 and abs(x.mean()) > 0 and (x.std() / abs(x.mean())) > cv_thresh) if informative: out["r"] = float(np.corrcoef(x, y)[0, 1]) out["informative"] = True out["acc"] = np.nan if np.isnan(acc_mape) else 0.5 * (acc_mape + max(0.0, out["r"])) else: out["acc"] = acc_mape return out def assay_accuracy(acc_by_metric, weights=WEIGHTS): """Blend per-metric accuracies into (equal_mean, weighted_mean) for one assay. NaN metrics dropped; weights renormalized over available metrics.""" items = [(k, a) for k, a in acc_by_metric.items() if np.isfinite(a)] if not items: return np.nan, np.nan accs = np.array([a for _, a in items], dtype=float) w = np.array([weights[k] for k, _ in items], dtype=float) return float(accs.mean()), float((w * accs).sum() / w.sum()) def drop_orphan_indels(d): """INDELS_SIM carries 66 sample names = 33 documented + 33 orphaned. The documented run has '_haploid_' in the name and matches 01b-wgs-simulation_indels_experiment's samples.tsv (33 rows), its reads in 004/wgs_reads_long_indels/ (33 files) and its results/ dirs (33). The other 33 ('5X_hi_err_...' with no '_haploid_') are an earlier superseded run of the SAME 33 parameter combos: they left GT files + snipe qc rows behind but have no samples.tsv entry, no reads and no results dir. Keeping them double-counts every INDELS design, so they are dropped here. """ orphan = (d["experiment_name"] == "INDELS_SIM") & (~d["sample_base"].str.contains("_haploid_")) return d[~orphan] def load(path=SRC, restrict_k1=True): d = pd.read_csv(path, sep="\t") d = drop_orphan_indels(d) d.loc[d["mean_depth_of_reference_coverage"] > DEPTH_GUARD, "mean_depth_of_reference_coverage"] = np.nan s, g = d["fraction_of_reference_covered_by_sample_(%)"], d["qualimap_covbreadth_1X"] if s.dropna().quantile(0.95) <= BREADTH_UNIT_THRESH and g.dropna().quantile(0.95) > BREADTH_UNIT_THRESH: d["fraction_of_reference_covered_by_sample_(%)"] = s * 100.0 d = attach_readlevel_mapping(d) # adds 'mapping_gt' if restrict_k1: d = d[d["k1"].isin(K1S)] return d.copy() def build_atomic(df): """(assay, k1, k2_ext, scale, metric) -> mape, r, informative, n, acc.""" rows = [] keys = ["k1", "edgemer_extension_length", "scale"] for a in ASSAYS: da = df[df["experiment_name"] == a] for (k1, k2, sc), g in da.groupby(keys): for key, scol, gcol in METRICS: res = metric_accuracy(g[gcol], g[scol]) rows.append(dict(assay=a, k1=int(k1), k2_ext=int(k2), scale=int(sc), metric=key, **res)) return pd.DataFrame(rows) def all_combos(assays=ASSAYS): out = [] for r in range(1, len(assays) + 1): out.extend(combinations(assays, r)) return out def build_composite(atomic, combos=None): """(combo, k1, k2_ext, scale) -> 6 score columns. `combos` defaults to all 31 subsets of ASSAYS; pass an explicit list of assay-tuples to score other groupings.""" if combos is None: combos = all_combos() acc = {(r.k1, r.k2_ext, r.scale, r.assay, r.metric): r.acc for r in atomic.itertuples(index=False)} configs = (atomic[["k1", "k2_ext", "scale"]] .drop_duplicates().itertuples(index=False, name=None)) configs = list(configs) rows = [] for combo in combos: label = "+".join(SHORT[a] for a in combo) for (k1, k2, sc) in configs: per_metric = {} for key in METRIC_KEYS: per_metric[key] = geomean([acc.get((k1, k2, sc, a, key), np.nan) for a in combo]) eq_assay, w_assay = [], [] for a in combo: abm = {key: acc.get((k1, k2, sc, a, key), np.nan) for key in METRIC_KEYS} e, w = assay_accuracy(abm, WEIGHTS) eq_assay.append(e); w_assay.append(w) rows.append(dict( combo=label, n_assays=len(combo), k1=k1, k2_ext=k2, scale=sc, score_seq_change=per_metric["seq_change"], score_depth=per_metric["depth"], score_mapping=per_metric["mapping"], score_breadth=per_metric["breadth"], score_equal=geomean(eq_assay), score_seq_weighted=geomean(w_assay))) return pd.DataFrame(rows) SCORE_COLS = ["score_seq_change","score_depth","score_mapping","score_breadth", "score_equal","score_seq_weighted"] def best_config_per_combo(comp, rank_by="score_equal"): """Per combo: the argmax-config row (ranked by rank_by) with all scores, plus the config that separately maximizes score_seq_weighted (seqw_k1/k2_ext/scale).""" idx = comp.groupby("combo")[rank_by].idxmax() cols = ["combo", "n_assays", "k1", "k2_ext", "scale"] + SCORE_COLS out = comp.loc[idx, cols].sort_values(["n_assays", "combo"]).reset_index(drop=True) idxw = comp.groupby("combo")["score_seq_weighted"].idxmax() seqw = (comp.loc[idxw, ["combo", "k1", "k2_ext", "scale"]] .rename(columns={"k1": "seqw_k1", "k2_ext": "seqw_k2_ext", "scale": "seqw_scale"})) return out.merge(seqw, on="combo", how="left") def headline_scores(comp): """Per combo, best value of each score flavor (independently maximized).""" g = comp.groupby(["combo", "n_assays"])[SCORE_COLS].max().reset_index() return g.sort_values(["n_assays", "combo"]).reset_index(drop=True) def robustness(comp, rank_by="score_equal"): """Global best-compromise config (max mean rank_by across combos) + per-combo score there vs each combo's own best. Returns (global_best_config dict, table).""" mean_by_cfg = comp.groupby(["k1", "k2_ext", "scale"])[rank_by].mean() k1, k2, sc = mean_by_cfg.idxmax() gbc = {"k1": int(k1), "k2_ext": int(k2), "scale": int(sc)} at_cfg = comp[(comp.k1 == k1) & (comp.k2_ext == k2) & (comp.scale == sc)] at_cfg = at_cfg.set_index("combo")[rank_by] own = comp.groupby("combo")[rank_by].max() tab = pd.DataFrame({"combo": own.index, "score_at_global_best": at_cfg.reindex(own.index).values, "own_best_score": own.values}) tab["tuning_gain"] = tab["own_best_score"] - tab["score_at_global_best"] return gbc, tab.sort_values("tuning_gain", ascending=False).reset_index(drop=True) def main(): df = load() atomic = build_atomic(df) comp = build_composite(atomic) best = best_config_per_combo(comp) head = headline_scores(comp) gbc, rob = robustness(comp) atomic.to_csv(os.path.join(HERE, "accuracy_atomic.tsv"), sep="\t", index=False) comp.to_csv(os.path.join(HERE, "composite_scores.tsv"), sep="\t", index=False) best.to_csv(os.path.join(HERE, "best_config_per_combo.tsv"), sep="\t", index=False) head.to_csv(os.path.join(HERE, "headline_scores.tsv"), sep="\t", index=False) rob.to_csv(os.path.join(HERE, "robustness.tsv"), sep="\t", index=False) with open(os.path.join(HERE, "global_best_config.txt"), "w") as f: f.write(f"global_best_config (max mean score_equal across combos): {gbc}\n") print("wrote accuracy_atomic / composite_scores / best_config_per_combo / headline_scores / robustness") print("global_best_config:", gbc) if __name__ == "__main__": main()