Add new Excel report BEA25P077_RP.xlsx with multiple worksheets and styles

Pilot/code/generate_sequencing_plan.py

@@ -0,0 +1,88 @@
+"""
+Generate sample_sequencing_plan.xlsx
+Labels every sample as Already sequenced / To be sequenced / Skip
+at the chosen concentration threshold, with colour-coded rows.
+"""
+
+import numpy as np
+import pandas as pd
+from openpyxl import load_workbook
+from openpyxl.styles import PatternFill
+
+THRESHOLD  = 0.25
+OUT_FILE   = "sample_sequencing_plan.xlsx"
+
+PILOT_IDS  = {("A", n) for n in [1, 2, 3, 4, 5, 6]} | \
+             {("B", n) for n in [1, 2, 15, 16, 17, 18]}
+
+FILLS = {
+    "Already sequenced": PatternFill("solid", fgColor="C6EFCE"),  # green
+    "To be sequenced":   PatternFill("solid", fgColor="FFEB9C"),  # yellow
+    "Skip":              PatternFill("solid", fgColor="FFC7CE"),  # red/pink
+}
+
+# ── load data ─────────────────────────────────────────────────────────────────
+patients = pd.read_excel("patients.xlsx", usecols=["Number", "Exf"])
+patients = patients.dropna(subset=["Number"])
+patients["Number"] = patients["Number"].astype(int)
+patients["Exf"]    = patients["Exf"].astype(int)
+
+def parse_conc(v):
+    if isinstance(v, str): return np.nan
+    try: return float(v)
+    except: return np.nan
+
+qa = pd.read_excel("BEA25P077_RP.xlsx", sheet_name="A-samples")
+qb = pd.read_excel("BEA25P077_RP.xlsx", 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"] = pd.to_numeric(qa["ID"].str.extract(r"^(\d+)A")[0], errors="coerce")
+qb["Number"] = pd.to_numeric(qb["ID"].str.extract(r"^(\d+)B")[0], errors="coerce")
+qa = qa.dropna(subset=["Number"]); qa["Number"] = qa["Number"].astype(int); qa["tissue"] = "A"
+qb = qb.dropna(subset=["Number"]); qb["Number"] = qb["Number"].astype(int); qb["tissue"] = "B"
+qa = qa.merge(patients, on="Number", how="left")
+qb = qb.merge(patients, on="Number", how="left")
+
+# ── build plan ────────────────────────────────────────────────────────────────
+rows = []
+for _, r in pd.concat([qa, qb], ignore_index=True).iterrows():
+    tissue, num, conc = r["tissue"], r["Number"], r["conc_val"]
+    if (tissue, num) in PILOT_IDS:
+        status = "Already sequenced"
+    elif pd.isna(conc) or conc < THRESHOLD:
+        status = "Skip"
+    else:
+        status = "To be sequenced"
+    rows.append({
+        "Sample":       f"{num}{tissue}_G2",
+        "Patient":      num,
+        "Tissue":       tissue,
+        "Conc (ng/ul)": conc,
+        "Status":       status,
+    })
+
+df = pd.DataFrame(rows).sort_values(["Tissue", "Patient"])
+
+# ── write Excel ───────────────────────────────────────────────────────────────
+df.to_excel(OUT_FILE, index=False)
+
+wb = load_workbook(OUT_FILE)
+ws = wb.active
+
+# auto-width columns
+for col_cells in ws.columns:
+    width = max(len(str(c.value)) if c.value is not None else 0 for c in col_cells)
+    ws.column_dimensions[col_cells[0].column_letter].width = width + 4
+
+# colour every row based on Status (Status is the last column)
+status_col = df.columns.get_loc("Status") + 1   # 1-indexed
+for row in ws.iter_rows(min_row=2):
+    status = row[status_col - 1].value
+    fill = FILLS.get(status)
+    if fill:
+        for cell in row:
+            cell.fill = fill
+
+wb.save(OUT_FILE)
+print(f"Saved {OUT_FILE}  ({len(df)} samples)")
+print(df["Status"].value_counts().to_string())

Pilot/code/miRNA_high_expression.py

@@ -0,0 +1,589 @@
+#!/usr/bin/env python3
+import argparse
+import csv
+import os
+import re
+import ssl
+import statistics
+import sys
+import urllib.error
+import urllib.parse
+import urllib.request
+import zipfile
+import xml.etree.ElementTree as ET
+from collections import defaultdict
+
+
+def col_to_index(col):
+    idx = 0
+    for c in col:
+        idx = idx * 26 + (ord(c.upper()) - ord("A") + 1)
+    return idx
+
+
+def parse_shared_strings(zf):
+    try:
+        xml = zf.read("xl/sharedStrings.xml")
+    except KeyError:
+        return []
+    root = ET.fromstring(xml)
+    ns = {"a": root.tag.split("}")[0].strip("{")}
+    shared = []
+    for si in root.findall(".//a:si", ns):
+        text_parts = []
+        for t in si.findall(".//a:t", ns):
+            text_parts.append(t.text or "")
+        shared.append("".join(text_parts))
+    return shared
+
+
+def list_sheets(zf):
+    wb = ET.fromstring(zf.read("xl/workbook.xml"))
+    ns = {"a": wb.tag.split("}")[0].strip("{")}
+    sheets = []
+    for sh in wb.findall(".//a:sheets/a:sheet", ns):
+        sheets.append(
+            (sh.attrib["name"], sh.attrib["{http://schemas.openxmlformats.org/officeDocument/2006/relationships}id"])
+        )
+
+    rels = ET.fromstring(zf.read("xl/_rels/workbook.xml.rels"))
+    relmap = {}
+    for rel in rels.findall(".//{http://schemas.openxmlformats.org/package/2006/relationships}Relationship"):
+        relmap[rel.attrib["Id"]] = rel.attrib["Target"]
+
+    return [(name, "xl/" + relmap[rid]) for name, rid in sheets]
+
+
+def read_xlsx_sheet(path, sheet_name):
+    with zipfile.ZipFile(path) as zf:
+        shared = parse_shared_strings(zf)
+        sheets = dict(list_sheets(zf))
+        if sheet_name not in sheets:
+            raise ValueError(f"Sheet '{sheet_name}' not found. Available: {', '.join(sheets.keys())}")
+
+        root = ET.fromstring(zf.read(sheets[sheet_name]))
+        ns = {"a": root.tag.split("}")[0].strip("{")}
+        row_dicts = []
+        max_col_idx = 0
+
+        for row in root.findall(".//a:sheetData/a:row", ns):
+            row_cells = {}
+            for c in row.findall("a:c", ns):
+                ref = c.attrib.get("r", "")
+                col = "".join(ch for ch in ref if ch.isalpha())
+                if not col:
+                    continue
+                v = c.find("a:v", ns)
+                value = v.text if v is not None else ""
+                if c.attrib.get("t") == "s":
+                    try:
+                        value = shared[int(value)]
+                    except Exception:
+                        pass
+                row_cells[col] = value
+                max_col_idx = max(max_col_idx, col_to_index(col))
+            if row_cells:
+                row_dicts.append(row_cells)
+
+        if not row_dicts:
+            return []
+
+        rows = []
+        for row_cells in row_dicts:
+            row_list = [""] * max_col_idx
+            for col, val in row_cells.items():
+                row_list[col_to_index(col) - 1] = val
+            rows.append(row_list)
+
+        header = rows[0]
+        ncols = len(header)
+        while ncols > 0 and header[ncols - 1] == "":
+            ncols -= 1
+
+        normalized = []
+        for row in rows:
+            if len(row) < ncols:
+                row = row + [""] * (ncols - len(row))
+            elif len(row) > ncols:
+                row = row[:ncols]
+            normalized.append(row)
+        return normalized
+
+
+def score_counts(values, metric):
+    if not values:
+        return 0.0
+    if metric == "mean":
+        return sum(values) / len(values)
+    if metric == "median":
+        return statistics.median(values)
+    if metric == "sum":
+        return sum(values)
+    if metric == "max":
+        return max(values)
+    raise ValueError(f"Unknown metric: {metric}")
+
+
+def canonical_mirna(name):
+    s = name.strip().lower()
+    s = re.sub(r"^[a-z]{3}-", "", s)
+    return s
+
+
+def _detect_delimiter(sample):
+    if "\t" in sample:
+        return "\t"
+    if "," in sample:
+        return ","
+    return "\t"
+
+
+def read_delimited_rows(path):
+    with open(path, "r", newline="") as f:
+        sample = f.readline()
+        delim = _detect_delimiter(sample)
+        f.seek(0)
+        try:
+            reader = csv.reader(f, delimiter=delim)
+            rows = list(reader)
+        except csv.Error:
+            f.seek(0)
+            reader = csv.reader(
+                f,
+                delimiter=delim,
+                quoting=csv.QUOTE_NONE,
+                escapechar="\\",
+            )
+            rows = list(reader)
+    return rows
+
+
+def iter_delimited_rows(path):
+    with open(path, "r", newline="") as f:
+        sample = f.readline()
+        delim = _detect_delimiter(sample)
+        f.seek(0)
+        try:
+            reader = csv.reader(f, delimiter=delim)
+            for row in reader:
+                yield row
+        except csv.Error:
+            f.seek(0)
+            reader = csv.reader(
+                f,
+                delimiter=delim,
+                quoting=csv.QUOTE_NONE,
+                escapechar="\\",
+            )
+            for row in reader:
+                yield row
+
+
+def read_targets(path):
+    rows = read_delimited_rows(path)
+    return parse_targets_rows(rows)
+
+
+def normalize_header(col):
+    return re.sub(r"[^a-z0-9]+", "", col.strip().lower())
+
+
+def normalize_species(value):
+    v = re.sub(r"[^a-z]+", " ", value.strip().lower())
+    v = re.sub(r"\s+", " ", v).strip()
+    if v in {"human", "homo sapiens", "homo sapienst"}:
+        return "homo sapiens"
+    return v
+
+
+def download_url_to_file(url, dest_path, allow_insecure=False):
+    os.makedirs(os.path.dirname(dest_path) or ".", exist_ok=True)
+    req = urllib.request.Request(url, headers={"User-Agent": "Mozilla/5.0"})
+    context = ssl._create_unverified_context() if allow_insecure else None
+    if context:
+        response = urllib.request.urlopen(req, context=context)
+    else:
+        response = urllib.request.urlopen(req)
+    with response as resp, open(dest_path, "wb") as f:
+        while True:
+            chunk = resp.read(1024 * 1024)
+            if not chunk:
+                break
+            f.write(chunk)
+
+
+def find_mirtarbase_mti_url(download_page_url):
+    req = urllib.request.Request(download_page_url, headers={"User-Agent": "Mozilla/5.0"})
+    with urllib.request.urlopen(req) as resp:
+        html = resp.read().decode("utf-8", errors="ignore")
+    links = re.findall(r'href=["\\\']([^"\\\']+)["\\\']', html, flags=re.IGNORECASE)
+    scored = []
+    for link in links:
+        if "mirtarbase" not in link.lower():
+            continue
+        if "mti" not in link.lower():
+            continue
+        if not re.search(r"\\.(txt|tsv|csv|xls|xlsx|zip)$", link, flags=re.IGNORECASE):
+            continue
+        scored.append(link)
+    if not scored:
+        return None
+    scored.sort(key=lambda s: (0 if s.lower().endswith(".txt") else 1, len(s)))
+    return urllib.parse.urljoin(download_page_url, scored[0])
+
+
+def read_mirtarbase(
+    path,
+    species="Homo sapiens",
+    support=None,
+    sheet=None,
+    include_mirnas=None,
+):
+    species_norm = normalize_species(species) if species else None
+    ext = os.path.splitext(path)[1].lower()
+    if ext == ".xlsx":
+        if sheet is None:
+            with zipfile.ZipFile(path) as zf:
+                sheets = list_sheets(zf)
+                if not sheets:
+                    raise ValueError("miRTarBase .xlsx has no sheets.")
+                sheet = sheets[0][0]
+        rows = read_xlsx_sheet(path, sheet)
+        if not rows:
+            return {"raw": defaultdict(set), "canonical": defaultdict(set)}
+        header = rows[0]
+        row_iter = iter(rows[1:])
+    else:
+        row_iter = iter_delimited_rows(path)
+        try:
+            header = next(row_iter)
+        except StopIteration:
+            return {"raw": defaultdict(set), "canonical": defaultdict(set)}
+
+    idx_mirna = None
+    idx_target = None
+    idx_species_target = None
+    idx_species_mirna = None
+    idx_support = None
+
+    for i, col in enumerate(header):
+        key = normalize_header(col)
+        if key in {"mirna", "mirnaname", "mirnaid", "mirname"} and idx_mirna is None:
+            idx_mirna = i
+        if key in {"targetgene", "target", "gene", "genesymbol"} and idx_target is None:
+            idx_target = i
+        if key in {"speciestargetgene", "speciestarget", "targetspecies"} and idx_species_target is None:
+            idx_species_target = i
+        if key in {"speciesmirna", "speciessourcemiRNA", "speciesmir"} and idx_species_mirna is None:
+            idx_species_mirna = i
+        if key in {"supporttype", "support", "evidence"} and idx_support is None:
+            idx_support = i
+
+    if idx_mirna is None or idx_target is None:
+        raise ValueError(
+            "miRTarBase format not recognized. Expected columns like 'miRNA' and 'Target Gene'."
+        )
+
+    mapping = defaultdict(set)
+    support_norm = support.strip().lower() if support else None
+
+    include_canonical = None
+    if include_mirnas:
+        include_canonical = {canonical_mirna(m) for m in include_mirnas}
+
+    for row in row_iter:
+        if len(row) <= max(idx_mirna, idx_target):
+            continue
+        mirna = row[idx_mirna].strip()
+        target = row[idx_target].strip()
+        if not mirna or not target:
+            continue
+
+        if include_canonical is not None and canonical_mirna(mirna) not in include_canonical:
+            continue
+
+        if species_norm:
+            species_value = ""
+            if idx_species_target is not None and idx_species_target < len(row):
+                species_value = row[idx_species_target]
+            elif idx_species_mirna is not None and idx_species_mirna < len(row):
+                species_value = row[idx_species_mirna]
+            if species_value:
+                if normalize_species(species_value) != species_norm:
+                    continue
+
+        if support_norm and idx_support is not None and idx_support < len(row):
+            if support_norm not in row[idx_support].lower():
+                continue
+
+        mapping[mirna].add(target)
+
+    canonical = defaultdict(set)
+    for mirna, targets in mapping.items():
+        canonical[canonical_mirna(mirna)].update(targets)
+
+    return {"raw": mapping, "canonical": canonical}
+def parse_targets_rows(rows):
+    if not rows:
+        return {}
+
+    header = rows[0]
+    idx_mirna = None
+    idx_target = None
+    for i, col in enumerate(header):
+        key = normalize_header(col)
+        if key in {"mirna", "mir", "mirname", "mirid"} and idx_mirna is None:
+            idx_mirna = i
+        if key in {"target", "gene", "protein", "symbol", "targetgene", "targets"} and idx_target is None:
+            idx_target = i
+
+    if idx_mirna is None or idx_target is None:
+        idx_mirna, idx_target = 0, 1
+
+    mapping = defaultdict(set)
+    for row in rows[1:]:
+        if len(row) <= max(idx_mirna, idx_target):
+            continue
+        mirna = row[idx_mirna].strip()
+        target_cell = row[idx_target].strip()
+        if not mirna or not target_cell:
+            continue
+        targets = [target_cell]
+        if ";" in target_cell:
+            targets = [t.strip() for t in target_cell.split(";") if t.strip()]
+        elif "," in target_cell:
+            targets = [t.strip() for t in target_cell.split(",") if t.strip()]
+
+        for target in targets:
+            mapping[mirna].add(target)
+
+    canonical = defaultdict(set)
+    for mirna, targets in mapping.items():
+        canonical[canonical_mirna(mirna)].update(targets)
+
+    return {"raw": mapping, "canonical": canonical}
+
+
+def select_high_expression(rows, metric, top_n, min_score, quantile):
+    header = rows[0]
+    sample_cols = header[1:]
+    entries = []
+
+    for row in rows[1:]:
+        if not row or not row[0]:
+            continue
+        counts = []
+        for val in row[1:]:
+            try:
+                counts.append(float(val))
+            except Exception:
+                counts.append(0.0)
+        score = score_counts(counts, metric)
+        entries.append((row[0], counts, score))
+
+    entries.sort(key=lambda x: x[2], reverse=True)
+    scores = [e[2] for e in entries]
+    selected = entries
+
+    if top_n is not None:
+        selected = entries[: max(0, top_n)]
+    elif quantile is not None:
+        if not scores:
+            selected = []
+        else:
+            idx = int(round((len(scores) - 1) * quantile))
+            thresh = sorted(scores)[idx]
+            selected = [e for e in entries if e[2] >= thresh]
+    elif min_score is not None:
+        selected = [e for e in entries if e[2] >= min_score]
+    else:
+        selected = entries[:50]
+
+    return sample_cols, selected
+
+
+def write_high_expression(path, sample_cols, selected, metric):
+    os.makedirs(os.path.dirname(path), exist_ok=True)
+    with open(path, "w", newline="") as f:
+        writer = csv.writer(f, delimiter="\t")
+        writer.writerow(["miRNA", f"{metric}_score", *sample_cols])
+        for mirna, counts, score in selected:
+            writer.writerow([mirna, f"{score:.6g}", *counts])
+
+
+def write_targets(outdir, selected, target_map, metric):
+    edges_path = os.path.join(outdir, "miRNA_protein_edges.tsv")
+    summary_path = os.path.join(outdir, "protein_targets_summary.tsv")
+    os.makedirs(outdir, exist_ok=True)
+
+    target_stats = defaultdict(lambda: {"miRNAs": set(), "score_sum": 0.0})
+    with open(edges_path, "w", newline="") as f:
+        writer = csv.writer(f, delimiter="\t")
+        writer.writerow(["miRNA", "target", f"{metric}_score", "edge_weight"])
+        for mirna, _counts, score in selected:
+            targets = set()
+            targets.update(target_map["raw"].get(mirna, set()))
+            targets.update(target_map["canonical"].get(canonical_mirna(mirna), set()))
+            for target in sorted(targets):
+                writer.writerow([mirna, target, f"{score:.6g}", f"{score:.6g}"])
+                target_stats[target]["miRNAs"].add(mirna)
+                target_stats[target]["score_sum"] += score
+
+    with open(summary_path, "w", newline="") as f:
+        writer = csv.writer(f, delimiter="\t")
+        writer.writerow(["target", "miRNA_count", f"{metric}_score_sum"])
+        for target, stats in sorted(
+            target_stats.items(), key=lambda item: (-len(item[1]["miRNAs"]), item[0])
+        ):
+            writer.writerow([target, len(stats["miRNAs"]), f"{stats['score_sum']:.6g}"])
+
+    return edges_path, summary_path
+
+
+def main():
+    parser = argparse.ArgumentParser(description="Select high-expression miRNAs and build a target map.")
+    parser.add_argument("--counts", default="miRNA_counts.xlsx", help="Path to miRNA counts .xlsx file.")
+    parser.add_argument("--sheet", default="Mature", help="Sheet name to use (e.g., Mature or Hairpin).")
+    parser.add_argument("--metric", choices=["mean", "median", "sum", "max"], default="mean")
+    parser.add_argument("--top", type=int, default=None, help="Select top N miRNAs by metric.")
+    parser.add_argument(
+        "--min-score",
+        type=float,
+        default=5000.0,
+        help="Select miRNAs with metric >= value. Default: 5000.",
+    )
+    parser.add_argument("--quantile", type=float, default=None, help="Select miRNAs with metric >= quantile (0-1).")
+    parser.add_argument("--targets", default=None, help="TSV/CSV file with miRNA-to-target mappings.")
+    parser.add_argument("--mirtarbase", default=None, help="miRTarBase TSV/CSV/XLSX file to build targets.")
+    parser.add_argument(
+        "--fetch-mirtarbase",
+        action="store_true",
+        help="Auto-download miRTarBase MTI file if missing (best-effort).",
+    )
+    parser.add_argument(
+        "--mirtarbase-url",
+        default=None,
+        help="Direct URL to miRTarBase MTI file (overrides auto-detection).",
+    )
+    parser.add_argument(
+        "--mirtarbase-download-page",
+        default="http://mirtarbase.cuhk.edu.cn/php/download.php",
+        help="miRTarBase download page URL to discover MTI file.",
+    )
+    parser.add_argument(
+        "--mirtarbase-insecure",
+        action="store_true",
+        help="Allow insecure SSL (disable certificate verification) for downloads.",
+    )
+    parser.add_argument(
+        "--mirtarbase-species",
+        default="Homo sapiens",
+        help="Species filter for miRTarBase (default: Homo sapiens).",
+    )
+    parser.add_argument(
+        "--mirtarbase-support",
+        default=None,
+        help="Optional miRTarBase support filter (e.g., Strong).",
+    )
+    parser.add_argument(
+        "--mirtarbase-sheet",
+        default=None,
+        help="Sheet name for miRTarBase .xlsx (defaults to first sheet).",
+    )
+    parser.add_argument("--outdir", default="out", help="Output directory.")
+    args = parser.parse_args()
+
+    if not os.path.exists(args.counts):
+        print(f"Counts file not found: {args.counts}", file=sys.stderr)
+        return 2
+
+    try:
+        rows = read_xlsx_sheet(args.counts, args.sheet)
+    except Exception as exc:
+        print(f"Failed to read sheet: {exc}", file=sys.stderr)
+        return 2
+
+    if not rows:
+        print("No data found in sheet.", file=sys.stderr)
+        return 2
+
+    sample_cols, selected = select_high_expression(
+        rows, args.metric, args.top, args.min_score, args.quantile
+    )
+    if not selected:
+        print("No miRNAs selected with the current criteria.", file=sys.stderr)
+        return 2
+
+    outdir = args.outdir
+    high_path = os.path.join(outdir, "high_expression_miRNAs.tsv")
+    write_high_expression(high_path, sample_cols, selected, args.metric)
+
+    edges_path = summary_path = None
+    target_map = None
+    if args.mirtarbase or args.fetch_mirtarbase:
+        if args.mirtarbase is None:
+            args.mirtarbase = os.path.join(args.outdir, "miRTarBase_MTI.txt")
+
+        if not os.path.exists(args.mirtarbase) and args.fetch_mirtarbase:
+            try:
+                url = args.mirtarbase_url
+                if url is None:
+                    url = find_mirtarbase_mti_url(args.mirtarbase_download_page)
+                if url is None:
+                    print(
+                        "Could not auto-detect miRTarBase MTI file URL from the download page.",
+                        file=sys.stderr,
+                    )
+                    return 2
+                print(f"Downloading miRTarBase MTI file from: {url}")
+                try:
+                    download_url_to_file(url, args.mirtarbase, allow_insecure=args.mirtarbase_insecure)
+                except urllib.error.URLError as exc:
+                    reason = getattr(exc, "reason", None)
+                    if isinstance(reason, ssl.SSLError) and url.startswith("https://"):
+                        http_url = "http://" + url[len("https://") :]
+                        print(f"SSL failed, retrying over HTTP: {http_url}")
+                        download_url_to_file(
+                            http_url, args.mirtarbase, allow_insecure=args.mirtarbase_insecure
+                        )
+                    else:
+                        raise
+            except (urllib.error.URLError, ValueError, OSError) as exc:
+                print(f"Failed to download miRTarBase file: {exc}", file=sys.stderr)
+                return 2
+
+        if not os.path.exists(args.mirtarbase):
+            print(f"miRTarBase file not found: {args.mirtarbase}", file=sys.stderr)
+            return 2
+
+        try:
+            target_map = read_mirtarbase(
+                args.mirtarbase,
+                species=args.mirtarbase_species,
+                support=args.mirtarbase_support,
+                sheet=args.mirtarbase_sheet,
+                include_mirnas=[m for m, _c, _s in selected],
+            )
+        except Exception as exc:
+            print(f"Failed to read miRTarBase file: {exc}", file=sys.stderr)
+            return 2
+    elif args.targets:
+        if not os.path.exists(args.targets):
+            print(f"Targets file not found: {args.targets}", file=sys.stderr)
+            return 2
+        target_map = read_targets(args.targets)
+
+    if target_map:
+        edges_path, summary_path = write_targets(outdir, selected, target_map, args.metric)
+
+    print(f"Wrote high-expression list: {high_path}")
+    if edges_path:
+        print(f"Wrote target edges: {edges_path}")
+        print(f"Wrote target summary: {summary_path}")
+    else:
+        print("No target map generated (provide --mirtarbase or --targets).")
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())

