Skills › Scientific Agent Skills › pyopenms › digest_protein.py
digest_protein.py
104 lines · skills/pyopenms/scripts/digest_protein.py at 5f2bc03
1#!/usr/bin/env python3
2"""
3In-Silico Protein Digestion
5Digest protein sequences (FASTA or a single sequence) with a configurable
6protease, producing theoretical peptides with masses and m/z. Useful for
7targeted method design and search-space estimation.
9Usage:
10 python digest_protein.py proteins.fasta --out peptides.csv
11 python digest_protein.py --sequence MKWVTFISLLLLFSSAYS --enzyme Trypsin --missed 2
12 python digest_protein.py proteins.fasta --min-length 7 --max-length 40 --charges 1 2 3
13"""
15import argparse
16import os
17import sys
19try:
20 import pyopenms as ms
21except ImportError:
22 print("Error: pyopenms not installed. Install with: uv pip install pyopenms")
23 sys.exit(1)
25PROTON = 1.0072764665789
28def read_fasta(path):
29 entries = []
30 fe = ms.FASTAFile()
31 seqs = []
32 fe.load(path, seqs)
33 for s in seqs:
34 entries.append((s.identifier, s.sequence))
35 return entries
38def digest(seq_str, enzyme, missed, min_len, max_len):
39 dig = ms.ProteaseDigestion()
40 dig.setEnzyme(enzyme)
41 dig.setMissedCleavages(missed)
42 out = []
43 dig.digest(ms.AASequence.fromString(seq_str), out, min_len, max_len)
44 return out
47def main():
48 parser = argparse.ArgumentParser(description="In-silico protein digestion.")
49 parser.add_argument("fasta", nargs="?", help="FASTA file (optional if --sequence given)")
50 parser.add_argument("--sequence", help="Single protein sequence")
51 parser.add_argument("--enzyme", default="Trypsin", help="Protease (default Trypsin)")
52 parser.add_argument("--missed", type=int, default=2, help="Max missed cleavages (default 2)")
53 parser.add_argument("--min-length", type=int, default=6, help="Min peptide length (default 6)")
54 parser.add_argument("--max-length", type=int, default=40, help="Max peptide length (default 40)")
55 parser.add_argument("--charges", type=int, nargs="+", default=[1, 2], help="m/z charge states")
56 parser.add_argument("--out", help="Output CSV of peptides")
57 args = parser.parse_args()
58 if any(z <= 0 for z in args.charges):
59 parser.error("--charges must contain positive charge magnitudes")
61 proteins = []
62 if args.fasta:
63 if not os.path.exists(args.fasta):
64 print(f"Error: file not found: {args.fasta}")
65 sys.exit(1)
66 proteins = read_fasta(args.fasta)
67 elif args.sequence:
68 proteins = [("input", args.sequence)]
69 else:
70 parser.error("provide a FASTA file or --sequence")
72 rows = []
73 seen = set()
74 for prot_id, seq in proteins:
75 for pep in digest(seq, args.enzyme, args.missed, args.min_length, args.max_length):
76 pep_str = pep.toString()
77 key = (prot_id, pep_str)
78 if key in seen:
79 continue
80 seen.add(key)
81 mono = pep.getMonoWeight()
82 row = {"protein": prot_id, "peptide": pep_str, "length": pep.size(),
83 "mono_mass": round(mono, 5)}
84 for z in args.charges:
85 row[f"mz_z{z}"] = round((mono + z * PROTON) / z, 5)
86 rows.append(row)
88 print(f"Proteins: {len(proteins)} Unique peptides: {len(rows)}")
89 for r in rows[:10]:
90 print(f" {r['peptide']} ({r['length']} aa, {r['mono_mass']} Da)")
91 if len(rows) > 10:
92 print(f" ... and {len(rows) - 10} more")
94 if args.out:
95 import csv
96 with open(args.out, "w", newline="") as fh:
97 w = csv.DictWriter(fh, fieldnames=list(rows[0].keys()) if rows else ["protein", "peptide"])
98 w.writeheader()
99 w.writerows(rows)
100 print(f"Wrote {args.out}")
103if __name__ == "__main__":
104 main()