Add an Identification QC page (score distribution, precursor error, PSMs per file)

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_qc.py"), title="Identification QC", icon="🩺"),
             st.Page(Path("content", "results_abundance.py"), title="Protein Abundance", icon="📋"),
             st.Page(Path("content", "results_download.py"), title="Download", icon="⬇️"),
         ],

content/results_qc.py

@@ -0,0 +1,84 @@
+"""Identification QC Page."""
+from pathlib import Path
+
+import pandas as pd
+import streamlit as st
+
+from src.common.common import page_setup, show_fig
+from src.common.qc_plots import psm_qc_table, precursor_error_histogram, psms_per_file, score_distribution
+from src.common.results_helpers import get_workflow_dir
+from src.workflow.ParameterManager import ParameterManager
+
+params = page_setup()
+
+st.title("Identification QC")
+st.markdown(
+    """
+Quality checks of the identifications across all runs: how well targets separate
+from decoys, how well the precursor masses are calibrated, and whether one run
+identified clearly fewer spectra than the others.
+"""
+)
+
+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"
+workflow_params = ParameterManager(workflow_dir, "TOPP Workflow").get_parameters_from_json()
+fdr = min(float(workflow_params.get("psm-fdr-percent", 1.0)), 100.0) / 100.0
+mzml_names = [Path(f).name for f in workflow_params.get("mzML-files", []) or []]
+
+
+@st.cache_data(show_spinner="Reading identifications...")
+def load_qc_table(paths: tuple[str, ...], mtimes: tuple[float, ...], mzml_names: tuple[str, ...]) -> pd.DataFrame:
+    # mtimes only key the cache, so a re-run of the workflow invalidates it
+    return pd.concat([psm_qc_table(Path(p), list(mzml_names)) for p in paths], ignore_index=True)
+
+
+def qc_table(directory: Path) -> pd.DataFrame | None:
+    files = sorted(directory.glob("*.idXML"))
+    if not files:
+        return None
+    return load_qc_table(
+        tuple(str(f) for f in files), tuple(f.stat().st_mtime for f in files), tuple(mzml_names)
+    )
+
+
+# All PSMs before filtering, so decoys and the FDR cutoff are visible
+all_psms = qc_table(results_dir / "percolator_results")
+if all_psms is None:
+    all_psms = qc_table(results_dir / "comet_results")
+filtered = qc_table(results_dir / "psm_filter")
+
+if all_psms is None and filtered is None:
+    st.info("No identification results available yet. Please run the workflow first.")
+    st.page_link("content/workflow_run.py", label="Go to Run", icon="🚀")
+    st.stop()
+
+if all_psms is not None and not all_psms.empty:
+    st.subheader("Score distribution")
+    st.caption(
+        "All PSMs before FDR filtering. Decoys model the incorrect matches: targets to the "
+        "right of the decoy peak are the correct identifications. The dashed line is the "
+        "lowest target score that passes the PSM FDR."
+    )
+    if not all_psms["decoy"].any():
+        st.warning("No decoy PSMs found, so the score distribution cannot show the incorrect matches.")
+    show_fig(score_distribution(all_psms, fdr), "qc-score-distribution")
+
+if filtered is not None and not filtered.empty:
+    st.subheader("Precursor mass error")
+    st.caption(
+        f"Target PSMs at the PSM FDR ({fdr * 100:.3g} %). A median far from 0 ppm points to a "
+        "calibration offset; a distribution as wide as the precursor tolerance suggests the "
+        "tolerance could be tighter."
+    )
+    show_fig(precursor_error_histogram(filtered), "qc-precursor-error")
+
+    st.subheader("PSMs per file")
+    st.caption(f"Target PSMs at the PSM FDR ({fdr * 100:.3g} %).")
+    show_fig(psms_per_file(filtered), "qc-psms-per-file")
+elif all_psms is not None:
+    st.info("No filtered PSMs yet: the precursor error and PSMs per file appear once filtering has run.")

src/common/qc_plots.py

