Astral PSM gap: feature-tolerance fix + iter19/20/21/22b parity cleanups (+4,574 PSMs @ 1% FDR)

benchmark/parity-fixtures/astral_mods_rust.txt

@@ -0,0 +1,22 @@
+# Rust-format mods.txt for the Astral ProteoBench Module 8 benchmark.
+#
+# Mirrors the Java mods.txt at
+#   benchmark/data/Astral_ProteoBench_Module_8/mods.txt
+# but with NUMERIC mass deltas instead of Java's composition strings
+# (C2H3N1O1, etc.) which msgf-rust's --mod parser does not yet support
+# (`see crates/msgf-rust/src/bin/msgf-rust.rs --mod help`).
+#
+# Pass to msgf-rust as `--mod benchmark/parity-fixtures/astral_mods_rust.txt`
+# when running the Astral bench. Failing to pass it yields a 1% FDR
+# deficit of ~1% (iter22b: 31,006 → iter24: 31,390 with the mod) due to
+# Rust's built-in default (Carbamidomethyl-C + Oxidation-M only)
+# missing Acetyl-Prot-N-term, which 8.4% of Java's label-flip wins use.
+#
+# Audit + bench results: docs/parity-analysis/notes/2026-05-21-iter24-acetyl-mod-fix.md
+#
+# Format: <mass_delta_Da>,<residue|*>,<fix|opt>,<location>,<name>
+# Locations: any | N-term | C-term | Prot-N-term | Prot-C-term
+NumMods=3
+57.02146,C,fix,any,Carbamidomethyl
+15.99491,M,opt,any,Oxidation
+42.01057,*,opt,Prot-N-term,Acetyl

benchmark/parity/README.md

@@ -43,6 +43,26 @@ Outputs:
   placeholder (Rust's i32::MIN for GF-uncomputed PSMs) are excluded from
   per-feature stats to avoid skewing the mean.
 
+### Mod-file parity (CRITICAL)
+
+Rust's `msgf-rust` binary, when invoked without `--mod`, uses built-in
+defaults (Carbamidomethyl-C fixed + Oxidation-M variable). Java MS-GF+
+typically reads a `mods.txt` configured for the dataset, which for
+HCD-on-protein-N-term datasets (Astral ProteoBench Module 8 included)
+adds **Acetyl-Prot-N-term** as a third variable mod.
+
+If Rust runs without `--mod` against a Java run that uses
+`Acetyl-Prot-N-term`, **Rust silently misses all acetylated peptides**.
+Iter24 measured ~+384 PSMs at 1% FDR (gap to Java 13.5% → 12.4%) just
+from passing the equivalent mods file.
+
+Use the Rust-compatible fixture at
+`benchmark/parity-fixtures/astral_mods_rust.txt` (numeric mass deltas;
+Rust's `--mod` parser does NOT yet support Java composition strings
+like `C2H3N1O1`) when reproducing the Astral parity harness against Java.
+
+Doc: `docs/parity-analysis/notes/2026-05-21-iter24-acetyl-mod-fix.md`.
+
 ## Dependencies
 
 Python stdlib only. No pandas / numpy / scipy. Matches the project's

benchmark/parity/analyze_rust_java_pin_diff.py

@@ -229,6 +229,10 @@ def main():
     # ── per-scan top-1 buckets ─────────────────────────────────────────────
     buckets = Counter()
     bucket_examples = defaultdict(list)
+    # iter38 (diagnostic upgrade): collect ALL scans in the 3 non-converging
+    # buckets along with both engines' top-1 rows. Written to non_converging.csv
+    # to enable per-feature analysis of the residual PSM-divergence buckets.
+    non_converging: list = []  # (scan, bucket, j_row, r_row) tuples
     for scan in all_scans:
         j = best_by_lnspec(java_by_scan.get(scan, []))
         r = best_by_lnspec(rust_by_scan.get(scan, []))
@@ -262,18 +266,21 @@ def main():
                     bucket_examples["both_target_diff_peptide"].append(
                         (scan, j["_peptide_residues"], r["_peptide_residues"])
                     )
+                non_converging.append((scan, "both_target_diff_peptide", j, r))
         elif j_lab == 1 and r_lab == -1:
             buckets["java_target_rust_decoy"] += 1
             if len(bucket_examples["java_target_rust_decoy"]) < 5:
                 bucket_examples["java_target_rust_decoy"].append(
                     (scan, j["_peptide_residues"], r["_peptide_residues"])
                 )
+            non_converging.append((scan, "java_target_rust_decoy", j, r))
         elif j_lab == -1 and r_lab == 1:
             buckets["rust_target_java_decoy"] += 1
             if len(bucket_examples["rust_target_java_decoy"]) < 5:
                 bucket_examples["rust_target_java_decoy"].append(
                     (scan, j["_peptide_residues"], r["_peptide_residues"])
                 )
+            non_converging.append((scan, "rust_target_java_decoy", j, r))
         else:
             buckets["both_decoy"] += 1
 
@@ -334,6 +341,34 @@ def main():
             ])
     print(f"Wrote {csv_path} ({len(per_row)} rows)", file=sys.stderr)
 
