Merge pull request #11 from OpenMS/f5/library-coverage-identifiable-denominator

Fix OpenSwath library coverage denominator

Author: nej1gotnochill

benchmark/openswath_metrics.py

@@ -22,7 +22,11 @@
               `_precursor_mz` metric names.
 - feature     one picked chromatographic feature (a top-level <feature>). The
               runner owns the feature count; this module does not duplicate it.
-- peptide     a distinct TraML <Peptide>.
+- peptide     a distinct TraML <Peptide>.  A declared peptide with at least
+              one transition is "identifiable"; one with none is an "orphan".
+              library_peptides counts BOTH (it is the declared population);
+              library_peptides_identifiable and library_peptides_orphan split
+              it, and library_coverage is measured over the identifiable half.
 - protein     a distinct TraML <Protein> reached through a peptide's
               ProteinRef.
 
@@ -63,6 +67,14 @@
 # scalars (overall_quality_sum rounds to 6 decimals).
 COVERAGE_DECIMALS = 6
 
+# library_peptides_identifiable is the base of library_coverage. A TraML
+# peptide with no transitions is still declared in the library (it counts
+# towards library_peptides, matching upstream TargetedExperiment::getSummary())
+# but it has no precursor, so no transition, so no native_id a feature could
+# carry: it can never be identified. Dividing by such an entry would cap
+# coverage at a ceiling fixed by library construction rather than by the tool,
+# so it is excluded from the ratio and reported separately instead.
+
 
 class OpenSwathMetricsError(Exception):
     """Raised for a structurally unusable TraML/featureXML, or a feature whose
@@ -163,10 +175,16 @@ def count_identifications(library, featurexml_path):
     """Join a featureXML against a parsed library and return the counts.
 
     Returns (all integers unless noted):
-        library_peptides, library_precursor_mz, library_transitions,
+        library_peptides, library_peptides_identifiable,
+        library_peptides_orphan, library_precursor_mz, library_transitions,
         identified_peptides, identified_precursor_mz, identified_proteins,
         transitions_used, library_coverage (float, identified_peptides /
-        library_peptides rounded to COVERAGE_DECIMALS)
+        library_peptides_identifiable rounded to COVERAGE_DECIMALS)
+
+    library_peptides stays the full declared count so an orphan library entry
+    stays visible; library_coverage is divided by the identifiable subset
+    only, because an orphan is structurally incapable of being identified and
+    would otherwise impose a ceiling no tool result can exceed.
 
     Raises OpenSwathMetricsError if the featureXML cannot be parsed or a
     feature references a transition the library does not define.
@@ -208,10 +226,18 @@ def count_identifications(library, featurexml_path):
             identified_proteins.update(peptide["proteins"])
 
     library_peptides = library["library_peptides"]
