Skills › Scientific Agent Skills › pyopenms › process_identifications.py
process_identifications.py
148 lines · skills/pyopenms/scripts/process_identifications.py at 5f2bc03
1#!/usr/bin/env python3
2"""
3Process Peptide/Protein Identifications
5Post-process search-engine results (idXML): optionally re-index peptides against
6a protein FASTA (assigns target/decoy + protein accessions), estimate FDR/q-values,
7filter by FDR threshold, peptide length, and best-hit-per-spectrum, then export a
8filtered idXML and a flat CSV of peptide hits.
10Usage:
11 python process_identifications.py search.idXML --out filtered.idXML --csv hits.csv
12 python process_identifications.py search.idXML --fasta db.fasta --fdr 0.01 --out filt.idXML
13 python process_identifications.py search.idXML --fdr 0.05 --min-length 7 --best-per-spectrum
14"""
16import argparse
17import os
18import math
19import sys
21try:
22 import pyopenms as ms
23except ImportError:
24 print("Error: pyopenms not installed. Install with: uv pip install pyopenms")
25 sys.exit(1)
28def index_peptides(fasta, prot_ids, pep_ids, decoy_string="DECOY_"):
29 fasta_entries = []
30 ms.FASTAFile().load(fasta, fasta_entries)
31 indexer = ms.PeptideIndexing()
32 p = indexer.getParameters()
33 p.setValue("decoy_string", decoy_string)
34 p.setValue("missing_decoy_action", "error")
35 indexer.setParameters(p)
36 result = indexer.run(fasta_entries, prot_ids, pep_ids)
37 if result[0] != ms.PeptideIndexing.ExitCodes.EXECUTION_OK:
38 raise ValueError(f"Peptide indexing failed: {result[0]}")
39 print(f"Indexed against {len(fasta_entries)} proteins")
43def filter_psm_qvalues(pep_ids, threshold):
44 """Target-decoy q-values for comparable top-hit PSMs from one search run."""
45 if not math.isfinite(threshold) or not 0 <= threshold <= 1:
46 raise ValueError("q-value threshold must be finite and between 0 and 1")
47 contexts = set()
48 labels = set()
49 for pid in pep_ids:
50 hits = pid.getHits()
51 if not hits:
52 continue
53 if not pid.getScoreType() or pid.getScoreType().lower() in {"q-value", "fdr"}:
54 raise ValueError("FDR estimation needs original search scores, not existing q-values/FDR")
55 contexts.add((pid.getIdentifier(), pid.getScoreType(), pid.isHigherScoreBetter()))
56 for hit in hits:
57 if not math.isfinite(hit.getScore()) or not hit.metaValueExists("target_decoy"):
58 raise ValueError("Every hit needs a finite score and target_decoy annotation")
59 if hit.getMetaValue("target_decoy") not in {"target", "decoy", "target+decoy"}:
60 raise ValueError("Invalid target_decoy annotation")
61 best = max(hits, key=lambda h: h.getScore()) if pid.isHigherScoreBetter() else min(hits, key=lambda h: h.getScore())
62 labels.add(best.getMetaValue("target_decoy"))
63 if len(contexts) != 1:
64 raise ValueError("Process one search run with one score type/direction at a time")
65 if "decoy" not in labels or not labels.intersection({"target", "target+decoy"}):
66 raise ValueError("FDR estimation requires both target and decoy top-scoring PSMs")
67 fdr = ms.FalseDiscoveryRate()
68 p = fdr.getParameters()
69 p.setValue("no_qvalues", "false")
70 p.setValue("use_all_hits", "false")
71 p.setValue("add_decoy_peptides", "false")
72 fdr.setParameters(p)
73 fdr.apply(pep_ids)
74 ms.IDFilter().filterHitsByScore(pep_ids, threshold)
77def main():
78 parser = argparse.ArgumentParser(description="Filter and export peptide identifications.")
79 parser.add_argument("input", help="Input idXML")
80 parser.add_argument("--decoy-prefix", default="DECOY_", help="FASTA decoy accession prefix for re-indexing")
81 parser.add_argument("--fasta", help="Protein FASTA (target+decoy) for re-indexing")
82 parser.add_argument("--out", help="Output filtered idXML")
83 parser.add_argument("--csv", help="Output CSV of peptide hits")
84 parser.add_argument("--fdr", type=float, help="FDR/q-value threshold (requires decoys)")
85 parser.add_argument("--min-length", type=int, help="Minimum peptide length")
86 parser.add_argument("--max-length", type=int, help="Maximum peptide length")
87 parser.add_argument("--best-per-spectrum", action="store_true",
88 help="Keep only the best hit per spectrum")
89 args = parser.parse_args()
91 if not os.path.exists(args.input):
92 print(f"Error: file not found: {args.input}")
93 sys.exit(1)
95 # In pyOpenMS 3.5+, idXML load/store require a PeptideIdentificationList
96 # (not a plain Python list) for the peptide IDs; protein IDs stay a list.
97 prot_ids = []
98 pep_ids = ms.PeptideIdentificationList()
99 ms.IdXMLFile().load(args.input, prot_ids, pep_ids)
100 print(f"Loaded {len(pep_ids)} spectra, "
101 f"{sum(len(p.getHits()) for p in pep_ids)} peptide hits")
103 if args.fasta:
104 index_peptides(args.fasta, prot_ids, pep_ids, args.decoy_prefix)
106 if args.fdr is not None:
107 filter_psm_qvalues(pep_ids, args.fdr)
108 print(f"Applied FDR filter (q <= {args.fdr})")
110 if args.best_per_spectrum:
111 ms.IDFilter().keepBestPeptideHits(pep_ids, False)
112 print("Kept best hit per spectrum")
114 if args.min_length is not None or args.max_length is not None:
115 lo = args.min_length or 0
116 hi = args.max_length or 0 # 0 = no upper bound in OpenMS
117 ms.IDFilter().filterPeptidesByLength(pep_ids, lo, hi)
118 print(f"Length filter ({lo}-{hi or 'inf'})")
120 # Drop now-empty identifications
121 ms.IDFilter().removeEmptyIdentifications(pep_ids)
122 remaining = sum(len(p.getHits()) for p in pep_ids)
123 print(f"Remaining: {len(pep_ids)} spectra, {remaining} peptide hits")
125 if args.out:
126 ms.IdXMLFile().store(args.out, prot_ids, pep_ids)
127 print(f"Wrote {args.out}")
129 if args.csv:
130 import csv
131 with open(args.csv, "w", newline="") as fh:
132 w = csv.writer(fh)
133 w.writerow(["rt", "mz", "sequence", "charge", "score", "score_type",
134 "target_decoy", "accessions"])
135 for pid in pep_ids:
136 rt = pid.getRT()
137 mz = pid.getMZ()
138 st = pid.getScoreType()
139 for hit in pid.getHits():
140 td = hit.getMetaValue("target_decoy") if hit.metaValueExists("target_decoy") else ""
141 accs = ";".join(a.decode() if isinstance(a, bytes) else a for a in hit.extractProteinAccessionsSet())
142 w.writerow([f"{rt:.2f}", f"{mz:.4f}", hit.getSequence().toString(),
143 hit.getCharge(), hit.getScore(), st, td, accs])
144 print(f"Wrote {args.csv}")
147if __name__ == "__main__":
148 main()