rpotter6298
371 lines
16 KiB
Python
371 lines
16 KiB
Python
"""
|
||
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()
|