@@ -0,0 +1,157 @@
+"""Identification QC: score distribution, precursor mass error and PSMs per file.
+
+``psm_qc_table`` reads one idXML into a flat pandas table (one row per PSM, top
+hit only); the plot functions take the concatenated tables of all runs and
+return Plotly figures, so the page only decides where to show them.
+"""
+from pathlib import Path
+
+import numpy as np
+import pandas as pd
+import plotly.graph_objects as go
+
+from src.common.results_helpers import extract_filename_from_idxml, load_idxml
+
+# Spacing of the 13C isotope peaks (Da): the precursor may have been picked on one of them
+C13_C12_DIFF = 1.0033548
+MAX_ISOTOPE_ERROR = 3
+
+# Percolator writes its scores as PSI-MS meta values on each hit
+PERCOLATOR_SCORE_KEYS = ("MS:1001492", "Percolator_score")
+Q_VALUE_KEYS = ("MS:1001491", "q-value", "Percolator_qvalue")
+
+TARGET_COLOR = "#2a78d6"
+DECOY_COLOR = "#eb6834"
+
+QC_COLUMNS = ["file", "decoy", "score", "score_type", "percolator_score", "q_value", "charge", "precursor_error_ppm"]
+
+
+def _meta(obj, keys):
+    for key in keys:
+        if obj.metaValueExists(key):
+            try:
+                return float(obj.getMetaValue(key))
+            except (TypeError, ValueError):
+                return np.nan
+    return np.nan
+
+
+def _is_decoy(hit, decoy_prefixes=("DECOY_", "rev_")) -> bool:
+    """A hit is a decoy if PeptideIndexer said so, else if every protein carries a decoy prefix."""
+    if hit.metaValueExists("target_decoy"):
+        value = hit.getMetaValue("target_decoy")
+        if isinstance(value, bytes):
+            value = value.decode()
+        return str(value) == "decoy"
+    accessions = [ev.getProteinAccession() for ev in hit.getPeptideEvidences()]
+    return bool(accessions) and all(a.startswith(decoy_prefixes) for a in accessions)
+
+
+def precursor_error_ppm(observed_mz: float, theoretical_mz: float, charge: int) -> float:
+    """Precursor error in ppm, after removing a 13C isotope error of up to three peaks."""
+    if not charge or theoretical_mz <= 0:
+        return np.nan
+    isotope = round((observed_mz - theoretical_mz) * charge / C13_C12_DIFF)
+    if abs(isotope) <= MAX_ISOTOPE_ERROR:
+        observed_mz -= isotope * C13_C12_DIFF / charge
+    return (observed_mz - theoretical_mz) / theoretical_mz * 1e6
+
+
+def psm_qc_table(idxml_path: Path, mzml_names=None) -> pd.DataFrame:
+    """One row per spectrum (its top hit) of an idXML file."""
+    _, peptides = load_idxml(idxml_path)
+    file = extract_filename_from_idxml(Path(idxml_path), mzml_names)
+    records = []
+    for pep in peptides:
+        hits = pep.getHits()
+        if not hits:
+            continue
+        hit = hits[0]
+        charge = hit.getCharge()
+        theoretical = hit.getSequence().getMZ(charge) if charge else np.nan
+        records.append({
+            "file": file,
+            "decoy": _is_decoy(hit),
+            "score": hit.getScore(),
+            "score_type": pep.getScoreType(),
+            "percolator_score": _meta(hit, PERCOLATOR_SCORE_KEYS),
+            "q_value": _meta(hit, Q_VALUE_KEYS),
+            "charge": charge,
+            "precursor_error_ppm": precursor_error_ppm(pep.getMZ(), theoretical, charge),
+        })
+    return pd.DataFrame(records, columns=QC_COLUMNS)
+
+
+def _score_axis(df: pd.DataFrame) -> tuple[pd.Series, str]:
+    """The score to plot, higher is better: Percolator's SVM score, else the search engine's."""
+    if df["percolator_score"].notna().all():
+        return df["percolator_score"], "Percolator score"
+    score_type = str(df["score_type"].iloc[0]) if len(df) else "score"
+    if "e-value" in score_type.lower() or "expect" in score_type.lower():
+        return -np.log10(df["score"].clip(lower=1e-300)), "-log10(E-value)"
+    return df["score"], score_type
+
+
+def score_distribution(df: pd.DataFrame, fdr: float | None = None) -> go.Figure:
+    """Overlaid histograms of target and decoy PSM scores.
+
+    Decoys should sit at low scores and targets show a second mode of correct
+    matches at high scores; the dashed line marks the lowest target score that
+    passes the PSM FDR.
+    """
+    scores, label = _score_axis(df)
+    fig = go.Figure()
+    for is_decoy, name, color in ((False, "Targets", TARGET_COLOR), (True, "Decoys", DECOY_COLOR)):
+        values = scores[df["decoy"] == is_decoy].dropna()
+        fig.add_trace(go.Histogram(x=values, name=f"{name} ({len(values):,})", marker_color=color, opacity=0.6, nbinsx=60))
+    if fdr is not None and fdr < 1.0:
+        passing = scores[(~df["decoy"]) & (df["q_value"] <= fdr)].dropna()
+        if len(passing):
+            fig.add_vline(
+                x=passing.min(), line_dash="dash", line_color="grey",
+                annotation_text=f"PSM FDR {fdr * 100:.3g} % (q ≤ {fdr:.3g})",
+            )
+    fig.update_layout(
+        barmode="overlay",
+        title="Score distribution of target and decoy PSMs",
+        xaxis_title=label,
+        yaxis_title="PSMs",
+        legend=dict(orientation="h", yanchor="bottom", y=1.02, x=0),
+    )
+    return fig
+
+
+def precursor_error_histogram(df: pd.DataFrame) -> go.Figure:
+    """Precursor mass error of target PSMs; the median shows the calibration offset."""
+    errors = df.loc[~df["decoy"], "precursor_error_ppm"].dropna()
+    fig = go.Figure(go.Histogram(x=errors, marker_color=TARGET_COLOR, nbinsx=60, name="Target PSMs"))
+    if len(errors):
+        median = float(errors.median())
+        fig.add_vline(x=median, line_dash="dash", line_color="grey", annotation_text=f"median {median:+.2f} ppm")
+        radius = max(1.0, float(errors.abs().quantile(0.99)))
+        fig.update_xaxes(range=[-radius * 1.05, radius * 1.05])
+    fig.update_layout(
+        title=f"Precursor mass error of target PSMs ({len(errors):,})",
+        xaxis_title="Precursor error (ppm), isotope errors removed",
+        yaxis_title="PSMs",
+        showlegend=False,
+    )
+    return fig
+
+
+def psms_per_file(df: pd.DataFrame) -> go.Figure:
+    """Target PSMs per run, so a failed or weak run stands out."""
+    counts = df.loc[~df["decoy"]].groupby("file").size()
+    files = sorted(df["file"].unique())
+    counts = counts.reindex(files, fill_value=0)
+    fig = go.Figure(go.Bar(
+        x=counts.values, y=counts.index, orientation="h", marker_color=TARGET_COLOR,
+        text=[f"{v:,}" for v in counts.values], textposition="auto",
+    ))
+    fig.update_layout(
+        title="Target PSMs per file",
+        xaxis_title="Target PSMs",
+        yaxis=dict(autorange="reversed"),
+        height=max(250, 40 * len(files) + 120),
+    )
+    return fig

