diff --git a/app.py b/app.py index 3e3b047..746e5f0 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_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="⬇️"), ], diff --git a/content/results_qc.py b/content/results_qc.py new file mode 100644 index 0000000..ab74a88 --- /dev/null +++ b/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.") diff --git a/src/common/qc_plots.py b/src/common/qc_plots.py new file mode 100644 index 0000000..199aa5c --- /dev/null +++ b/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 diff --git a/tests/test_qc_plots.py b/tests/test_qc_plots.py new file mode 100644 index 0000000..3a1ec81 --- /dev/null +++ b/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])