Add a Proteins page: protein table with coverage, linked peptides and PSMs

app.py

@@ -29,6 +29,7 @@
             st.Page(Path("content", "results_database_search.py"), title="Database Search", icon="🔬"),
             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_proteins.py"), title="Proteins", icon="🧩"),
             st.Page(Path("content", "results_abundance.py"), title="Protein Abundance", icon="📋"),
             st.Page(Path("content", "results_download.py"), title="Download", icon="⬇️"),
         ],

content/results_proteins.py

@@ -0,0 +1,118 @@
+"""Protein Identification Results Page."""
+from pathlib import Path
+
+import streamlit as st
+from openms_insight import StateManager
+
+from src.common.common import page_setup
+from src.common.insight_tables import show_insight_table
+from src.common.protein_table import build_tables, collect_psms, coverage_html, read_fasta
+from src.common.results_helpers import get_workflow_dir
+from src.WorkflowTest import WorkflowTest
+
+params = page_setup()
+wf = WorkflowTest()
+
+st.title("Proteins")
+st.markdown(
+    """
+Proteins identified by the FDR-filtered PSMs of all runs, with their sequence
+coverage. **Select a protein** to list its peptides and see its coverage map;
+**select a peptide** to list its PSMs.
+"""
+)
+
+workflow_dir = get_workflow_dir(st.session_state["workspace"])
+filter_dir = workflow_dir / "results" / "psm_filter"
+cache_dir = workflow_dir / "results" / "insight_cache"
+idxml_files = sorted(filter_dir.glob("*.idXML")) if filter_dir.exists() else []
+
+if not idxml_files:
+    st.info("No filtered PSMs available yet. Please run the workflow first.")
+    st.page_link("content/workflow_run.py", label="Go to Run", icon="🚀")
+    st.stop()
+
+fasta_path = None
+if wf.params.get("fasta-file"):
+    try:
+        fasta_path = Path(wf.file_manager.get_files([wf.params["fasta-file"]])[0])
+    except (ValueError, IndexError):
+        fasta_path = None
+if fasta_path is None or not fasta_path.is_file():
+    st.warning("The FASTA database of this run was not found, so coverage and descriptions cannot be shown.")
+
+mzml_names = [Path(f).name for f in wf.params.get("mzML-files", []) or []]
+source_files = [*idxml_files, *([fasta_path] if fasta_path and fasta_path.is_file() else [])]
+
+
+@st.cache_data(show_spinner="Building protein tables...")
+def load_tables(paths: tuple[str, ...], fasta: str | None, mtimes: tuple[float, ...], mzml_names: tuple[str, ...]):
+    # mtimes only key the cache, so a re-run of the workflow invalidates it
+    psms = collect_psms([Path(p) for p in paths], list(mzml_names))
+    fasta_entries = read_fasta(fasta) if fasta else {}
+    files = sorted(psms["file"].unique())
+    return (*build_tables(psms, fasta_entries, files), fasta_entries)
+
+
+proteins, peptides, psms_long, coverage, fasta_entries = load_tables(
+    tuple(str(p) for p in idxml_files),
+    str(fasta_path) if fasta_path and fasta_path.is_file() else None,
+    tuple(p.stat().st_mtime for p in source_files),
+    tuple(mzml_names),
+)
+
+if proteins.empty:
+    st.info("No target proteins were identified.")
+    st.stop()
+
+c1, c2, c3 = st.columns(3)
+c1.metric("Proteins", f"{len(proteins):,}")
+c2.metric("Peptides", f"{peptides['peptide'].nunique():,}")
+c3.metric("Target PSMs", f"{psms_long.drop_duplicates(['file', 'scan', 'peptide']).shape[0]:,}")
+
+state_manager = StateManager("proteins_page")
+# Columns the linked tables filter on: already shown by the selection above them
+HIDDEN = {"visible": False}
+
+
+def decimals(n: int) -> dict:
+    return {"formatter": "money", "formatterParams": {"thousand": False, "precision": n}}
+
+
+st.subheader("Proteins")
+show_insight_table(
+    proteins, name="proteins", cache_dir=cache_dir, source_files=source_files,
+    title="proteins", height=450, state_manager=state_manager,
+    interactivity={"protein": "accession"},
+    column_overrides={"coverage (%)": decimals(1)},
+)
+
+selected = state_manager.get_selection("protein")
+if selected in coverage and selected in fasta_entries:
+    description, sequence = fasta_entries[selected]
+    mask = coverage[selected]
+    st.subheader(f"Coverage of {selected}")
+    st.caption(
+        f"{description} · {len(sequence)} residues · {mask.sum()} covered ({100 * mask.mean():.1f} %). "
+        "Covered residues are highlighted; every occurrence of a peptide counts."
+    )
+    st.html(coverage_html(sequence, mask))
+
+st.subheader("Peptides of the selected protein")
+show_insight_table(
+    peptides, name="protein_peptides", cache_dir=cache_dir, source_files=source_files,
+    title="peptides", height=350, state_manager=state_manager,
+    filters={"protein": "accession"}, interactivity={"peptide": "peptide"},
+    column_overrides={"accession": HIDDEN},
+)
+
+st.subheader("PSMs of the selected peptide")
+show_insight_table(
+    psms_long, name="protein_psms", cache_dir=cache_dir, source_files=source_files,
+    title="psms", height=350, state_manager=state_manager,
+    filters={"protein": "accession", "peptide": "peptide"},
+    column_overrides={
+        "accession": HIDDEN, "peptide": HIDDEN,
+        "RT (s)": decimals(1), "m/z": decimals(4),
+    },
+)

