From 00fd9723cd160dc4bed9e55decba75c95a8768a3 Mon Sep 17 00:00:00 2001 From: Claude Date: Wed, 21 Oct 2026 06:00:25 +0000 Subject: [PATCH] Filtered PSMs: fragment mass error plot and b/y ion table for the selected PSM Third step towards feature parity with the ProSE web app's PSM view. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01D5Wo2FHPpi5a8X6GppNxEk --- content/results_filtered.py | 32 ++++++++ src/common/fragment_plots.py | 139 +++++++++++++++++++++++++++++++++++ tests/test_fragment_plots.py | 63 ++++++++++++++++ 3 files changed, 234 insertions(+) create mode 100644 src/common/fragment_plots.py create mode 100644 tests/test_fragment_plots.py diff --git a/content/results_filtered.py b/content/results_filtered.py index d6f31aa..347fcac 100644 --- a/content/results_filtered.py +++ b/content/results_filtered.py @@ -4,6 +4,9 @@ from src.common.common import page_setup from src.common.results_helpers import get_workflow_dir from openms_insight import Table, Heatmap, LinePlot, SequenceView, StateManager +from src.common.common import show_fig +from src.common.fragment_plots import fragment_error_plot, ion_table, ion_table_html, load_psm_fragments +from src.workflow.ParameterManager import ParameterManager params = page_setup() st.title("Filtered PSMs") @@ -75,6 +78,35 @@ st.subheader("MS2 Spectrum") line_plot(key=f"lineplot_{cache_id_prefix}", state_manager=state_manager, height=450, sequence_view_key=f"seqview_{cache_id_prefix}") +st.subheader("Fragment ions") +fragments_tsv = workflow_dir / "results" / "exports" / f"{cache_id_prefix}_fragment_ions.tsv" +psm_id = state_manager.get_selection("identification") +if not fragments_tsv.is_file(): + st.info("No fragment ion export for this run. Re-run the workflow to create it.") +elif psm_id is None: + st.info("Select a PSM in the table to see its fragment mass errors and ion table.") +else: + fragments = load_psm_fragments(fragments_tsv, int(psm_id)) + if fragments.empty: + st.info("No fragment ion of this PSM matched a peak within the tolerance.") + else: + # The search's fragment tolerance: user value, else the workflow's default for it + workflow_params = ParameterManager(workflow_dir, "TOPP Workflow").get_parameters_from_json() + comet = "CometAdapter-TMT" if workflow_params.get("analysis-mode", "LFQ") == "TMT" else "CometAdapter" + comet_params = {**workflow_params.get("_defaults", {}).get(comet, {}), **workflow_params.get(comet, {})} + tolerance = float(comet_params.get("fragment_mass_tolerance", 0.02)) + tolerance_ppm = comet_params.get("fragment_error_units", "Da") != "Da" + + first = fragments.iloc[0] + st.caption( + f"{first['sequence']}, charge {int(first['precursor_charge'])}. Matched with the search's fragment " + f"tolerance (±{tolerance:g} {'ppm' if tolerance_ppm else 'Da'}, shaded)." + ) + show_fig(fragment_error_plot(fragments, tolerance, tolerance_ppm), f"fragment-errors-{cache_id_prefix}-{psm_id}") + values, matched = ion_table(first["sequence"], int(first["precursor_charge"]), fragments) + st.markdown("**b/y ions** (theoretical m/z; matched ions highlighted, neutral losses not counted)") + st.html(ion_table_html(values, matched)) + st.markdown("---") st.markdown("**Next step:** View abundance quantification") st.page_link("content/results_abundance.py", label="Go to Abundance", icon="📋") diff --git a/src/common/fragment_plots.py b/src/common/fragment_plots.py new file mode 100644 index 0000000..1859f36 --- /dev/null +++ b/src/common/fragment_plots.py @@ -0,0 +1,139 @@ +"""Fragment-level evidence for one PSM: mass error plot and b/y ion table. + +Both read the run's ``_fragment_ions.tsv`` export (written by the +workflow with the same ion types and tolerance as the SequenceView), filtered +to the selected PSM. +""" +import html +from pathlib import Path + +import numpy as np +import pandas as pd +import plotly.graph_objects as go +import polars as pl +from pyopenms import AASequence, Residue + +B_COLOR = "#2a78d6" +Y_COLOR = "#eb6834" + + +def load_psm_fragments(tsv_path: Path, sequence_id: int) -> pd.DataFrame: + """The fragment ions matched in one PSM's spectrum, read lazily from the run's export.""" + return ( + pl.scan_csv(tsv_path, separator="\t") + .filter(pl.col("sequence_id") == sequence_id) + .collect() + .to_pandas() + ) + + +def _is_plain(ion: str, ion_type: str, number: int) -> bool: + """Whether an ion label is the plain fragment, not a neutral-loss one (``b4`` vs ``b4-H2O``).""" + return ion == f"{ion_type}{number}" + + +def fragment_error_plot(fragments: pd.DataFrame, tolerance: float, tolerance_ppm: bool) -> go.Figure: + """Mass error of every matched fragment against its m/z, in the search tolerance's unit. + + Dots are sized by peak intensity; the shaded band is the fragment + tolerance. Errors spread evenly around 0 within the band are what a + correct match looks like; a trend with m/z points to a calibration issue. + """ + unit = "ppm" if tolerance_ppm else "Da" + error_col = "mass_error_ppm" if tolerance_ppm else "mass_error_da" + fig = go.Figure() + fig.add_hrect(y0=-tolerance, y1=tolerance, fillcolor="#a1a1aa", opacity=0.15, line_width=0) + max_intensity = fragments["intensity"].max() if len(fragments) else 1.0 + for ion_type, color in (("b", B_COLOR), ("y", Y_COLOR)): + part = fragments[fragments["ion_type"] == ion_type] + if part.empty: + continue + labels = [f"{ion}{'+' * int(z)}" for ion, z in zip(part["ion"], part["charge"])] + fig.add_trace(go.Scatter( + x=part["observed_mz"], y=part[error_col], mode="markers", name=f"{ion_type} ions", + marker=dict(color=color, size=6 + 14 * np.sqrt(part["intensity"] / max_intensity), opacity=0.8), + text=labels, + hovertemplate="%{text}
m/z %{x:.4f}
error %{y:.3f} " + unit + "", + )) + fig.add_hline(y=0, line_color="grey", line_width=1) + fig.update_layout( + title=f"Fragment mass error ({len(fragments)} matched ions)", + xaxis_title="m/z", + yaxis_title=f"Mass error ({unit})", + yaxis_range=[-tolerance * 1.2, tolerance * 1.2], + legend=dict(orientation="h", yanchor="bottom", y=1.02, x=0), + ) + return fig + + +def ion_table(sequence: str, precursor_charge: int, fragments: pd.DataFrame) -> tuple[pd.DataFrame, pd.DataFrame]: + """Theoretical b and y m/z per position, and which of them a peak matched. + + One row per residue; b ions count from the N-terminus and y ions from the + C-terminus, at charges 1 up to ``precursor_charge - 1`` (at most 2). Only + plain ions count as matches, not their neutral losses. + + Returns: + (values, matched): the m/z table and a boolean table of the same shape. + """ + peptide = AASequence.fromString(sequence) + n = peptide.size() + charges = list(range(1, max(1, min(int(precursor_charge) - 1, 2)) + 1)) + plain = fragments[[ + _is_plain(i, t, k) for i, t, k in zip(fragments["ion"], fragments["ion_type"], fragments["ion_number"]) + ]] if len(fragments) else fragments + matched_keys = set(zip(plain["ion_type"], plain["ion_number"], plain["charge"])) if len(plain) else set() + + def column(ion_type, z): + return f"{ion_type}{'+' * z}" + + values = {"#": list(range(1, n + 1)), "residue": [peptide.getResidue(i).getOneLetterCode() for i in range(n)]} + matched = {"#": [False] * n, "residue": [False] * n} + for z in charges: + b_vals, b_hit, y_vals, y_hit = [], [], [], [] + for i in range(1, n + 1): + # b_i covers residues 1..i (no b_n), y_i the last i residues; row i shows b_i and y_(n-i+1) + if i < n: + b_vals.append(peptide.getPrefix(i).getMonoWeight(Residue.ResidueType.BIon, z) / z) + b_hit.append(("b", i, z) in matched_keys) + else: + b_vals.append(np.nan) + b_hit.append(False) + k = n - i + 1 + if k < n: + y_vals.append(peptide.getSuffix(k).getMonoWeight(Residue.ResidueType.YIon, z) / z) + y_hit.append(("y", k, z) in matched_keys) + else: + y_vals.append(np.nan) + y_hit.append(False) + values[column("b", z)], matched[column("b", z)] = b_vals, b_hit + values[column("y", z)], matched[column("y", z)] = y_vals, y_hit + values["#y"] = [n - i + 1 for i in range(1, n + 1)] + matched["#y"] = [False] * n + order = ["#"] + [column("b", z) for z in charges] + ["residue"] + [column("y", z) for z in charges] + ["#y"] + return pd.DataFrame(values)[order], pd.DataFrame(matched)[order] + + +def ion_table_html(values: pd.DataFrame, matched: pd.DataFrame) -> str: + """The ion table as HTML, matched b cells in blue and y cells in orange.""" + head = "".join(f"{html.escape('#' if c == '#y' else c)}" for c in values.columns) + rows = [] + for r in range(len(values)): + cells = [] + for c in values.columns: + v = values.iat[r, values.columns.get_loc(c)] + text = "" if isinstance(v, float) and np.isnan(v) else (f"{v:.4f}" if isinstance(v, float) else html.escape(str(v))) + cls = "" + if matched.iat[r, matched.columns.get_loc(c)]: + cls = ' class="hit-b"' if c.startswith("b") else ' class="hit-y"' + elif c == "residue": + cls = ' class="res"' + cells.append(f"{text}") + rows.append("" + "".join(cells) + "") + return ( + "" + f'{head}{"".join(rows)}
' + ) diff --git a/tests/test_fragment_plots.py b/tests/test_fragment_plots.py new file mode 100644 index 0000000..77498ce --- /dev/null +++ b/tests/test_fragment_plots.py @@ -0,0 +1,63 @@ +import pandas as pd +import pytest +from pyopenms import AASequence, Residue + +from src.common.fragment_plots import fragment_error_plot, ion_table, ion_table_html, load_psm_fragments + + +def _fragments(): + return pd.DataFrame({ + "sequence_id": [0, 0, 0, 0], + "sequence": ["PEPTIDEK"] * 4, + "precursor_charge": [2] * 4, + "ion": ["b2", "b2-H2O", "y3", "y3"], + "ion_type": ["b", "b", "y", "y"], + "ion_number": [2, 2, 3, 3], + "charge": [1, 1, 1, 1], + "theoretical_mz": [227.1, 209.1, 391.2, 391.2], + "observed_mz": [227.1, 209.1, 391.2, 391.21], + "intensity": [100.0, 50.0, 200.0, 20.0], + "mass_error_da": [0.001, -0.002, 0.0005, 0.01], + "mass_error_ppm": [4.0, -9.0, 1.2, 25.0], + "peak_id": [1, 2, 3, 4], + }) + + +def test_ion_table_values_and_matches(): + values, matched = ion_table("PEPTIDEK", 2, _fragments()) + assert list(values.columns) == ["#", "b+", "residue", "y+", "#y"] + peptide = AASequence.fromString("PEPTIDEK") + assert values.loc[1, "b+"] == pytest.approx(peptide.getPrefix(2).getMonoWeight(Residue.ResidueType.BIon, 1)) + # row i holds y_(n-i+1): y3 sits in row 6 (index 5) + assert values.loc[5, "#y"] == 3 + assert values.loc[5, "y+"] == pytest.approx(peptide.getSuffix(3).getMonoWeight(Residue.ResidueType.YIon, 1)) + assert matched.loc[1, "b+"] and matched.loc[5, "y+"] + assert matched["b+"].sum() == 1 # b2-H2O is a neutral loss, not counted + assert pd.isna(values.loc[7, "b+"]) and pd.isna(values.loc[0, "y+"]) + + +def test_ion_table_doubly_charged_columns(): + values, _ = ion_table("PEPTIDEK", 3, _fragments()) + assert list(values.columns) == ["#", "b+", "b++", "residue", "y+", "y++", "#y"] + + +def test_ion_table_html_highlights_matches(): + values, matched = ion_table("PEPTIDEK", 2, _fragments()) + out = ion_table_html(values, matched) + assert out.count('class="hit-b"') == 1 and out.count('class="hit-y"') == 1 + + +def test_fragment_error_plot_uses_tolerance_unit(): + fig = fragment_error_plot(_fragments(), 20, True) + assert fig.layout.yaxis.title.text == "Mass error (ppm)" + assert [t.name for t in fig.data] == ["b ions", "y ions"] + fig_da = fragment_error_plot(_fragments(), 0.02, False) + assert list(fig_da.data[0].y) == [0.001, -0.002] + + +def test_load_psm_fragments_filters_one_psm(tmp_path): + tsv = tmp_path / "run_fragment_ions.tsv" + frags = _fragments() + pd.concat([frags, frags.assign(sequence_id=1)]).to_csv(tsv, sep="\t", index=False) + assert len(load_psm_fragments(tsv, 1)) == 4 + assert load_psm_fragments(tsv, 7).empty