+    # iter38 diagnostic upgrade: dump per-feature values for the
+    # 3 non-converging buckets so future audits can characterize which
+    # features are driving the divergence. Each row has both engines'
+    # top-1 PSM values side-by-side for every NUMERIC_COL feature.
+    nc_path = args.out_dir / "non_converging.csv"
+    with open(nc_path, "w", newline="") as f:
+        w = csv.writer(f)
+        header = ["scan", "bucket", "java_peptide", "rust_peptide"]
+        for c in NUMERIC_COLS:
+            header.append(f"java_{c}")
+            header.append(f"rust_{c}")
+        w.writerow(header)
+        for scan, bucket, j, r in non_converging:
+            row = [
+                scan,
+                bucket,
+                j.get("_peptide_residues", ""),
+                r.get("_peptide_residues", ""),
+            ]
+            for c in NUMERIC_COLS:
+                jv = j.get(c, "")
+                rv = r.get(c, "")
+                row.append(jv)
+                row.append(rv)
+            w.writerow(row)
+    print(f"Wrote {nc_path} ({len(non_converging)} non-converging PSMs)",
+          file=sys.stderr)
+
     # ── write markdown report ──────────────────────────────────────────────
     md_path = args.out_dir / "report.md"
     with open(md_path, "w") as f:

docs/parity-analysis/diff/iter27/report.md