Pilot/code/mirna_threshold_analysis.py

@@ -0,0 +1,256 @@
+"""
+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()

Pilot/code/mirna_threshold_unified.py

@@ -0,0 +1,370 @@
+"""
+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()

Pilot/code/qc_threshold_analysis.py

@@ -0,0 +1,62 @@
+import pandas as pd
+import numpy as np
+
+QC_FILE = "BEA25P077_RP.xlsx"
+PATIENTS_FILE = "patients.xlsx"
+THRESHOLDS = [0, 0.1, 0.25, 0.5, 0.75, 1.0, 1.5, 2.0, 5.0]
+
+# Load patient 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)
+
+def parse_conc(val):
+    """Return float concentration, or NaN for 'Too low' / missing."""
+    if isinstance(val, str):
+        return np.nan
+    try:
+        return float(val)
+    except (TypeError, ValueError):
+        return np.nan
+
+def analyse_sheet(sheet_name, conc_col):
+    df = pd.read_excel(QC_FILE, sheet_name=sheet_name, usecols=["ID", conc_col])
+    df["conc_val"] = df[conc_col].apply(parse_conc)
+    # Extract patient number from ID like "4A_G2" or "4B_G2"
+    extracted = df["ID"].str.extract(r"^(\d+)[AB]")[0]
+    df["Number"] = pd.to_numeric(extracted, errors="coerce")
+    df = df.dropna(subset=["Number"])
+    df["Number"] = df["Number"].astype(int)
+    # Merge with patient diagnosis
+    df = df.merge(patients, on="Number", how="left")
+    return df
+
+sheets = {
+    "A (Schirmer strips)": ("A-samples", "conc (ng/ul)"),
+    "B (Lens tissues)":    ("B-samples", "Conc [ng/ul]"),
+}
+
+print(f"{'Threshold':>10}  {'Sheet':<22}  {'Total':>6}  {'Healthy (Exf=0)':>15}  {'Diseased (Exf=1)':>16}")
+print("-" * 80)
+
+for label, (sheet_name, conc_col) in sheets.items():
+    df = analyse_sheet(sheet_name, conc_col)
+    total_samples = len(df)
+    n_too_low = df["conc_val"].isna().sum()
+    print(f"\n  {label}  ({total_samples} samples total, {n_too_low} 'Too low')")
+    print(f"  {'Threshold':>10}  {'Passing':>7}  {'Healthy':>8}  {'Diseased':>9}  {'Unknown dx':>10}")
+    print(f"  {'-'*55}")
+    for t in THRESHOLDS:
+        if t == 0:
+            # threshold=0: include all samples (Too low counts as passing)
+            passing = df.copy()
+        else:
+            # Only samples with a numeric conc >= threshold pass
+            passing = df[df["conc_val"] >= t]
+        healthy  = (passing["Exf"] == 0).sum()
+        diseased = (passing["Exf"] == 1).sum()
+        unknown  = passing["Exf"].isna().sum()
+        print(f"  {t:>10.2f}  {len(passing):>7}  {healthy:>8}  {diseased:>9}  {unknown:>10}")
+
+print()