-    coverage = (len(identified_peptides) / library_peptides
-                if library_peptides else 0.0)
+    # Identifiable = declared peptides that at least one transition references.
+    # Intersecting with the declared set keeps a transition pointing at an
+    # undeclared peptide out of the denominator; the join above still fails
+    # loudly if a feature reaches such a peptide.
+    identifiable = {t["peptide"] for t in transitions.values()} & set(peptides)
+    identifiable_count = len(identifiable)
+    coverage = (len(identified_peptides) / identifiable_count
+                if identifiable_count else 0.0)
     return {
         "library_peptides": library_peptides,
+        "library_peptides_identifiable": identifiable_count,
+        "library_peptides_orphan": library_peptides - identifiable_count,
         "library_precursor_mz": library["library_precursor_mz"],
         "library_transitions": library["library_transitions"],
         "identified_peptides": len(identified_peptides),

benchmark/report/report_generate.py

@@ -201,6 +201,13 @@ def _fmt_metric(name, value):
             return f"{value:.2f} s"
         if name.endswith("peak_rss_kb"):
             return fmt_bytes(float(value) * 1024)
+        # library_coverage is a fraction of a population, not a bare scalar:
+        # render it as a percentage so it cannot be misread as an absolute
+        # count or an error magnitude. The numerator/denominator are shown as
+        # their own rows (identified_peptides / library_peptides_identifiable),
+        # so the fraction itself stays reconstructible from the table.
+        if name.endswith("library_coverage"):
+            return f"{value * 100:.1f}%"
         if name.endswith("_s"):
             return fmt_seconds(value)
         if abs(value) < 10:
@@ -1132,11 +1139,27 @@ def normalize_smoke(args):
 #
 # PSM / FDR counts are absent by design: DIA emits no spectrum-level
 # identifications. See benchmark/openswath_metrics.py for the definitions.
+#
+# library_coverage is measured over library_peptides_identifiable, NOT over
+# library_peptides. A declared peptide with no transitions has no precursor,
+# so no transition, so no native_id a feature could carry: it can never be
+# identified, and counting it in the ratio capped coverage at a ceiling fixed
+# by library construction rather than by the tool. library_peptides stays the
+# full declared count so such an orphan stays visible, and
+# library_peptides_orphan reports it outright.
+#
+# NOTE FOR TRENDS: this changed library_coverage's MEANING, not just its
+# accuracy. Historical values are identified/library_peptides and are NOT
+# comparable with values recorded after this change, which are
+# identified/library_peptides_identifiable. Runs recorded either side of the
+# change will plot as one series with a step at the change point.
 _OPENSWATH_CORRECTNESS_METRICS = {
     "actual_features": "identified_features",
     "overall_quality_sum": "feature_quality_sum",
     "total_intensity": "total_intensity",
     "library_peptides": "library_peptides",
+    "library_peptides_identifiable": "library_peptides_identifiable",
+    "library_peptides_orphan": "library_peptides_orphan",
     "library_precursor_mz": "library_precursor_mz",
     "library_transitions": "library_transitions",
     "identified_peptides": "identified_peptides",

benchmark/report/test_generic_renderer.py

@@ -20,6 +20,7 @@
 import json
 import math
 import os
+import re
 import sys
 import tempfile
 import shutil
@@ -1450,7 +1451,12 @@ def test_render_trends_single_run_and_missing_metric():
 
 def _library_correctness():
     """The correctness block run_openswath_benchmark.sh records after the
-    TraML <-> featureXML join (see benchmark/openswath_metrics.py)."""
+    TraML <-> featureXML join (see benchmark/openswath_metrics.py).
+
+    Mirrors the reference fixture: 7 declared peptides, one of which
+    (PEPTIDEA_Extra) carries no transitions, so 6 are identifiable and all 6
+    were identified -> library_coverage 1.0 over the identifiable base.
+    """
     return {
         "expected_features": 6,
         "actual_features": 6,
@@ -1459,13 +1465,15 @@ def _library_correctness():
         "total_intensity": 149891.59,
         "qc_charge_distribution": [[1, 19], [2, 19], [3, 19]],
         "library_peptides": 7,
+        "library_peptides_identifiable": 6,
+        "library_peptides_orphan": 1,
         "library_precursor_mz": 6,
         "library_transitions": 18,
         "identified_peptides": 6,
         "identified_precursor_mz": 6,
         "identified_proteins": 1,
         "transitions_used": 18,
-        "library_coverage": 0.857143,
+        "library_coverage": 1.0,
     }
 
 
@@ -1495,13 +1503,17 @@ def test_normalize_openswath_lifts_library_metrics():
         # same names, straight from the library join
         assert m["total_intensity"] == 149891.59, m
         assert m["library_peptides"] == 7, m
+        assert m["library_peptides_identifiable"] == 6, m
+        assert m["library_peptides_orphan"] == 1, m
         assert m["library_precursor_mz"] == 6, m
         assert m["library_transitions"] == 18, m
         assert m["identified_peptides"] == 6, m
         assert m["identified_precursor_mz"] == 6, m
         assert m["identified_proteins"] == 1, m
         assert m["transitions_used"] == 18, m
-        assert m["library_coverage"] == 0.857143, m
+        # coverage is over the identifiable base, so a run that finds every
+        # searchable peptide reads 1.0 while the orphan stays visible above
+        assert m["library_coverage"] == 1.0, m
 
         # not measurements -> not lifted
         assert "expected_features" not in m, m
@@ -1546,6 +1558,8 @@ def test_openswath_library_metrics_reach_trends():
         cols = {c["name"]: c["kind"] for c in series[0]["columns"]}
         for key in ("identified_peptides", "identified_precursor_mz",
                     "identified_proteins", "library_peptides",
+                    "library_peptides_identifiable",
+                    "library_peptides_orphan",
                     "library_coverage", "identified_features",
                     "feature_quality_sum", "total_intensity"):
             assert key in cols, (key, sorted(cols))
@@ -1554,8 +1568,10 @@ def test_openswath_library_metrics_reach_trends():
         # the run that predates the metrics keeps a None (rendered as a dash)
         assert series[0]["values"]["identified_peptides"] == [None, 6], \
             series[0]["values"]["identified_peptides"]
-        assert series[0]["values"]["library_coverage"] == [None, 0.857143], \
+        assert series[0]["values"]["library_coverage"] == [None, 1.0], \
             series[0]["values"]["library_coverage"]
+        assert series[0]["values"]["library_peptides_orphan"] == [None, 1], \
+            series[0]["values"]["library_peptides_orphan"]
 
         # and the values reach the rendered HTML through the generic path
         import datetime as _dt
@@ -1566,11 +1582,36 @@ def test_openswath_library_metrics_reach_trends():
             html = fh.read()
         assert "identified_peptides" in html, "metric column not rendered"
         assert "library_coverage" in html, "metric column not rendered"
-        # the cell text is whatever the shared formatter produces for the
-        # metric, not the raw JSON number (_fmt_metric applies its own
-        # precision rules), so assert against the formatter itself.
-        assert _fmt_metric("library_coverage", 0.857143) in html, \
-            "formatted metric value not rendered"
+        assert "library_peptides_identifiable" in html, "metric column not rendered"
+        assert "library_peptides_orphan" in html, "metric column not rendered"
+        # Assert against a LITERAL expected rendering, not a live
+        # _fmt_metric() call. Comparing against the formatter's own output
+        # passes for the wrong reason: any formatter change, including a total
+        # regression, would satisfy it. Instead the coverage CELL is extracted
+        # from the rendered row and compared to the expected literal, so the
+        # formatter is genuinely protected. library_coverage is a fraction of
+        # a population, so it must render as a percentage (1.0 -> "100.0%")
+        # and never as a bare "1.000" scalar.
+        cell = re.search(
+            r"<td class='mono'>library_coverage</td><td class='num'>([^<]*)</td>",
+            html)
+        assert cell, "library_coverage row not found in rendered HTML"
+        assert cell.group(1) == "100.0%", (
+            "library_coverage rendered as %r, expected '100.0%%'" % cell.group(1))
+        # same expectation for the comparison table, whose key is namespaced.
+        # That row is [key | baseline | current | delta]; the baseline run
+        # predates the library metrics, so its cell is an em dash and the
+        # current run's cell is the formatted value.
+        cell_cmp = re.search(
+            r"<td class='mono'>correctness\.library_coverage</td>"
+            r"<td class='num'>([^<]*)</td><td class='num'>([^<]*)</td>", html)
+        assert cell_cmp, "correctness.library_coverage row not found"
+        assert cell_cmp.group(1) == "—", (
+            "expected an em dash for the pre-metrics baseline, got %r"
+            % cell_cmp.group(1))
+        assert cell_cmp.group(2) == "100.0%", (
+            "correctness.library_coverage rendered as %r, expected '100.0%%'"
+            % cell_cmp.group(2))
         print("  PASS: OpenSwath library metrics reach the generic Trends table")
     finally:
         shutil.rmtree(tmpdir)

benchmark/test_openswath_metrics.py

@@ -189,11 +189,17 @@ def test_feature_join_counts():
         assert counts["identified_precursor_mz"] == 2, counts    # one m/z each
         assert counts["identified_proteins"] == 2, counts       # ProtA, ProtB
         assert counts["transitions_used"] == 3, counts
-        assert counts["library_peptides"] == 4, counts          # incl. orphan
+        assert counts["library_peptides"] == 4, counts          # incl. orphans
+        # Only PepA and PepB carry transitions; PepC and PepOrphan are both
+        # declared-but-transition-less, so 4 declared splits 2 + 2.
+        assert counts["library_peptides_identifiable"] == 2, counts  # PepA, PepB
+        assert counts["library_peptides_orphan"] == 2, counts   # PepC, PepOrphan
         assert counts["library_precursor_mz"] == 2, counts
         assert counts["library_transitions"] == 3, counts
-        # 2 of 4 library peptides
-        assert counts["library_coverage"] == 0.5, counts
+        # both IDENTIFIABLE peptides were found, so coverage is 1.0 over the
+        # identifiable base even though the library declares twice as many
+        # peptides as can actually be searched.
+        assert counts["library_coverage"] == 1.0, counts
     print("PASS: feature join (distinct peptides/precursors/proteins, coverage)")
 
 
@@ -235,6 +241,10 @@ def test_empty_and_missing_inputs():
         assert counts["identified_peptides"] == 0, counts
         assert counts["library_coverage"] == 0.0, counts
         assert counts["library_peptides"] == 2, counts
+        # PepB is declared but transition-less, so only 1 peptide is
+        # identifiable even though the library declares 2.
+        assert counts["library_peptides_identifiable"] == 1, counts
+        assert counts["library_peptides_orphan"] == 1, counts
 
         for bad in (os.path.join(tmp, "nope.featureXML"),
                     os.path.join(tmp, "broken.TraML")):
@@ -254,8 +264,11 @@ def test_empty_and_missing_inputs():
 
 def test_reference_fixture_join():
     """Regression guard on the real committed OpenSwath fixtures: the vendored
-    library has 7 peptides but only 6 precursor m/z values, because one library
-    peptide carries no transitions and therefore can never be identified."""
+    library declares 7 peptides but only 6 carry transitions, because
+    PEPTIDEA_Extra has no transitions and therefore can never be identified.
+    The run is perfect on every other axis (18/18 transitions, 6/6 precursor
+    m/z), so coverage must be 6/6 = 1.0 over the identifiable subset, while
+    the orphan stays visible as its own count instead of penalising it."""
     traml = os.path.join(FIXTURES, "OpenSwathWorkflow_1_input.TraML")
     feat = os.path.join(FIXTURES, "OpenSwathWorkflow_1_output.featureXML")
     lib = om.parse_traml(traml)
@@ -268,12 +281,91 @@ def test_reference_fixture_join():
     assert counts["identified_precursor_mz"] == 6, counts
     assert counts["identified_proteins"] == 1, counts
     assert counts["transitions_used"] == 18, counts
-    assert counts["library_coverage"] == 0.857143, counts
+    # declared / identifiable / orphan must all stay observable
+    assert counts["library_peptides"] == 7, counts
+    assert counts["library_peptides_identifiable"] == 6, counts
+    assert counts["library_peptides_orphan"] == 1, counts
+    assert counts["library_peptides_orphan"] == (
+        counts["library_peptides"] - counts["library_peptides_identifiable"]), counts
+    assert counts["library_coverage"] == 1.0, counts
     # No PSM/FDR key may ever appear: DIA emits no spectrum-level
     # identifications, so such a count would be fabricated.
     assert not [k for k in counts if "psm" in k.lower()], counts
     assert not [k for k in counts if "fdr" in k.lower()], counts
-    print("PASS: reference fixture join (7 peptides / 6 precursor m/z / 6 identified)")
+    print("PASS: reference fixture join (7 declared / 6 identifiable / 6 identified "
+          "-> coverage 1.0)")
+
+
+def test_orphan_peptide_does_not_penalise_coverage():
+    """F5 regression guard: a declared-but-transition-less peptide must not
+    lower library_coverage.
+
+    The fixture library carries one such orphan (PEPTIDEA_Extra). Two runs
+    that find every IDENTIFIABLE peptide must both report coverage 1.0,
+    whether or not the library also declares orphans. Before this guard the
+    metric divided by all declared peptides, so the same perfect run reported
+    6/7 = 0.857143 and looked like a 14% identification failure caused by the
+    tool. If this test ever fails, either the denominator changed again or the
+    orphan counts changed -- both require the module docstring to change too.
+    """
+    traml = os.path.join(FIXTURES, "OpenSwathWorkflow_1_input.TraML")
+    feat = os.path.join(FIXTURES, "OpenSwathWorkflow_1_output.featureXML")
+    lib = om.parse_traml(traml)
+    counts = om.count_identifications(lib, feat)
+
+    # every identifiable peptide was found
+    assert counts["identified_peptides"] == counts["library_peptides_identifiable"], \
+        counts
+    # so coverage is exactly 1.0 despite the orphan being declared
+    assert counts["library_peptides_orphan"] > 0, counts   # the orphan is present
+    assert counts["library_coverage"] == 1.0, counts
+
+    # Same result with an orphan-free library: 3 declared, 3 identifiable,
+    # all 3 identified -> coverage 1.0. Identifiable coverage is therefore
+    # independent of how many orphans the library happens to declare.
+    with tempfile.TemporaryDirectory() as tmp:
+        lib2 = os.path.join(tmp, "noorphan.TraML")
+        make_traml(
+            lib2,
+            proteins=[("ProtA", "uniprot_a")],
+            peptides=[("PepA", "AAAA", "2", "ProtA"),
+                      ("PepB", "BBBB", "2", "ProtA"),
+                      ("PepC", "CCCC", "2", "ProtA")],
+            transitions=[("t1", "PepA", "400.5"),
+                         ("t2", "PepB", "500.5"),
+                         ("t3", "PepC", "600.5")],
+        )
+        feat2 = os.path.join(tmp, "all.featureXML")
+        make_featurexml(feat2, [["t1"], ["t2"], ["t3"]])
+        c2 = om.count_identifications(om.parse_traml(lib2), feat2)
+        assert c2["library_peptides"] == 3, c2
+        assert c2["library_peptides_identifiable"] == 3, c2
+        assert c2["library_peptides_orphan"] == 0, c2
+        assert c2["library_coverage"] == 1.0, c2
+
+    # A genuine miss must still be penalised: drop one peptide's feature and
+    # coverage falls below 1.0. This is what stops the orphan exclusion from
+    # degenerating into "coverage is always 1.0".
+    with tempfile.TemporaryDirectory() as tmp:
+        lib3 = os.path.join(tmp, "miss.TraML")
+        make_traml(
+            lib3,
+            proteins=[("ProtA", "uniprot_a")],
+            peptides=[("PepA", "AAAA", "2", "ProtA"),
+                      ("PepB", "BBBB", "2", "ProtA"),
+                      ("PepC", "CCCC", "2", "ProtA")],
+            transitions=[("t1", "PepA", "400.5"),
+                         ("t2", "PepB", "500.5"),
+                         ("t3", "PepC", "600.5")],
+        )
+        feat3 = os.path.join(tmp, "partial.featureXML")
+        make_featurexml(feat3, [["t1"], ["t2"]])   # PepC not identified
+        c3 = om.count_identifications(om.parse_traml(lib3), feat3)
+        assert c3["identified_peptides"] == 2, c3
+        assert c3["library_peptides_identifiable"] == 3, c3
+        assert c3["library_coverage"] == 0.666667, c3
+
+    print("PASS: orphan peptide does not penalise coverage; real misses still do")
 
 
 def test_isobaric_peptides_share_one_precursor_mz():
@@ -325,6 +417,7 @@ def main():
     test_unmatched_transition_fails_loudly()
     test_empty_and_missing_inputs()
     test_reference_fixture_join()
+    test_orphan_peptide_does_not_penalise_coverage()
     test_isobaric_peptides_share_one_precursor_mz()
     print("\nAll openswath_metrics tests passed.")