@@ -0,0 +1,106 @@
+# Rust↔Java PIN diff report
+
+- Java:  `/srv/data/msgf-bench/bench-merged-results/astral-java.pin` (136271 rows, 121654 scans)
+- Rust:  `/srv/data/msgf-bench/bench-iter27-results/astral-rust-iter27.pin` (149685 rows, 121677 scans)
+- Total unique scans: 121681
+
+## Top-1-per-scan buckets
+
+| Bucket | Count | % of total |
+|---|---:|---:|
+| both_target_same_peptide | 45,870 | 37.70% |
+| both_target_diff_peptide | 18,892 | 15.53% |
+| java_target_rust_decoy | 16,623 | 13.66% |
+| rust_target_java_decoy | 13,632 | 11.20% |
+| java_only_target | 3 | 0.00% |
+| rust_only_target | 14 | 0.01% |
+| both_decoy | 26,633 | 21.89% |
+| java_only_decoy | 1 | 0.00% |
+| rust_only_decoy | 13 | 0.01% |
+| both_missing | 0 | 0.00% |
+| **total** | **121,681** | 100.00% |
+
+### Sample disagreements
+
+**both_target_diff_peptide** (first 5):
+
+- `(11, 'QKYFYR', 'QSELDRR')`
+- `(13, 'TPEEGEK', 'MTYTTR')`
+- `(14, 'MLQTPESR', 'MSNRNNNK')`
+- `(15, 'PEWTQR', 'SSRSSHR')`
+- `(29, 'MPPAGGPR', 'MAAQGAPR')`
+
+**rust_target_java_decoy** (first 5):
+
+- `(20, 'EQCDTK', 'QMNDEK')`
+- `(23, 'LHNYEDMLEKNK', 'FQARCCPLQNQK')`
+- `(25, 'RYMEYR', 'PYSYNHR')`
+- `(33, 'QMAASPAK', 'SSENPLR')`
+- `(37, 'WQDFTK', 'YCAERK')`
+
+**java_target_rust_decoy** (first 5):
+
+- `(21, 'NEEQSR', 'TEAPCGK')`
+- `(45, 'RCSEDK', 'EKDCDK')`
+- `(47, 'VWFSQIEYIVLR', 'FQLLEKYEPLNR')`
+- `(48, 'VEEQEK', 'TAQEWK')`
+- `(58, 'ELYETR', 'FMELQK')`
+
+**java_only_target** (first 3):
+
+- `(13201, 'MCYGYGCGCGSFCR')`
+- `(14402, 'MDPNCSCATGGSCSCASSCKCK')`
+- `(30441, 'MDLSCSCATGGSCTCASSCK')`
+
+**rust_only_target** (first 5):
+
+- `(106652, 'DCPSCK')`
+- `(113243, 'LYYALK')`
+- `(113285, 'MKRGQR')`
+- `(114265, 'LIVLKEK')`
+- `(114482, 'YNEGCR')`
+
+## Per-feature diff (agreement bucket: same scan + peptide, both target)
+
+_50,466 PSMs in agreement bucket._
+
+Sorted by mean |Δ| (Rust - Java), descending:
+
+| Feature | n | mean Δ | median Δ | stdev | p5 | p95 | mean \|Δ\| | mean rel Δ | %frac \|relΔ\|>1% |
+|---|---:|---:|---:|---:|---:|---:|---:|---:|---:|
+| MS2IonCurrent | 50,466 | +252.7 | +0 | 7909 | -4 | +4.1 | 254.1 | +0.000258 | 0.2% |
+| DeNovoScore | 50,466 | -13.09 | -13 | 12.28 | -32 | +4 | 14.89 | -0.1788 | 98.1% |
+| RawScore | 50,466 | -2.891 | -2 | 10.48 | -21 | +12 | 8.562 | +0.1294 | 96.4% |
+| lnEValue | 50,466 | -6.758 | -6.881 | 4.087 | -13.12 | +0.1029 | 7.073 | -11.17 | 99.6% |
+| lnSpecEValue | 50,466 | -2.189 | -2.325 | 4.115 | -8.584 | +4.792 | 3.67 | -0.1291 | 96.4% |
+| MeanRelErrorTop7 | 50,466 | +1.338 | +1.462 | 2.744 | -3.254 | +5.475 | 2.269 | -6.267 | 99.5% |
+| NumMatchedMainIons | 50,466 | -1.296 | -1 | 2.36 | -5 | +2 | 1.961 | -0.147 | 80.2% |
+| MeanErrorTop7 | 50,466 | -1.695 | -1.466 | 1.89 | -5.094 | +0.7979 | 1.941 | -0.2931 | 99.1% |
+| StdevRelErrorTop7 | 50,466 | -1.009 | -0.418 | 2.676 | -6.548 | +2.398 | 1.722 | -0.08367 | 98.1% |
+| StdevErrorTop7 | 50,466 | -0.451 | -0.3355 | 1.703 | -3.885 | +2.145 | 1.154 | -0.07068 | 98.5% |
+| longest_y | 50,466 | -0.1274 | +0 | 0.848 | -1 | +0 | 0.22 | -0.009756 | 9.5% |
+| matchedIonRatio | 50,466 | +0.04442 | +0.02381 | 0.2151 | -0.2848 | +0.425 | 0.1716 | +0.01366 | 96.9% |
+| lnDeltaSpecEValue | 50,466 | -0.1313 | +0 | 1.255 | +0 | +0 | 0.1313 | +nan | — |
+| longest_y_pct | 50,466 | -0.06723 | -0.05494 | 0.06946 | -0.1667 | -0.005 | 0.07206 | -0.1015 | 100.0% |
+| longest_b | 50,466 | +0.03836 | +0 | 0.3911 | +0 | +0 | 0.06432 | +0.02031 | 3.7% |
+| ExplainedIonCurrentRatio | 50,466 | -0.004182 | -1e-08 | 0.01876 | -0.02749 | +4.7e-07 | 0.004671 | -0.03985 | 19.5% |
+| CTermIonCurrentRatio | 50,466 | -0.00432 | -5.85e-09 | 0.01888 | -0.02808 | +3.7e-07 | 0.004644 | -0.05176 | 18.6% |
+| NTermIonCurrentRatio | 50,466 | +0.0001374 | +0 | 0.002249 | -1e-07 | +0.0005453 | 0.0003215 | +0.02227 | 10.2% |
+| enzC | 50,466 | +0 | +0 | 0.008903 | +0 | +0 | 7.926e-05 | -4.061e-05 | 0.0% |
+| dm | 50,466 | +1.85e-06 | +1.125e-07 | 2.801e-05 | -3.922e-05 | +4.957e-05 | 2.229e-05 | +0.0008993 | 60.7% |
+| absdm | 50,466 | +1.519e-06 | +8.8e-07 | 2.801e-05 | -4.196e-05 | +4.836e-05 | 2.228e-05 | +0.0097 | 60.7% |
+| enzN | 50,466 | +1.982e-05 | +0 | 0.004451 | +0 | +0 | 1.982e-05 | +0 | 0.0% |
+| isotope_error | 50,466 | +0 | +0 | 0 | +0 | +0 | 0 | +0 | 0.0% |
+| peplen | 50,466 | +0 | +0 | 0 | +0 | +0 | 0 | +0 | 0.0% |
+| charge2 | 50,466 | +0 | +0 | 0 | +0 | +0 | 0 | +0 | 0.0% |
+| charge3 | 50,466 | +0 | +0 | 0 | +0 | +0 | 0 | +0 | 0.0% |
+| charge4 | 50,466 | +0 | +0 | 0 | +0 | +0 | 0 | +0 | 0.0% |
+| enzInt | 50,466 | +0 | +0 | 0 | +0 | +0 | 0 | +0 | 0.0% |
+| IsolationWindowEfficiency | 50,466 | +0 | +0 | 0 | +0 | +0 | 0 | +nan | — |
+
+## Notes
+
+- Δ = (Rust value) - (Java value).
+- `mean rel Δ` = mean of (Δ / |java|) over PSMs with |java| > 1e-12.
+- `%frac |relΔ|>1%` = fraction of PSMs where the relative diff exceeds 1%.
+- Agreement bucket restricts to scans + peptides present as target on BOTH sides; this strips ranking-flip + retention-only effects so the table measures FEATURE divergence specifically.