src/common/insight_tables.py

@@ -54,7 +54,7 @@ def _cache_key(source_files, extra: str) -> str:
     return hashlib.sha1("|".join(parts).encode()).hexdigest()[:10]
 
 
-def _to_polars(df: pd.DataFrame) -> tuple[pl.DataFrame, list[dict]]:
+def _to_polars(df: pd.DataFrame, column_overrides: dict | None = None) -> tuple[pl.DataFrame, list[dict]]:
     """Convert to polars with Tabulator-safe field names and column definitions.
 
     Tabulator reads a "." in a field name as a nested lookup, so such columns get
@@ -73,6 +73,7 @@ def _to_polars(df: pd.DataFrame) -> tuple[pl.DataFrame, list[dict]]:
         if pd.api.types.is_numeric_dtype(df[column]):
             col_def["sorter"] = "number"
             col_def["hozAlign"] = "right"
+        col_def.update((column_overrides or {}).get(column, {}))
         column_definitions.append(col_def)
     df = df.rename(columns=renames)
     for column in df.columns:
@@ -91,6 +92,8 @@ def show_insight_table(
     title: str | None = None,
     height: int = 500,
     page_size: int = 100,
+    state_manager=None,
+    column_overrides: dict | None = None,
     **table_kwargs,
 ) -> None:
     """Render ``df`` as a server-side paginated openms-insight table.
