From aed714c18967ee9236e9e40c2f02ec4fb3dd9f19 Mon Sep 17 00:00:00 2001 From: Claude Date: Sun, 25 Oct 2026 06:04:18 +0000 Subject: [PATCH] Add a Samples page: PSMs, peptides and proteins per sample, by condition Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01D5Wo2FHPpi5a8X6GppNxEk --- app.py | 1 + content/results_samples.py | 142 +++++++++++++++++++++++++++++++++++ src/common/samples.py | 147 +++++++++++++++++++++++++++++++++++++ tests/test_samples.py | 108 +++++++++++++++++++++++++++ 4 files changed, 398 insertions(+) create mode 100644 content/results_samples.py create mode 100644 src/common/samples.py create mode 100644 tests/test_samples.py diff --git a/app.py b/app.py index 3e3b047..c1abf04 100644 --- a/app.py +++ b/app.py @@ -30,6 +30,7 @@ st.Page(Path("content", "results_rescoring.py"), title="Rescoring", icon="๐Ÿ“ˆ"), st.Page(Path("content", "results_filtered.py"), title="Filtered PSMs", icon="๐ŸŽฏ"), st.Page(Path("content", "results_abundance.py"), title="Protein Abundance", icon="๐Ÿ“‹"), + st.Page(Path("content", "results_samples.py"), title="Samples", icon="๐Ÿงช"), st.Page(Path("content", "results_download.py"), title="Download", icon="โฌ‡๏ธ"), ], # Steps run in order; each one works on the previous step's output. diff --git a/content/results_samples.py b/content/results_samples.py new file mode 100644 index 0000000..4c0ffd4 --- /dev/null +++ b/content/results_samples.py @@ -0,0 +1,142 @@ +"""Samples Page: PSMs, peptides and proteins per sample, grouped by condition.""" +from pathlib import Path + +import plotly.express as px +import streamlit as st + +from src.common.common import page_setup, show_fig +from src.common.insight_tables import show_insight_table +from src.common.results_helpers import get_abundance_data, get_sample_group_map, get_workflow_dir +from src.common.samples import COUNT_COLUMNS, channel_samples, condition_summary, lfq_samples +from src.workflow.ParameterManager import ParameterManager + +params = page_setup() +st.title("Samples") + +st.markdown( + """ +How many PSMs, peptides and proteins each sample contributed, grouped by the +conditions set on the Configure page. A sample far below the others in its +condition points to a failed run or a loading problem. +""" +) + +if "workspace" not in st.session_state: + st.warning("Please initialize your workspace first.") + st.stop() + +workflow_dir = get_workflow_dir(st.session_state["workspace"]) +results_dir = workflow_dir / "results" +cache_dir = results_dir / "insight_cache" +workflow_params = ParameterManager(workflow_dir, "TOPP Workflow").get_parameters_from_json() +analysis_mode = workflow_params.get("analysis-mode", "LFQ") +fdr = float(workflow_params.get("psm-fdr-percent", 1.0)) + + +def _mtimes(paths) -> tuple: + return tuple(Path(p).stat().st_mtime if Path(p).exists() else 0.0 for p in paths) + + +@st.cache_data(show_spinner="Counting per sample...") +def _lfq_samples(workflow_dir: str, files: tuple, groups: tuple, mtimes: tuple): + """lfq_samples, cached on the group assignments and the results' modification times.""" + return lfq_samples(workflow_dir, list(files), dict(groups)) + + +if analysis_mode == "LFQ": + filter_dir = results_dir / "psm_filter" + if not filter_dir.is_dir() or not any(filter_dir.glob("*.idXML")): + st.info("No results available yet. Please run the workflow first.") + st.page_link("content/workflow_run.py", label="Go to Run", icon="๐Ÿš€") + st.stop() + mzml_files = [str(f) for f in workflow_params.get("mzML-files", [])] + group_map = { + key[len("mzML-group-"):]: value + for key, value in workflow_params.items() + if key.startswith("mzML-group-") and value + } + sources = sorted(filter_dir.glob("*.idXML")) + sorted((results_dir / "quant_results").glob("*.csv")) + samples = _lfq_samples( + str(workflow_dir), tuple(mzml_files), tuple(sorted(group_map.items())), _mtimes(sources) + ) + st.caption( + f"PSMs, peptides and proteins are the targets identified at {fdr:g} % PSM FDR, from the " + "top hit of each spectrum. Quantified proteins and peptides have an intensity above 0 in " + "the sample's ProteomicsLFQ result." + ) +else: + abundance = get_abundance_data(st.session_state["workspace"]) + if abundance is None: + st.info("No quantification result available yet. Please run the workflow first.") + st.page_link("content/workflow_run.py", label="Go to Run", icon="๐Ÿš€") + st.stop() + pivot_df, _, group_map = abundance + sample_group_map = get_sample_group_map(st.session_state["workspace"], pivot_df, group_map) + # The report's first columns are the protein and its description; the channels follow + sample_columns = [c for c in pivot_df.columns[4:] if c != "PeptideSequence"] + samples = channel_samples(pivot_df, sample_columns, sample_group_map) + sources = sorted((results_dir / "quant_results").glob("*.csv")) + st.caption( + "Each sample is an isobaric channel. Quantified proteins have an intensity above 0 in the " + "channel. Identifications are made per file, not per channel: see the Overview page." + ) + +if samples.empty: + st.info("No samples in the last run.") + st.stop() + +def _table_height(rows: int) -> int: + """Room for every row up to a page of 12, plus the header and pager.""" + return 140 + 34 * min(rows, 12) + + +counts = [c for c in COUNT_COLUMNS if samples[c].notna().any()] +group_key = "|".join(f"{s}={c}" for s, c in zip(samples["Sample"], samples["Condition"])) + +st.subheader("Per sample") +metric = st.selectbox("Count", counts, key="samples-metric") +fig = px.bar( + samples, + x="Sample", + y=metric, + color="Condition", + category_orders={"Sample": samples["Sample"].tolist()}, +) +# One full-width bar per sample: px gives each condition its own offset group, +# which would leave room beside every bar for the other conditions +fig.update_traces(offsetgroup=None, alignmentgroup=None) +fig.update_layout(xaxis_title=None, legend_title_text="Condition", bargap=0.2) +# Slant the sample names once they no longer fit side by side +label_chars = sum(len(s) for s in samples["Sample"]) +fig.update_xaxes(tickangle=-35 if label_chars > 90 else 0, automargin=True) +show_fig(fig, f"samples_{metric.lower().replace(' ', '_')}") + +show_insight_table( + samples[["Sample", "Condition"] + counts], + name="samples_per_sample", + cache_dir=cache_dir, + source_files=sources, + extra_key=group_key, + title="samples", + height=_table_height(len(samples)), + layout="fitColumns", +) + +st.subheader("Per condition") +summary = condition_summary(samples) +# Means go as text, since the table stores floats as float32 and would show 36.8 as 36.79999923706055 +for col in summary.columns[2:]: + summary[col] = summary[col].map(lambda v: "" if v != v else f"{v:g}") +show_insight_table( + summary, + name="samples_per_condition", + cache_dir=cache_dir, + source_files=sources, + extra_key=group_key, + title="conditions", + height=_table_height(len(summary)), + layout="fitColumns", +) +if (samples["Condition"] == "Unassigned").any(): + st.caption("Samples without a condition are shown as Unassigned. Set conditions on the Configure page.") + st.page_link("content/workflow_configure.py", label="Go to Configure", icon="โš™๏ธ") diff --git a/src/common/samples.py b/src/common/samples.py new file mode 100644 index 0000000..9ae1cc8 --- /dev/null +++ b/src/common/samples.py @@ -0,0 +1,147 @@ +"""Per-sample counts of a finished run: PSMs, peptides and proteins, by condition. + +Reads only what the workflow already wrote: the IDFilter idXMLs (identifications +at the run's PSM FDR) and the quantification result. It makes no Streamlit +calls, so the page caches it and tests call it. +""" +from pathlib import Path + +import pandas as pd + +from src.common.results_helpers import extract_filename_from_idxml, load_idxml + +UNASSIGNED = "Unassigned" + +SAMPLE_COLUMNS = [ + "Sample", + "Condition", + "Target PSMs", + "Peptides", + "Proteins", + "Quantified proteins", + "Quantified peptides", +] + +# Columns the bar chart can show, in the order the page offers them +COUNT_COLUMNS = SAMPLE_COLUMNS[2:] + + +def _is_decoy(hit) -> bool: + if not hit.metaValueExists("target_decoy"): + return False + value = hit.getMetaValue("target_decoy") + if isinstance(value, bytes): + value = value.decode() + return str(value).startswith("decoy") + + +def identification_counts(idxml_path) -> dict: + """Target PSMs, distinct peptides and distinct proteins of one idXML. + + Only the top hit of each spectrum counts, and decoy hits are left out. + Proteins are the accessions the target peptides map to. + """ + _, peptides = load_idxml(idxml_path) + psms = 0 + sequences, proteins = set(), set() + for pep in peptides: + hits = pep.getHits() + if not hits or _is_decoy(hits[0]): + continue + hit = hits[0] + psms += 1 + sequences.add(hit.getSequence().toString()) + proteins.update(ev.getProteinAccession() for ev in hit.getPeptideEvidences()) + return {"Target PSMs": psms, "Peptides": len(sequences), "Proteins": len(proteins)} + + +def lfq_quantified_counts(quant_csv) -> dict: + """Quantified proteins and peptides per input file name, from ProteomicsLFQ's CSV. + + A protein (group) or peptide counts for a file when its intensity there is above 0. + """ + try: + df = pd.read_csv(quant_csv, usecols=["ProteinName", "PeptideSequence", "Reference", "Intensity"]) + except (OSError, ValueError): + return {} + df = df[pd.to_numeric(df["Intensity"], errors="coerce") > 0] + counts = df.groupby("Reference").agg( + proteins=("ProteinName", "nunique"), peptides=("PeptideSequence", "nunique") + ) + return { + Path(str(ref)).name: {"Quantified proteins": int(row.proteins), "Quantified peptides": int(row.peptides)} + for ref, row in counts.iterrows() + } + + +def _condition(group_map: dict, *names) -> str: + for name in names: + value = group_map.get(name) + if value is not None and str(value).strip(): + return str(value).strip() + return UNASSIGNED + + +def _order(df: pd.DataFrame) -> pd.DataFrame: + """Samples grouped by condition (unassigned last), in their original order within one.""" + conditions = sorted(c for c in df["Condition"].unique() if c != UNASSIGNED) + if UNASSIGNED in set(df["Condition"]): + conditions.append(UNASSIGNED) + df = df.copy() + df["Condition"] = pd.Categorical(df["Condition"], categories=conditions, ordered=True) + df = df.sort_values("Condition", kind="stable").reset_index(drop=True) + df["Condition"] = df["Condition"].astype(str) + return df + + +def lfq_samples(workflow_dir, mzml_files, group_map: dict) -> pd.DataFrame: + """One row per input file (SAMPLE_COLUMNS); files the run did not reach show empty cells. + + ``group_map`` maps an input file name to its condition, as set on the Configure page. + """ + results = Path(workflow_dir) / "results" + names = [Path(f).name for f in mzml_files] + filter_dir = results / "psm_filter" + idxmls = {} + if filter_dir.is_dir(): + idxmls = {extract_filename_from_idxml(p, names): p for p in sorted(filter_dir.glob("*.idXML"))} + quant_csvs = sorted((results / "quant_results").glob("*.csv")) if (results / "quant_results").is_dir() else [] + quantified = lfq_quantified_counts(quant_csvs[0]) if quant_csvs else {} + + rows = [] + for name in names: + ids = identification_counts(idxmls[name]) if name in idxmls else {} + rows.append({ + "Sample": Path(name).stem, + "Condition": _condition(group_map, name, Path(name).stem), + **{col: ids.get(col) for col in ("Target PSMs", "Peptides", "Proteins")}, + **{col: quantified.get(name, {}).get(col) for col in ("Quantified proteins", "Quantified peptides")}, + }) + return _order(pd.DataFrame(rows, columns=SAMPLE_COLUMNS)) + + +def channel_samples(pivot_df: pd.DataFrame, sample_columns, sample_group_map: dict) -> pd.DataFrame: + """One row per isobaric channel: the proteins quantified in it (intensity above 0). + + Identifications are per file, not per channel, so those columns stay empty. + """ + rows = [] + for col in sample_columns: + values = pd.to_numeric(pivot_df[col], errors="coerce") + rows.append({ + "Sample": col, + "Condition": _condition(sample_group_map, col), + "Quantified proteins": int((values > 0).sum()), + }) + return _order(pd.DataFrame(rows, columns=SAMPLE_COLUMNS)) + + +def condition_summary(samples: pd.DataFrame) -> pd.DataFrame: + """Per condition: the number of samples and the mean of each count over them.""" + counts = [c for c in COUNT_COLUMNS if samples[c].notna().any()] + numeric = samples[["Condition"] + counts].copy() + numeric[counts] = numeric[counts].apply(pd.to_numeric, errors="coerce") + summary = numeric.groupby("Condition", sort=False).agg( + Samples=("Condition", "size"), **{f"Mean {c}": (c, "mean") for c in counts} + ) + return summary.round(1).reset_index() diff --git a/tests/test_samples.py b/tests/test_samples.py new file mode 100644 index 0000000..8270312 --- /dev/null +++ b/tests/test_samples.py @@ -0,0 +1,108 @@ +"""Tests for the per-sample counts (Samples page).""" +import pandas as pd +import pyopenms as poms + +# src.common.samples is imported inside the tests, not at collection: importing +# it pulls in src.workflow.ParameterManager with the real streamlit, and +# test_parameter_defaults needs to import it with a mock. + + +def _idxml(path, hits): + """hits: (sequence, protein accession, "target" or "decoy").""" + prot = poms.ProteinIdentification() + prot.setIdentifier("run") + peptides = [] + for seq, accession, td in hits: + pep = poms.PeptideIdentification() + pep.setIdentifier("run") + hit = poms.PeptideHit() + hit.setSequence(poms.AASequence.fromString(seq)) + hit.setCharge(2) + hit.setMetaValue("target_decoy", td) + evidence = poms.PeptideEvidence() + evidence.setProteinAccession(accession) + hit.setPeptideEvidences([evidence]) + pep.setHits([hit]) + peptides.append(pep) + protein_hits = [] + for accession in sorted({h[1] for h in hits}): + protein_hit = poms.ProteinHit() + protein_hit.setAccession(accession) + protein_hits.append(protein_hit) + prot.setHits(protein_hits) + if hasattr(poms, "PeptideIdentificationList"): + container = poms.PeptideIdentificationList() + for pep in peptides: + container.push_back(pep) + peptides = container + path.parent.mkdir(parents=True, exist_ok=True) + poms.IdXMLFile().store(str(path), [prot], peptides) + + +def _quant_csv(path, rows): + """rows: (protein, peptide, reference, intensity).""" + path.parent.mkdir(parents=True, exist_ok=True) + pd.DataFrame(rows, columns=["ProteinName", "PeptideSequence", "Reference", "Intensity"]).to_csv(path, index=False) + + +def test_identification_counts_skip_decoys(tmp_path): + from src.common.samples import identification_counts + + _idxml(tmp_path / "a.idXML", [ + ("PEPTIDE", "P1", "target"), + ("PEPTIDE", "P1", "target"), + ("PEPTIDEK", "P2", "target"), + ("EDITPEP", "DECOY_P1", "decoy"), + ]) + assert identification_counts(tmp_path / "a.idXML") == {"Target PSMs": 3, "Peptides": 2, "Proteins": 2} + + +def test_lfq_samples_group_by_condition(tmp_path): + from src.common.samples import lfq_samples + + results = tmp_path / "results" + _idxml(results / "psm_filter" / "ctrl1_filter.idXML", [("PEPTIDE", "P1", "target")]) + _idxml(results / "psm_filter" / "case1_filter.idXML", [ + ("PEPTIDE", "P1", "target"), ("PEPTIDEK", "P2", "target"), + ]) + _idxml(results / "psm_filter" / "case2_filter.idXML", [("PEPTIDE", "P1", "target")]) + _quant_csv(results / "quant_results" / "quant.csv", [ + ("P1", "PEPTIDE", "ctrl1.mzML", 10.0), + ("P1", "PEPTIDE", "case1.mzML", 5.0), + ("P2", "PEPTIDEK", "case1.mzML", 3.0), + ("P2", "PEPTIDEK", "case2.mzML", 0.0), + ]) + files = ["in/ctrl1.mzML", "in/case1.mzML", "in/case2.mzML", "in/blank.mzML"] + groups = {"ctrl1.mzML": "control", "case1.mzML": "case", "case2.mzML": "case"} + + df = lfq_samples(tmp_path, files, groups) + + # Conditions sorted, unassigned last; file order kept within a condition + assert df["Sample"].tolist() == ["case1", "case2", "ctrl1", "blank"] + assert df["Condition"].tolist() == ["case", "case", "control", "Unassigned"] + case1 = df.iloc[0] + assert (case1["Target PSMs"], case1["Peptides"], case1["Proteins"]) == (2, 2, 2) + assert (case1["Quantified proteins"], case1["Quantified peptides"]) == (2, 2) + # An intensity of 0 is not quantified + assert pd.isna(df.iloc[1]["Quantified proteins"]) + # A file the run did not reach has empty cells + assert df.iloc[3][["Target PSMs", "Quantified proteins"]].isna().all() + + +def test_channel_samples_and_condition_summary(): + from src.common.samples import channel_samples, condition_summary + + pivot = pd.DataFrame({ + "protein": ["P1", "P2", "P3"], + "sample1[126]": [1.0, 2.0, 0.0], + "sample2[127N]": [1.0, None, None], + "sample3[127C]": [1.0, 1.0, 1.0], + }) + groups = {"sample1[126]": "A", "sample2[127N]": "A", "sample3[127C]": "B"} + df = channel_samples(pivot, ["sample1[126]", "sample2[127N]", "sample3[127C]"], groups) + assert df["Quantified proteins"].tolist() == [2, 1, 3] + assert df["Target PSMs"].isna().all() + + summary = condition_summary(df) + assert summary.columns.tolist() == ["Condition", "Samples", "Mean Quantified proteins"] + assert summary.set_index("Condition")["Mean Quantified proteins"].to_dict() == {"A": 1.5, "B": 3.0}