diff --git a/app.py b/app.py index 3e3b047..f30589c 100644 --- a/app.py +++ b/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="โฌ‡๏ธ"), ], diff --git a/content/results_proteins.py b/content/results_proteins.py new file mode 100644 index 0000000..93d5532 --- /dev/null +++ b/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), + }, +) diff --git a/src/common/insight_tables.py b/src/common/insight_tables.py index 14cdcad..776a35e 100644 --- a/src/common/insight_tables.py +++ b/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) diff --git a/src/common/protein_table.py b/src/common/protein_table.py new file mode 100644 index 0000000..f8842bd --- /dev/null +++ b/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'{residue}' if mask[i] else residue) + blocks.append("".join(chars)) + lines.append(f'{row_start + 1:>5} ' + " ".join(blocks)) + return ( + "" + '
' + "\n".join(lines) + "
" + ) + diff --git a/tests/test_insight_table.py b/tests/test_insight_table.py index 6730d01..ea418ba 100644 --- a/tests/test_insight_table.py +++ b/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 diff --git a/tests/test_protein_table.py b/tests/test_protein_table.py new file mode 100644 index 0000000..5df4102 --- /dev/null +++ b/tests/test_protein_table.py @@ -0,0 +1,71 @@ +import numpy as np +import pandas as pd + +from src.common.protein_table import ( + build_tables, + coverage_html, + coverage_mask, + gene_name, + protein_name, + read_fasta, +) + + +def test_coverage_mask_counts_every_occurrence(): + mask = coverage_mask("PEPTIDEKAAAPEPTIDEK", ["PEPTIDEK", "AA"]) + assert mask.tolist() == [True] * 8 + [True, True, True] + [True] * 8 + assert coverage_mask("ACDEFG", ["XYZ"]).sum() == 0 + + +def test_uniprot_description_fields(): + desc = "Serum albumin OS=Homo sapiens OX=9606 GN=ALB PE=1 SV=2" + assert gene_name(desc) == "ALB" + assert protein_name(desc) == "Serum albumin" + assert gene_name("no gene here") == "" + + +def test_read_fasta_skips_decoys(tmp_path): + fasta = tmp_path / "db.fasta" + fasta.write_text(">sp|P1|A_HUMAN Protein A GN=A\nPEPTIDEKAAA\n>rev_sp|P1|A_HUMAN\nAAAKEDITPEP\n") + entries = read_fasta(fasta) + assert list(entries) == ["sp|P1|A_HUMAN"] + assert entries["sp|P1|A_HUMAN"][1] == "PEPTIDEKAAA" + + +def test_build_tables(): + fasta = { + "P1": ("Protein one GN=ONE", "PEPTIDEKSAMPLERAAAA"), + "P2": ("Protein two GN=TWO", "SAMPLERGGGG"), + } + psms = pd.DataFrame({ + "file": ["a.mzML", "a.mzML", "b.mzML"], + "scan": [1, 2, 3], + "rt": [10.0, 20.0, 30.0], + "mz": [400.0, 500.0, 400.0], + "charge": [2, 2, 3], + "peptide": ["PEPTIDEK", "SAMPLER", "PEPTIDEK"], + "sequence": ["PEPTIDEK", "SAMPLER", "PEPTIDEK"], + "score": [0.01, 0.02, 0.03], + "q_value": [0.001, 0.002, 0.003], + "proteins": [["P1"], ["P1", "P2"], ["P1"]], + }) + proteins, peptides, psms_long, coverage = build_tables(psms, fasta, ["a.mzML", "b.mzML"]) + + p1 = proteins.set_index("accession").loc["P1"] + assert p1["peptides"] == 2 + assert p1["unique peptides"] == 1 # SAMPLER is shared with P2 + assert p1["PSMs"] == 3 + assert p1["found in"] == "2 of 2" + assert p1["coverage (%)"] == round(100 * 15 / 19, 1) + assert p1["PSMs a"] == 2 and p1["PSMs b"] == 1 + assert proteins.iloc[0]["accession"] == "P1" # most PSMs first + + assert set(peptides.loc[peptides["accession"] == "P2", "peptide"]) == {"SAMPLER"} + assert len(psms_long) == 4 # the shared PSM appears under both proteins + assert coverage["P2"].tolist() == [True] * 7 + [False] * 4 + + +def test_coverage_html_marks_covered_residues(): + out = coverage_html("ACDEFGHIKLM", np.array([True, True] + [False] * 9)) + assert out.count('') == 2 + assert "ACDEFGHIKL M" in out.replace('A', "A").replace('C', "C")