tests/test_qc_plots.py

@@ -0,0 +1,106 @@
+import math
+
+import pandas as pd
+import pyopenms as poms
+import pytest
+
+from src.common.qc_plots import (
+    C13_C12_DIFF,
+    precursor_error_histogram,
+    precursor_error_ppm,
+    psm_qc_table,
+    psms_per_file,
+    score_distribution,
+)
+
+try:  # pyOpenMS >= 3.5.0
+    from pyopenms import PeptideIdentificationList
+except ImportError:  # pyOpenMS <= 3.4.x
+    PeptideIdentificationList = None
+
+
+def _psm(sequence, charge, mz_offset_ppm, decoy, svm, q, isotope=0):
+    seq = poms.AASequence.fromString(sequence)
+    theo = seq.getMZ(charge)
+    hit = poms.PeptideHit()
+    hit.setSequence(seq)
+    hit.setCharge(charge)
+    hit.setScore(0.01)
+    hit.setMetaValue("target_decoy", "decoy" if decoy else "target")
+    hit.setMetaValue("MS:1001492", svm)
+    hit.setMetaValue("MS:1001491", q)
+    pep = poms.PeptideIdentification()
+    pep.setScoreType("pep")
+    pep.setMZ(theo * (1 + mz_offset_ppm * 1e-6) + isotope * C13_C12_DIFF / charge)
+    pep.setRT(100.0)
+    pep.setHits([hit])
+    return pep
+
+
+def _write_idxml(path, peptides):
+    if PeptideIdentificationList:
+        pep_list = PeptideIdentificationList()
+        for p in peptides:
+            pep_list.push_back(p)
+    else:
+        pep_list = peptides
+    poms.IdXMLFile().store(str(path), [poms.ProteinIdentification()], pep_list)
+
+
+def test_precursor_error_removes_isotope_error():
+    theo = 500.0
+    assert precursor_error_ppm(theo * (1 + 5e-6), theo, 2) == pytest.approx(5.0, abs=1e-6)
+    shifted = theo * (1 + 5e-6) + C13_C12_DIFF / 2
+    assert precursor_error_ppm(shifted, theo, 2) == pytest.approx(5.0, abs=1e-6)
+    assert math.isnan(precursor_error_ppm(theo, theo, 0))
+
+
+def test_psm_qc_table_reads_hits(tmp_path):
+    path = tmp_path / "run1_comet_perc.idXML"
+    _write_idxml(path, [
+        _psm("PEPTIDEK", 2, 3.0, False, 2.5, 0.001),
+        _psm("PEPTIDER", 2, -4.0, False, 1.0, 0.02, isotope=1),
+        _psm("KEDITPEP", 2, 0.0, True, -1.0, 0.5),
+    ])
+    df = psm_qc_table(path, ["run1.mzML"])
+    assert list(df["file"].unique()) == ["run1.mzML"]
+    assert df["decoy"].tolist() == [False, False, True]
+    assert df["percolator_score"].tolist() == [2.5, 1.0, -1.0]
+    assert df["q_value"].tolist() == pytest.approx([0.001, 0.02, 0.5])
+    assert df["precursor_error_ppm"].tolist()[:2] == pytest.approx([3.0, -4.0], abs=1e-3)
+
+
+def test_figures_build(tmp_path):
+    df = pd.DataFrame({
+        "file": ["a.mzML", "a.mzML", "b.mzML", "b.mzML"],
+        "decoy": [False, True, False, False],
+        "score": [0.01, 0.9, 0.02, 0.03],
+        "score_type": ["pep"] * 4,
+        "percolator_score": [2.0, -1.0, 1.5, 0.5],
+        "q_value": [0.001, 0.5, 0.005, 0.02],
+        "charge": [2, 2, 3, 2],
+        "precursor_error_ppm": [1.0, 5.0, -2.0, 0.5],
+    })
+    fig = score_distribution(df, fdr=0.01)
+    assert fig.layout.xaxis.title.text == "Percolator score"
+    assert [t.name for t in fig.data] == ["Targets (3)", "Decoys (1)"]
+    # the cutoff line sits at the lowest passing target score
+    assert fig.layout.shapes[0].x0 == 1.5
+
+    err = precursor_error_histogram(df)
+    assert len(err.data[0].x) == 3  # decoys left out
+
+    bars = psms_per_file(df)
+    assert list(bars.data[0].y) == ["a.mzML", "b.mzML"]
+    assert list(bars.data[0].x) == [1, 2]
+
+
+def test_score_distribution_falls_back_to_evalue():
+    df = pd.DataFrame({
+        "file": ["a.mzML"] * 2, "decoy": [False, True], "score": [1e-5, 1.0],
+        "score_type": ["expect"] * 2, "percolator_score": [float("nan")] * 2,
+        "q_value": [float("nan")] * 2, "charge": [2, 2], "precursor_error_ppm": [0.0, 0.0],
+    })
+    fig = score_distribution(df)
+    assert fig.layout.xaxis.title.text == "-log10(E-value)"
+    assert list(fig.data[0].x) == pytest.approx([5.0])