@@ -115,6 +118,11 @@ def show_insight_table(
         title: Title, also used as the CSV download's file name.
         height: Component height in pixels.
         page_size: Rows per page.
+        state_manager: Shared ``StateManager`` for tables linked through
+            ``filters``/``interactivity``; defaults to openms-insight's own.
+        column_overrides: Extra Tabulator column settings per column name,
+            e.g. a ``formatter``. openms-insight stores floats as float32, so
+            a column of rounded values needs one to display without noise.
         **table_kwargs: Passed on to ``Table`` (e.g. ``initial_sort``).
     """
     cache_dir = Path(cache_dir) if cache_dir else _default_cache_dir()
@@ -127,7 +135,7 @@ def show_insight_table(
         for stale in cache_dir.glob(f"{name}_" + "?" * 10):
             shutil.rmtree(stale, ignore_errors=True)
         cache_dir.mkdir(parents=True, exist_ok=True)
-        pl_df, column_definitions = _to_polars(df)
+        pl_df, column_definitions = _to_polars(df, column_overrides)
         Table(
             cache_id=cache_id,
             data=pl_df.lazy(),
@@ -140,4 +148,4 @@ def show_insight_table(
         )
 
     table = Table(cache_id=cache_id, cache_path=str(cache_dir))
-    table(key=f"insight_table_{name}", height=height)
+    table(key=f"insight_table_{name}", height=height, state_manager=state_manager)

src/common/protein_table.py

@@ -0,0 +1,175 @@
+"""Protein-level identification tables: proteins, their peptides and PSMs, and sequence coverage.
+
+The tables are built from the FDR-filtered PSMs of every run and the target
+FASTA database, so they exist for LFQ and TMT alike. Decoy hits are left out.
+"""
+import html
+import re
+from pathlib import Path
+
+import numpy as np
+import pandas as pd
+from pyopenms import FASTAFile
+
+from src.common.results_helpers import extract_filename_from_idxml, extract_scan_from_ref, load_idxml
+
+DECOY_PREFIXES = ("DECOY_", "decoy_", "rev_", "REV_", "XXX_", "reversed_")
+Q_VALUE_KEYS = ("MS:1001491", "q-value")
+
+PSM_COLUMNS = ["file", "scan", "rt", "mz", "charge", "peptide", "sequence", "score", "q_value", "proteins"]
+
+
+def _is_decoy_accession(accession: str) -> bool:
+    return accession.startswith(DECOY_PREFIXES)
+
+
+def read_fasta(path) -> dict[str, tuple[str, str]]:
+    """Accession -> (description, sequence) of the target entries of a FASTA file."""
+    entries = []
+    FASTAFile().load(str(path), entries)
+    return {
+        e.identifier: (e.description, e.sequence)
+        for e in entries
+        if not _is_decoy_accession(e.identifier)
+    }
+
+
+def gene_name(description: str) -> str:
+    """The gene from a UniProt description (``GN=``), or an empty string."""
+    match = re.search(r"\bGN=(\S+)", description or "")
+    return match.group(1) if match else ""
+
+
+def protein_name(description: str) -> str:
+    """A UniProt description without its ``OS=``/``OX=``/``GN=``/``PE=``/``SV=`` fields."""
+    return re.split(r"\s+[A-Z]{2}=", description or "", maxsplit=1)[0]
+
+
+def collect_psms(idxml_paths, mzml_names=None) -> pd.DataFrame:
+    """Top hits of all runs, one row per target PSM, with their target proteins."""
+    records = []
+    for path in idxml_paths:
+        file = extract_filename_from_idxml(Path(path), mzml_names)
+        for pep in load_idxml(path)[1]:
+            hits = pep.getHits()
+            if not hits:
+                continue
+            hit = hits[0]
+            proteins = sorted({
+                ev.getProteinAccession() for ev in hit.getPeptideEvidences()
+                if not _is_decoy_accession(ev.getProteinAccession())
+            })
+            if not proteins:
+                continue  # a decoy hit
+            ref = pep.getMetaValue("spectrum_reference") if pep.metaValueExists("spectrum_reference") else ""
+            if isinstance(ref, bytes):
+                ref = ref.decode()
+            q_value = np.nan
+            for key in Q_VALUE_KEYS:
+                if hit.metaValueExists(key):
+                    q_value = float(hit.getMetaValue(key))
+                    break
+            seq = hit.getSequence()
+            records.append({
+                "file": file,
+                "scan": extract_scan_from_ref(ref),
+                "rt": pep.getRT(),
+                "mz": pep.getMZ(),
+                "charge": hit.getCharge(),
+                "peptide": seq.toString(),
+                "sequence": seq.toUnmodifiedString(),
+                "score": hit.getScore(),
+                "q_value": q_value,
+                "proteins": proteins,
+            })
+    return pd.DataFrame(records, columns=PSM_COLUMNS)
+
+
+def coverage_mask(protein_sequence: str, peptides) -> np.ndarray:
+    """Per residue, whether some peptide covers it (every occurrence of each peptide counts)."""
+    mask = np.zeros(len(protein_sequence), dtype=bool)
+    for peptide in set(peptides):
+        if not peptide:
+            continue
+        start = protein_sequence.find(peptide)
+        while start != -1:
+            mask[start:start + len(peptide)] = True
+            start = protein_sequence.find(peptide, start + 1)
+    return mask
+
+
+def build_tables(psms: pd.DataFrame, fasta: dict[str, tuple[str, str]], files: list[str]):
+    """Proteins, peptide-to-protein and PSM-to-protein tables.
+
+    Peptides and PSMs are in long format (one row per protein they map to), so
+    selecting a protein filters them by its accession.
+
+    Returns:
+        (proteins, peptides, psms_long, coverage) where ``coverage`` maps an
+        accession to its residue mask.
+    """
+    long = psms.explode("proteins").rename(columns={"proteins": "accession"})
+    long["shared"] = psms["proteins"].str.len().reindex(long.index).gt(1)
+
+    rows, coverage = [], {}
+    for accession, group in long.groupby("accession", sort=False):
+        description, sequence = fasta.get(accession, ("", ""))
+        mask = coverage_mask(sequence, group["sequence"])
+        coverage[accession] = mask
+        per_file = group.groupby("file").size()
+        row = {
+            "accession": accession,
+            "gene": gene_name(description),
+            "description": protein_name(description),
+            "coverage (%)": round(100 * mask.mean(), 1) if len(mask) else np.nan,
+            "peptides": group["sequence"].nunique(),
+            "unique peptides": group.loc[~group["shared"], "sequence"].nunique(),
+            "PSMs": len(group),
+            "found in": f"{per_file.size} of {len(files)}",
+            "length": len(sequence) or np.nan,
+        }
+        for file in files:
+            row[f"PSMs {Path(file).stem}"] = int(per_file.get(file, 0))
+        rows.append(row)
+    proteins = pd.DataFrame(rows)
+    if not proteins.empty:
+        proteins = proteins.sort_values(["PSMs", "peptides"], ascending=False, ignore_index=True)
+        # the description last: it is the widest column, and the one to truncate
+        proteins = proteins[[c for c in proteins.columns if c != "description"] + ["description"]]
+        if not proteins["gene"].astype(bool).any():
+            proteins = proteins.drop(columns="gene")
+
+    peptides = (
+        long.groupby(["accession", "peptide"], sort=False)
+        .agg(PSMs=("file", "size"), files=("file", "nunique"),
+             best_q=("q_value", "min"), charges=("charge", lambda c: ",".join(map(str, sorted(set(c))))),
+             shared=("shared", "first"))
+        .reset_index()
+        .rename(columns={"best_q": "best q-value", "files": "found in files"})
+        .sort_values(["accession", "PSMs"], ascending=[True, False], ignore_index=True)
+    )
+    psms_long = long[["accession", "peptide", "file", "scan", "rt", "mz", "charge", "score", "q_value"]].rename(
+        columns={"rt": "RT (s)", "mz": "m/z", "q_value": "q-value"}
+    ).reset_index(drop=True)
+    return proteins, peptides, psms_long, coverage
+
+
+def coverage_html(sequence: str, mask: np.ndarray, block: int = 10, per_row: int = 50) -> str:
+    """The protein sequence in blocks of ``block`` residues, covered residues highlighted."""
+    lines = []
+    for row_start in range(0, len(sequence), per_row):
+        blocks = []
+        for start in range(row_start, min(row_start + per_row, len(sequence)), block):
+            chars = []
+            for i in range(start, min(start + block, len(sequence))):
+                residue = html.escape(sequence[i])
+                chars.append(f'<span class="cov">{residue}</span>' if mask[i] else residue)
+            blocks.append("".join(chars))
+        lines.append(f'<span class="pos">{row_start + 1:>5}</span> ' + " ".join(blocks))
+    return (
+        "<style>.covmap{font-family:monospace;font-size:14px;line-height:1.7;white-space:pre}"
+        ".covmap .cov{background:#d7f5ea;color:#065f46;font-weight:600}"
+        ".covmap .pos{color:#a1a1aa}</style>"
+        '<div class="covmap">' + "\n".join(lines) + "</div>"
+    )
+

tests/test_insight_table.py

@@ -99,3 +99,12 @@ def test_table_without_source_files_is_cached_by_content(tmp_path, monkeypatch):
     assert rendered[2] != rendered[0]
     # The cache of the previous version of the table is removed
     assert [p.name for p in tmp_path.iterdir()] == [rendered[2]]
+
+
+def test_to_polars_applies_column_overrides():
+    df = pd.DataFrame({"a.b": [1.25], "c": ["x"]})
+    _, col_defs = _to_polars(df, {"a.b": {"formatter": "money"}, "c": {"visible": False}})
+    by_title = {d["title"]: d for d in col_defs}
+    assert by_title["a.b"]["formatter"] == "money"
+    assert by_title["a.b"]["field"] == "a_b"
+    assert by_title["c"]["visible"] is False