docs/parity-analysis/diff/iter29/report.md

@@ -0,0 +1,106 @@
+# Rust↔Java PIN diff report
+
+- Java:  `/srv/data/msgf-bench/bench-merged-results/astral-java.pin` (136271 rows, 121654 scans)
+- Rust:  `/srv/data/msgf-bench/bench-iter29-results/astral-rust-iter29.pin` (149577 rows, 121677 scans)
+- Total unique scans: 121681
+
+## Top-1-per-scan buckets
+
+| Bucket | Count | % of total |
+|---|---:|---:|
+| both_target_same_peptide | 45,856 | 37.69% |
+| both_target_diff_peptide | 18,893 | 15.53% |
+| java_target_rust_decoy | 16,636 | 13.67% |
+| rust_target_java_decoy | 13,633 | 11.20% |
+| java_only_target | 3 | 0.00% |
+| rust_only_target | 14 | 0.01% |
+| both_decoy | 26,632 | 21.89% |
+| java_only_decoy | 1 | 0.00% |
+| rust_only_decoy | 13 | 0.01% |
+| both_missing | 0 | 0.00% |
+| **total** | **121,681** | 100.00% |
+
+### Sample disagreements
+
+**both_target_diff_peptide** (first 5):
+
+- `(11, 'QKYFYR', 'QSELDRR')`
+- `(13, 'TPEEGEK', 'MTYTTR')`
+- `(14, 'MLQTPESR', 'MSNRNNNK')`
+- `(15, 'PEWTQR', 'SSRSSHR')`
+- `(32, 'NEYNDR', 'MQEEEK')`
+
+**rust_target_java_decoy** (first 5):
+
+- `(20, 'EQCDTK', 'QMNDEK')`
+- `(23, 'LHNYEDMLEKNK', 'FQARCCPLQNQK')`
+- `(25, 'RYMEYR', 'PYSYNHR')`
+- `(33, 'QMAASPAK', 'SSENPLR')`
+- `(37, 'WQDFTK', 'YCAERK')`
+
+**java_target_rust_decoy** (first 5):
+
+- `(21, 'NEEQSR', 'TEAPCGK')`
+- `(29, 'MPPAGGPR', 'HARAGCK')`
+- `(45, 'RCSEDK', 'EKDCDK')`
+- `(47, 'VWFSQIEYIVLR', 'FQLLEKYEPLNR')`
+- `(48, 'VEEQEK', 'TAQEWK')`
+
+**java_only_target** (first 3):
+
+- `(13201, 'MCYGYGCGCGSFCR')`
+- `(14402, 'MDPNCSCATGGSCSCASSCKCK')`
+- `(30441, 'MDLSCSCATGGSCTCASSCK')`
+
+**rust_only_target** (first 5):
+
+- `(106652, 'DCPSCK')`
+- `(113243, 'LYYALK')`
+- `(113285, 'MKRGQR')`
+- `(114265, 'LIVLKEK')`
+- `(114482, 'YNEGCR')`
+
+## Per-feature diff (agreement bucket: same scan + peptide, both target)
+
+_50,450 PSMs in agreement bucket._
+
+Sorted by mean |Δ| (Rust - Java), descending:
+
+| Feature | n | mean Δ | median Δ | stdev | p5 | p95 | mean \|Δ\| | mean rel Δ | %frac \|relΔ\|>1% |
+|---|---:|---:|---:|---:|---:|---:|---:|---:|---:|
+| MS2IonCurrent | 50,450 | +251.6 | +0 | 7906 | -4 | +4.082 | 253 | +0.0002487 | 0.2% |
+| RawScore | 50,450 | -2.898 | -2 | 10.48 | -21 | +12 | 8.564 | +0.1287 | 96.4% |
+| lnEValue | 50,450 | -2.857 | -3.348 | 3.852 | -8.314 | +4.252 | 4.041 | -6.94 | 98.9% |
+| lnSpecEValue | 50,450 | +1.712 | +1.199 | 3.895 | -3.763 | +8.952 | 3.127 | +0.04634 | 95.5% |
+| MeanRelErrorTop7 | 50,450 | +1.344 | +1.463 | 2.737 | -3.233 | +5.477 | 2.267 | -6.28 | 99.5% |
+| NumMatchedMainIons | 50,450 | -1.295 | -1 | 2.36 | -5 | +2 | 1.961 | -0.1469 | 80.2% |
+| MeanErrorTop7 | 50,450 | -1.695 | -1.466 | 1.887 | -5.089 | +0.7928 | 1.941 | -0.2929 | 99.1% |
+| StdevRelErrorTop7 | 50,450 | -1.009 | -0.4194 | 2.67 | -6.535 | +2.392 | 1.719 | -0.08391 | 98.1% |
+| DeNovoScore | 50,450 | +0.3781 | +0 | 4.477 | -3 | +6 | 1.417 | +0.007059 | 34.3% |
+| StdevErrorTop7 | 50,450 | -0.4525 | -0.3365 | 1.7 | -3.882 | +2.141 | 1.152 | -0.07113 | 98.5% |
+| longest_y | 50,450 | -0.1275 | +0 | 0.848 | -1 | +0 | 0.2199 | -0.009755 | 9.5% |
+| matchedIonRatio | 50,450 | +0.04452 | +0.02381 | 0.2151 | -0.2848 | +0.425 | 0.1716 | +0.01379 | 96.9% |
+| lnDeltaSpecEValue | 50,450 | -0.09411 | +0 | 0.8988 | +0 | +0 | 0.09411 | +nan | — |
+| longest_y_pct | 50,450 | -0.06726 | -0.05494 | 0.06945 | -0.1667 | -0.004926 | 0.07208 | -0.1015 | 100.0% |
+| longest_b | 50,450 | +0.03845 | +0 | 0.3913 | +0 | +0 | 0.06438 | +0.02039 | 3.7% |
+| ExplainedIonCurrentRatio | 50,450 | -0.004183 | -1e-08 | 0.01876 | -0.0275 | +4.7e-07 | 0.004672 | -0.03985 | 19.5% |
+| CTermIonCurrentRatio | 50,450 | -0.004321 | -6e-09 | 0.01888 | -0.02808 | +3.7e-07 | 0.004646 | -0.05177 | 18.6% |
+| NTermIonCurrentRatio | 50,450 | +0.0001379 | +0 | 0.002249 | -1e-07 | +0.0005492 | 0.0003216 | +0.02266 | 10.2% |
+| enzC | 50,450 | +0 | +0 | 0.008904 | +0 | +0 | 7.929e-05 | -4.062e-05 | 0.0% |
+| dm | 50,450 | +1.836e-06 | +8e-08 | 2.801e-05 | -3.928e-05 | +4.959e-05 | 2.23e-05 | +0.0007948 | 60.7% |
+| absdm | 50,450 | +1.541e-06 | +9.185e-07 | 2.802e-05 | -4.196e-05 | +4.836e-05 | 2.229e-05 | +0.00948 | 60.7% |
+| enzN | 50,450 | +1.982e-05 | +0 | 0.004452 | +0 | +0 | 1.982e-05 | +0 | 0.0% |
+| isotope_error | 50,450 | +0 | +0 | 0 | +0 | +0 | 0 | +0 | 0.0% |
+| peplen | 50,450 | +0 | +0 | 0 | +0 | +0 | 0 | +0 | 0.0% |
+| charge2 | 50,450 | +0 | +0 | 0 | +0 | +0 | 0 | +0 | 0.0% |
+| charge3 | 50,450 | +0 | +0 | 0 | +0 | +0 | 0 | +0 | 0.0% |
+| charge4 | 50,450 | +0 | +0 | 0 | +0 | +0 | 0 | +0 | 0.0% |
+| enzInt | 50,450 | +0 | +0 | 0 | +0 | +0 | 0 | +0 | 0.0% |
+| IsolationWindowEfficiency | 50,450 | +0 | +0 | 0 | +0 | +0 | 0 | +nan | — |
+
+## Notes
+
+- Δ = (Rust value) - (Java value).
+- `mean rel Δ` = mean of (Δ / |java|) over PSMs with |java| > 1e-12.
+- `%frac |relΔ|>1%` = fraction of PSMs where the relative diff exceeds 1%.
+- Agreement bucket restricts to scans + peptides present as target on BOTH sides; this strips ranking-flip + retention-only effects so the table measures FEATURE divergence specifically.