rpotter6298
257 lines
11 KiB
Python
257 lines
11 KiB
Python
"""
|
||
miRNA threshold analysis
|
||
------------------------
|
||
For each tissue (A = Schirmer strips, B = Lens tissues):
|
||
- Derives concentration thresholds from the pilot group's actual measured values
|
||
- At each threshold shows: # samples passing from full cohort (healthy / diseased)
|
||
and the mean ± SD total miRNA counts from pilot samples that pass that cut-off
|
||
Outputs: summary Excel table + multi-panel figure
|
||
"""
|
||
|
||
import re
|
||
import numpy as np
|
||
import pandas as pd
|
||
import matplotlib.pyplot as plt
|
||
import matplotlib.gridspec as gridspec
|
||
from matplotlib.lines import Line2D
|
||
|
||
# ── file paths ────────────────────────────────────────────────────────────────
|
||
QC_FILE = "BEA25P077_RP.xlsx"
|
||
PATIENTS_FILE = "patients.xlsx"
|
||
COUNTS_FILE = "miRNA_counts.xlsx"
|
||
OUT_EXCEL = "mirna_threshold_summary.xlsx"
|
||
OUT_FIG = "mirna_threshold_analysis.png"
|
||
|
||
# ── colours ───────────────────────────────────────────────────────────────────
|
||
COL_HEALTHY = "#4393c3" # blue
|
||
COL_DISEASED = "#d6604d" # red
|
||
COL_PILOT_A = "#2ca02c" # green
|
||
COL_PILOT_B = "#9467bd" # purple
|
||
|
||
# ── 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):
|
||
"""Pull the patient number from IDs like '4A_G2'."""
|
||
extracted = series.str.extract(rf"^(\d+){letter}")[0]
|
||
return pd.to_numeric(extracted, 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 = pd.read_excel(COUNTS_FILE, sheet_name="Mature", index_col=0)
|
||
totals = counts.sum(axis=0).rename("total_miRNA") # total per sample
|
||
|
||
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").dropna().astype(int)
|
||
qb["Number"] = extract_num(qb["ID"], "B").dropna().astype(int)
|
||
|
||
qa = qa.dropna(subset=["Number"]).merge(patients, on="Number", how="left")
|
||
qb = qb.dropna(subset=["Number"]).merge(patients, on="Number", how="left")
|
||
|
||
|
||
# ── build pilot lookup: number → (conc, total_miRNA) ─────────────────────────
|
||
def build_pilot(qc_df, tissue_letter):
|
||
"""Join QC concentrations with miRNA totals for pilot samples."""
|
||
pilot_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))
|
||
conc_row = qc_df[qc_df["Number"] == num]
|
||
if conc_row.empty:
|
||
continue
|
||
pilot_rows.append({
|
||
"Number": num,
|
||
"sample_id": col,
|
||
"conc_val": conc_row["conc_val"].values[0],
|
||
"total_miRNA": totals[col],
|
||
"Exf": conc_row["Exf"].values[0],
|
||
})
|
||
return pd.DataFrame(pilot_rows).sort_values("conc_val")
|
||
|
||
|
||
pilot_A = build_pilot(qa, "A")
|
||
pilot_B = build_pilot(qb, "B")
|
||
|
||
|
||
# ── derive thresholds from pilot concentrations ───────────────────────────────
|
||
def pilot_thresholds(pilot_df):
|
||
"""
|
||
Use the sorted unique measured concentrations as tier cut-offs,
|
||
prepend 0 (all samples) and append a step just above the max.
|
||
"""
|
||
measured = sorted(pilot_df["conc_val"].dropna().unique())
|
||
# round to 2 dp to keep labels clean
|
||
thresholds = [0.0] + [round(v, 2) for v in measured]
|
||
return thresholds
|
||
|
||
|
||
thresholds_A = pilot_thresholds(pilot_A)
|
||
thresholds_B = pilot_thresholds(pilot_B)
|
||
|
||
|
||
# ── per-threshold summary ─────────────────────────────────────────────────────
|
||
def threshold_summary(qc_df, pilot_df, thresholds):
|
||
rows = []
|
||
for t in thresholds:
|
||
# full cohort
|
||
if t == 0:
|
||
cohort = qc_df.copy()
|
||
else:
|
||
cohort = qc_df[qc_df["conc_val"] >= t]
|
||
n_total = len(cohort)
|
||
n_healthy = (cohort["Exf"] == 0).sum()
|
||
n_diseased = (cohort["Exf"] == 1).sum()
|
||
|
||
# pilot subset
|
||
if t == 0:
|
||
pilot_pass = pilot_df.copy()
|
||
else:
|
||
pilot_pass = pilot_df[pilot_df["conc_val"] >= t]
|
||
n_pilot = len(pilot_pass)
|
||
mean_mirna = pilot_pass["total_miRNA"].mean() if n_pilot else np.nan
|
||
sd_mirna = pilot_pass["total_miRNA"].std() if n_pilot > 1 else np.nan
|
||
|
||
rows.append({
|
||
"threshold (ng/µl)": t,
|
||
"cohort_total": n_total,
|
||
"cohort_healthy": n_healthy,
|
||
"cohort_diseased": n_diseased,
|
||
"pilot_n": n_pilot,
|
||
"pilot_mean_miRNA": mean_mirna,
|
||
"pilot_sd_miRNA": sd_mirna,
|
||
})
|
||
return pd.DataFrame(rows)
|
||
|
||
|
||
summary_A = threshold_summary(qa, pilot_A, thresholds_A)
|
||
summary_B = threshold_summary(qb, pilot_B, thresholds_B)
|
||
|
||
|
||
# ── save Excel ────────────────────────────────────────────────────────────────
|
||
with pd.ExcelWriter(OUT_EXCEL, engine="openpyxl") as writer:
|
||
for df, sheet in [(summary_A, "A - Schirmer strips"),
|
||
(summary_B, "B - Lens tissues")]:
|
||
df.to_excel(writer, sheet_name=sheet, index=False)
|
||
ws = writer.sheets[sheet]
|
||
# widen columns
|
||
for col_cells in ws.columns:
|
||
max_len = max(len(str(c.value)) if c.value else 0 for c in col_cells)
|
||
ws.column_dimensions[col_cells[0].column_letter].width = max_len + 4
|
||
print(f"Excel saved → {OUT_EXCEL}")
|
||
|
||
|
||
# ── figure ────────────────────────────────────────────────────────────────────
|
||
fig = plt.figure(figsize=(16, 10))
|
||
fig.suptitle("miRNA QC Threshold Analysis", fontsize=14, fontweight="bold", y=0.98)
|
||
|
||
# 2 tissues × 3 columns: scatter | cohort counts | expected miRNA
|
||
gs = gridspec.GridSpec(2, 3, figure=fig, hspace=0.45, wspace=0.38,
|
||
left=0.07, right=0.97, top=0.92, bottom=0.08)
|
||
|
||
|
||
def plot_tissue(row, tissue_label, pilot_df, summary_df, pilot_col):
|
||
thresholds = summary_df["threshold (ng/µl)"].tolist()
|
||
x_labels = [str(t) for t in thresholds]
|
||
x = np.arange(len(thresholds))
|
||
bar_w = 0.38
|
||
|
||
# ── panel 1: scatter concentration vs total miRNA ─────────────────────────
|
||
ax1 = fig.add_subplot(gs[row, 0])
|
||
healthy_mask = pilot_df["Exf"] == 0
|
||
diseased_mask = pilot_df["Exf"] == 1
|
||
|
||
ax1.scatter(pilot_df.loc[healthy_mask, "conc_val"],
|
||
pilot_df.loc[healthy_mask, "total_miRNA"] / 1e6,
|
||
color=COL_HEALTHY, edgecolors="k", linewidths=0.5,
|
||
s=70, zorder=3, label="Healthy")
|
||
ax1.scatter(pilot_df.loc[diseased_mask, "conc_val"],
|
||
pilot_df.loc[diseased_mask, "total_miRNA"] / 1e6,
|
||
color=COL_DISEASED, edgecolors="k", linewidths=0.5,
|
||
s=70, zorder=3, label="Diseased")
|
||
|
||
# annotate patient numbers
|
||
for _, r in pilot_df.iterrows():
|
||
if pd.notna(r["conc_val"]):
|
||
ax1.annotate(str(int(r["Number"])),
|
||
(r["conc_val"], r["total_miRNA"] / 1e6),
|
||
textcoords="offset points", xytext=(4, 3),
|
||
fontsize=7, color="dimgray")
|
||
|
||
# vertical lines at each threshold (skip 0)
|
||
for t in thresholds[1:]:
|
||
ax1.axvline(t, color="gray", lw=0.7, ls="--", alpha=0.5)
|
||
|
||
ax1.set_xlabel("Concentration (ng/µl)", fontsize=9)
|
||
ax1.set_ylabel("Total miRNA counts (×10⁶)", fontsize=9)
|
||
ax1.set_title(f"{tissue_label}\nPilot: conc vs miRNA yield", fontsize=9)
|
||
ax1.legend(fontsize=8, framealpha=0.7)
|
||
ax1.tick_params(labelsize=8)
|
||
|
||
# ── panel 2: cohort counts per threshold ──────────────────────────────────
|
||
ax2 = fig.add_subplot(gs[row, 1])
|
||
ax2.bar(x - bar_w / 2, summary_df["cohort_healthy"],
|
||
width=bar_w, color=COL_HEALTHY, label="Healthy", alpha=0.85)
|
||
ax2.bar(x + bar_w / 2, summary_df["cohort_diseased"],
|
||
width=bar_w, color=COL_DISEASED, label="Diseased", alpha=0.85)
|
||
ax2.set_xticks(x)
|
||
ax2.set_xticklabels(x_labels, rotation=45, ha="right", fontsize=8)
|
||
ax2.set_xlabel("Min. concentration threshold (ng/µl)", fontsize=9)
|
||
ax2.set_ylabel("N samples passing", fontsize=9)
|
||
ax2.set_title(f"{tissue_label}\nCohort samples at each threshold", fontsize=9)
|
||
ax2.legend(fontsize=8, framealpha=0.7)
|
||
ax2.tick_params(labelsize=8)
|
||
# add count labels
|
||
for bar in ax2.patches:
|
||
h = bar.get_height()
|
||
if h > 0:
|
||
ax2.text(bar.get_x() + bar.get_width() / 2, h + 0.3,
|
||
str(int(h)), ha="center", va="bottom", fontsize=7)
|
||
|
||
# ── panel 3: expected miRNA from pilot ────────────────────────────────────
|
||
ax3 = fig.add_subplot(gs[row, 2])
|
||
means = summary_df["pilot_mean_miRNA"] / 1e6
|
||
sds = summary_df["pilot_sd_miRNA"].fillna(0) / 1e6
|
||
ns = summary_df["pilot_n"]
|
||
|
||
ax3.bar(x, means, width=0.55, color=pilot_col, alpha=0.8, label="Mean ± SD")
|
||
ax3.errorbar(x, means, yerr=sds, fmt="none", color="k",
|
||
capsize=4, linewidth=1.2, zorder=4)
|
||
|
||
# label n= above each bar
|
||
for i, (m, n) in enumerate(zip(means, ns)):
|
||
ax3.text(i, (m + sds.iloc[i]) + 0.05, f"n={int(n)}",
|
||
ha="center", va="bottom", fontsize=7, color="dimgray")
|
||
|
||
ax3.set_xticks(x)
|
||
ax3.set_xticklabels(x_labels, rotation=45, ha="right", fontsize=8)
|
||
ax3.set_xlabel("Min. concentration threshold (ng/µl)", fontsize=9)
|
||
ax3.set_ylabel("Total miRNA counts (×10⁶)", fontsize=9)
|
||
ax3.set_title(f"{tissue_label}\nExpected miRNA yield (pilot mean ± SD)", fontsize=9)
|
||
ax3.tick_params(labelsize=8)
|
||
|
||
|
||
plot_tissue(0, "A — Schirmer strips", pilot_A, summary_A, COL_PILOT_A)
|
||
plot_tissue(1, "B — Lens tissues", pilot_B, summary_B, COL_PILOT_B)
|
||
|
||
fig.savefig(OUT_FIG, dpi=150, bbox_inches="tight")
|
||
print(f"Figure saved → {OUT_FIG}")
|
||
plt.show()
|