""" miRNA threshold analysis — unified thresholds ---------------------------------------------- Round-number concentration thresholds applied to both tissues. Each threshold panel shows cumulative cohort samples passing (healthy / diseased), annotated with the approximate expected miRNA yield from pilot samples whose concentration falls in that bucket (>= threshold, < next threshold), combining both tissues. Outputs: summary Excel + figure """ import re import numpy as np import pandas as pd import matplotlib.pyplot as plt from matplotlib.patches import Patch # ── config ──────────────────────────────────────────────────────────────────── QC_FILE = "BEA25P077_RP.xlsx" PATIENTS_FILE = "patients.xlsx" COUNTS_FILE = "miRNA_counts.xlsx" OUT_EXCEL = "mirna_threshold_unified_summary.xlsx" OUT_FIG = "mirna_threshold_unified.png" THRESHOLDS = [0.0, 0.25, 1.0] YIELD_THRESHOLDS = THRESHOLDS YIELD_BUCKET_MAP = {t: t for t in THRESHOLDS} COL_HEALTHY = "#4393c3" COL_DISEASED = "#d6604d" # ── helpers ─────────────────────────────────────────────────────────────────── def parse_conc(val): if isinstance(val, str): return np.nan try: return float(val) except (TypeError, ValueError): return np.nan def extract_num(series, letter): return pd.to_numeric( series.str.extract(rf"^(\d+){letter}")[0], errors="coerce" ) # ── load data ───────────────────────────────────────────────────────────────── patients = pd.read_excel(PATIENTS_FILE, usecols=["Number", "Exf"]) patients = patients.dropna(subset=["Number"]) patients["Number"] = patients["Number"].astype(int) patients["Exf"] = patients["Exf"].astype(int) counts_mat = pd.read_excel(COUNTS_FILE, sheet_name="Mature", index_col=0) totals = counts_mat.sum(axis=0).rename("total_miRNA") qa = pd.read_excel(QC_FILE, sheet_name="A-samples") qb = pd.read_excel(QC_FILE, sheet_name="B-samples") qa["conc_val"] = qa["conc (ng/ul)"].apply(parse_conc) qb["conc_val"] = qb["Conc [ng/ul]"].apply(parse_conc) qa["Number"] = extract_num(qa["ID"], "A") qb["Number"] = extract_num(qb["ID"], "B") qa = qa.dropna(subset=["Number"]).copy(); qa["Number"] = qa["Number"].astype(int) qb = qb.dropna(subset=["Number"]).copy(); qb["Number"] = qb["Number"].astype(int) qa = qa.merge(patients, on="Number", how="left") qb = qb.merge(patients, on="Number", how="left") # ── pilot dataframe (both tissues) ─────────────────────────────────────────── def build_pilot(qc_df, tissue_letter): rows = [] for col in totals.index: m = re.match(rf"s(\d+){tissue_letter}_G2", col) if not m: continue num = int(m.group(1)) row = qc_df[qc_df["Number"] == num] if row.empty: continue rows.append({ "Number": num, "sample_id": col, "tissue": tissue_letter, "conc_val": row["conc_val"].values[0], "total_miRNA": totals[col], "Exf": row["Exf"].values[0], }) return pd.DataFrame(rows) pilot = pd.concat([build_pilot(qa, "A"), build_pilot(qb, "B")], ignore_index=True) # ── bucketed yield: pilot samples in [threshold_i, threshold_i+1) ───────────── def assign_bucket(conc): if pd.isna(conc): return -1 # "Too low" / invalid excluded from all buckets for i in range(len(YIELD_THRESHOLDS) - 1, -1, -1): if conc >= YIELD_THRESHOLDS[i]: return i return -1 pilot["bucket"] = pilot["conc_val"].apply(assign_bucket) bucket_yield = {} for i, t in enumerate(YIELD_THRESHOLDS): grp = pilot[pilot["bucket"] == i] n = len(grp) mean = grp["total_miRNA"].mean() if n else np.nan sd = grp["total_miRNA"].std() if n > 1 else np.nan mn = grp["total_miRNA"].min() if n else np.nan mx = grp["total_miRNA"].max() if n else np.nan lbl = f"≥{t}" if i == len(YIELD_THRESHOLDS) - 1 else f"{t} – <{YIELD_THRESHOLDS[i+1]}" bucket_yield[t] = { "mean": mean, "sd": sd, "min": mn, "max": mx, "n": n, "label": lbl, "first": i == 0, "last": i == len(THRESHOLDS) - 1, "samples": ", ".join(grp["sample_id"].tolist()), } def fmt_yield(yld): n, mean, sd, mn, mx = yld["n"], yld["mean"], yld["sd"], yld["min"], yld["max"] if n == 0 or np.isnan(mean): return "no pilot data" if yld["first"]: return f"<{mx/1e6:.2f}M (n={n})" if yld["last"]: return f">{mn/1e6:.2f}M (n={n})" # middle buckets: mean ± SD sd_str = f" ± {sd/1e6:.2f}M" if not np.isnan(sd) else "" return f"~{mean/1e6:.2f}M{sd_str} (n={n})" # ── cumulative cohort counts per threshold ──────────────────────────────────── def cohort_counts(qc_df): rows = [] valid = qc_df[qc_df["conc_val"].notna()] # "Too low" / invalid always excluded for t in THRESHOLDS: passing = valid if t == 0 else valid[valid["conc_val"] >= t] rows.append({ "threshold": t, "total": len(passing), "healthy": (passing["Exf"] == 0).sum(), "diseased": (passing["Exf"] == 1).sum(), }) return pd.DataFrame(rows) counts_A = cohort_counts(qa) counts_B = cohort_counts(qb) # ── pilot failures per threshold (already-run samples excluded at each cutoff) ─ def pilot_failures(pilot_tissue_df): """For each threshold, count pilot samples that FAIL (conc < threshold or NaN).""" rows = [] for t in THRESHOLDS: if t == 0: failing = pilot_tissue_df.iloc[0:0] # nothing fails at 0 else: failing = pilot_tissue_df[ pilot_tissue_df["conc_val"].isna() | (pilot_tissue_df["conc_val"] < t) ] rows.append({ "threshold": t, "fail_healthy": (failing["Exf"] == 0).sum(), "fail_diseased": (failing["Exf"] == 1).sum(), }) return pd.DataFrame(rows) pilot_A = pilot[pilot["tissue"] == "A"] pilot_B = pilot[pilot["tissue"] == "B"] pilot_fail_A = pilot_failures(pilot_A) pilot_fail_B = pilot_failures(pilot_B) def pilot_passes(pilot_tissue_df): """For each threshold, count pilot samples that PASS (conc >= threshold).""" rows = [] valid = pilot_tissue_df[pilot_tissue_df["conc_val"].notna()] for t in THRESHOLDS: passing = valid if t == 0 else valid[valid["conc_val"] >= t] rows.append({ "threshold": t, "pass_healthy": (passing["Exf"] == 0).sum(), "pass_diseased": (passing["Exf"] == 1).sum(), }) return pd.DataFrame(rows) pilot_pass_A = pilot_passes(pilot_A) pilot_pass_B = pilot_passes(pilot_B) # ── save Excel ──────────────────────────────────────────────────────────────── bucket_rows = [ {"threshold": t, "bucket": d["label"], "pilot_n": d["n"], "mean_miRNA": d["mean"], "sd_miRNA": d["sd"], "samples": d["samples"]} for t, d in bucket_yield.items() ] with pd.ExcelWriter(OUT_EXCEL, engine="openpyxl") as writer: counts_A.to_excel(writer, sheet_name="A cohort counts", index=False) counts_B.to_excel(writer, sheet_name="B cohort counts", index=False) pd.DataFrame(bucket_rows).to_excel(writer, sheet_name="Bucketed yield", index=False) for ws in writer.sheets.values(): for col_cells in ws.columns: w = max(len(str(c.value)) if c.value else 0 for c in col_cells) ws.column_dimensions[col_cells[0].column_letter].width = w + 4 print(f"Excel saved → {OUT_EXCEL}") # ── figure ──────────────────────────────────────────────────────────────────── # 4 bars per threshold: A-healthy | A-diseased | B-healthy | B-diseased # Each bar is stacked into 3 segments (bottom → top): # 1. New samples passing — solid colour # 2. Pilot samples passing — same colour, hatched "////" # 3. Pilot samples failing — light grey, hatched "xxxx" (already run, excluded) fig, ax = plt.subplots(figsize=(18, 6)) fig.suptitle("miRNA QC Threshold Analysis", fontsize=13, fontweight="bold", y=1.01) n_thresh = len(THRESHOLDS) group_w = 0.72 bar_w = group_w / 4 offsets = np.array([-1.5, -0.5, 0.5, 1.5]) * bar_w x = np.arange(n_thresh) COL_PILOT_FAIL = "#bbbbbb" # muted grey for already-run failures # For each of the 4 bar positions we need three value arrays (length = n_thresh): # new_pass, pilot_pass, pilot_fail def split_counts(counts_df, ppass_df, pfail_df, diagnosis): """Return (new_pass, pilot_pass, pilot_fail) arrays across thresholds.""" total = counts_df[diagnosis].values pp = ppass_df[f"pass_{diagnosis}"].values pf = pfail_df[f"fail_{diagnosis}"].values new_pass = total - pp # non-pilot samples passing return new_pass, pp, pf specs = [ # (counts_df, ppass_df, pfail_df, diag, color, tissue_hatch, label_base) (counts_A, pilot_pass_A, pilot_fail_A, "healthy", COL_HEALTHY, False, "A — Healthy"), (counts_A, pilot_pass_A, pilot_fail_A, "diseased", COL_DISEASED, False, "A — Diseased"), (counts_B, pilot_pass_B, pilot_fail_B, "healthy", COL_HEALTHY, True, "B — Healthy"), (counts_B, pilot_pass_B, pilot_fail_B, "diseased", COL_DISEASED, True, "B — Diseased"), ] # Track whether we've added each legend entry already legend_added = set() for (c_df, pp_df, pf_df, diag, color, is_B, lbl), offset in zip(specs, offsets): new_pass, pilot_pass, pilot_fail = split_counts(c_df, pp_df, pf_df, diag) xpos = x + offset tissue_hatch = "////" if is_B else "" # ── segment 1: new samples passing ─────────────────────────────────────── key1 = f"new_{diag}_{'B' if is_B else 'A'}" ax.bar(xpos, new_pass, width=bar_w, color=color, hatch=tissue_hatch, alpha=0.85, edgecolor="white" if is_B else color, linewidth=0.5, zorder=3, label=lbl if key1 not in legend_added else "_nolegend_") legend_added.add(key1) # ── segment 2: pilot samples passing (stacked, same colour, denser hatch) ─ pilot_hatch = "////////" if is_B else "////" ax.bar(xpos, pilot_pass, width=bar_w, bottom=new_pass, color=color, hatch=pilot_hatch, alpha=0.55, edgecolor="white", linewidth=0.5, zorder=3, label="Pilot — passing (already run)" if "pilot_pass" not in legend_added else "_nolegend_") legend_added.add("pilot_pass") # ── segment 3: pilot samples failing (stacked above, grey) ─────────────── pilot_fail_bottom = new_pass + pilot_pass ax.bar(xpos, pilot_fail, width=bar_w, bottom=pilot_fail_bottom, color=COL_PILOT_FAIL, hatch="xxxx", alpha=0.7, edgecolor="white", linewidth=0.5, zorder=3, label="Pilot — failing (already run)" if "pilot_fail" not in legend_added else "_nolegend_") legend_added.add("pilot_fail") # yield annotations y_max = max( (counts_A[["healthy","diseased"]] + pilot_fail_A[["fail_healthy","fail_diseased"]].values).max().max(), (counts_B[["healthy","diseased"]] + pilot_fail_B[["fail_healthy","fail_diseased"]].values).max().max(), ) annot_y = y_max * 1.08 for i, t in enumerate(THRESHOLDS): yld = bucket_yield[YIELD_BUCKET_MAP[t]] txt = fmt_yield(yld) ax.text(i, annot_y, txt, ha="center", va="bottom", fontsize=6.5, color="#444444", bbox=dict(boxstyle="round,pad=0.25", fc="lightyellow", ec="goldenrod", alpha=0.8, lw=0.6)) ax.set_xticks(x) ax.set_xticklabels([str(t) for t in THRESHOLDS], fontsize=9) ax.set_xlabel("Min. concentration threshold (ng/µl)", fontsize=10) ax.set_ylabel("N samples", fontsize=10) ax.set_ylim(0, annot_y * 1.22) ax.tick_params(labelsize=9) legend_handles = [ Patch(facecolor="gray", edgecolor="gray", hatch="", label="A — Schirmer strips"), Patch(facecolor="gray", edgecolor="white", hatch="////", label="B — Lens tissues"), Patch(facecolor=COL_HEALTHY, label="Healthy"), Patch(facecolor=COL_DISEASED, label="Diseased"), Patch(facecolor="gray", edgecolor="white", hatch="////", alpha=0.55, label="Pilot — passing"), Patch(facecolor=COL_PILOT_FAIL,edgecolor="white", hatch="xxxx", alpha=0.7, label="Pilot — failing"), ] ax.legend(handles=legend_handles, fontsize=7.5, framealpha=0.85, loc="upper right", ncol=2) ax.text(0.01, 0.99, "Bars: new samples passing (solid) + pilot passing (light hatch) + pilot failing (grey, above)\n" "Annotations: approx. expected miRNA yield (pilot bucket mean, A + B combined)", transform=ax.transAxes, fontsize=6, va="top", color="gray", linespacing=1.5) # ── per-threshold report table below x-axis ─────────────────────────────────── # Uses blended transform: data coords for x, axes fraction for y (negative = below axis) from matplotlib.transforms import blended_transform_factory trans = blended_transform_factory(ax.transData, ax.transAxes) table_rows = [] for i in range(len(THRESHOLDS)): ah_tot = int(counts_A["healthy"].iloc[i]) bh_tot = int(counts_B["healthy"].iloc[i]) ag_tot = int(counts_A["diseased"].iloc[i]) bg_tot = int(counts_B["diseased"].iloc[i]) ah_rem = ah_tot - int(pilot_pass_A["pass_healthy"].iloc[i]) bh_rem = bh_tot - int(pilot_pass_B["pass_healthy"].iloc[i]) ag_rem = ag_tot - int(pilot_pass_A["pass_diseased"].iloc[i]) bg_rem = bg_tot - int(pilot_pass_B["pass_diseased"].iloc[i]) table_rows.append((ah_tot, bh_tot, ag_tot, bg_tot, ah_rem, bh_rem, ag_rem, bg_rem)) for i, t in enumerate(THRESHOLDS): ah_tot, bh_tot, ag_tot, bg_tot, ah_rem, bh_rem, ag_rem, bg_rem = table_rows[i] def cell(val, prev=None): if prev is None: return f"{val:2d} " d = val - prev return f"{val:2d} (Δ{d:+d})" if d != 0 else f"{val:2d} (Δ0) " if i == 0: line1 = (f"Total: AH: {cell(ah_tot)} | BH: {cell(bh_tot)} | " f"AG: {cell(ag_tot)} | BG: {cell(bg_tot)}") line2 = (f"Remaining: AH: {cell(ah_rem)} | BH: {cell(bh_rem)} | " f"AG: {cell(ag_rem)} | BG: {cell(bg_rem)}") else: p = table_rows[i - 1] line1 = (f"Total: AH: {cell(ah_tot, p[0])} | BH: {cell(bh_tot, p[1])} | " f"AG: {cell(ag_tot, p[2])} | BG: {cell(bg_tot, p[3])}") line2 = (f"Remaining: AH: {cell(ah_rem, p[4])} | BH: {cell(bh_rem, p[5])} | " f"AG: {cell(ag_rem, p[6])} | BG: {cell(bg_rem, p[7])}") total_rem = ah_rem + bh_rem + ag_rem + bg_rem line3 = f"Total remaining to sequence: {total_rem}" table_txt = line1 + "\n" + line2 + "\n" + line3 ax.text(i, -0.18, table_txt, transform=trans, ha="center", va="top", fontsize=7, family="monospace", color="#222222", bbox=dict(boxstyle="round,pad=0.4", fc="white", ec="#cccccc", lw=0.8), clip_on=False) fig.subplots_adjust(bottom=0.28) fig.savefig(OUT_FIG, dpi=150, bbox_inches="tight") print(f"Figure saved → {OUT_FIG}") plt.show()