From ae3ac01628c40d9fd7621f123ef4fce687e2f12d Mon Sep 17 00:00:00 2001 From: Longping Fu Date: Tue, 4 Aug 2026 20:35:48 -0400 Subject: [PATCH 1/2] Route all spectrum reading through open_spectra; unpin AlphaFold models Four global rules were violated across the active analysis scripts. None changes any output today; each removes a way the scripts break tomorrow. pyteomics.mzml (4 scripts) streamed whole calibrated mzML files to pull a handful of scans, re-deriving precursor m/z, filter string and RT from raw dict paths. Replaced with indexed reader.get_spectrum(). The two extract_spectrum_data() helpers keep their contract exactly: None for a missing scan, and precursor keys left absent for MS1 -- get_spectrum reports precursor_mz as 0.0 rather than None there, so the guard is a truthiness check, matching what callers' .get(..., 0) fallbacks expect. Direct MzMLReader construction (9 scripts, 13 sites) became open_spectra, the documented drop-in. With no spectra cache built yet it returns an MzMLReader, so behaviour is identical now and the cache speedup arrives for free once one exists. sys.path.insert of /Users/longpingfu/Downloads/... (4 scripts) pinned this machine's layout for two packages that are pip-installed editable; verified both import without them. Removed, along with the imports left unused. Hard-coded AF-{acc}-F1-model_v6.pdb paths (4 PyMOL scripts) now resolve the newest cached model, the pattern Figure6F_EWSR1_S274_pymol.py already uses -- that script was pinned to model_v4 while the rest of the figure used v6. All three accessions resolve to the same v6 file they named, so the panels render unchanged. Sort is numeric on the version, not lexical. Verified: every active script compiles; open_spectra and both rewritten helpers exercised against a real calibrated mzML (MS1, HCD, EThcD and a missing scan); each PyMOL resolver block executed as it sits in the file. Co-Authored-By: Claude Opus 5 (1M context) --- data_analysis/Figure6F_HOXA13_pymol.py | 22 ++++++++- data_analysis/Figure6F_HYOU1_pymol.py | 22 ++++++++- data_analysis/Figure6F_pymol.py | 22 ++++++++- data_analysis/Figure6F_ray_modes.py | 22 ++++++++- data_analysis/annotate_edem1_spectra.py | 2 +- .../annotate_ethcd_secretory_sites.py | 2 +- data_analysis/annotate_figure6d_all_psms.py | 2 +- data_analysis/annotate_figure6d_spectra.py | 4 +- .../annotate_low_prob_hcd_spectra.py | 43 ++++++----------- .../annotate_nonsecretory_ogalnac.py | 48 +++++-------------- data_analysis/batch_annotate_spectra.py | 33 +++++-------- .../create_spring_spectrum_candidates_v2.py | 16 +++---- data_analysis/fix_ogalnac_site_quant.py | 2 +- data_analysis/generate_ogalnac_spectra.py | 2 +- data_analysis/generate_selected_spectra.py | 4 -- data_analysis/regenerate_ogalnac_spectra.py | 2 +- data_analysis/regenerate_tyr_spectra_v5.py | 2 +- data_analysis/regenerate_with_localization.py | 29 ++++------- 18 files changed, 145 insertions(+), 134 deletions(-) diff --git a/data_analysis/Figure6F_HOXA13_pymol.py b/data_analysis/Figure6F_HOXA13_pymol.py index b23f9b3..de54fe0 100644 --- a/data_analysis/Figure6F_HOXA13_pymol.py +++ b/data_analysis/Figure6F_HOXA13_pymol.py @@ -19,7 +19,27 @@ # Configuration # ============================================ -STRUCTURE_PATH = "/Volumes/cos-lab-rwu60/Longping/OGlycoTM_Final_Version/data_source/alphafold_structures/AF-P31271-F1-model_v6.pdb" +STRUCTURE_DIR = "/Volumes/cos-lab-rwu60/Longping/OGlycoTM_Final_Version/data_source/alphafold_structures" +PROTEIN_ID = "P31271" + + +def resolve_structure(protein_id): + """Newest cached AlphaFold model for this accession. + + Never hand-write an AlphaFold filename or version: v4/v5 files are deleted upstream and v6 + abolished the -F2-/-F3- fragment scheme, so a pinned `model_v6` path is a latent breakage. + Take the newest cached file, or resolve one with mzml_utils.structure.fetch_structure() + (as data_analysis/pymol_site_panels.py does). + """ + import glob + import re + hits = glob.glob(os.path.join(STRUCTURE_DIR, f"AF-{protein_id}-F1-model_v*.pdb")) + if not hits: + raise SystemExit(f"no cached AlphaFold model for {protein_id} in {STRUCTURE_DIR}") + return max(hits, key=lambda p: int(re.search(r"model_v(\d+)", p).group(1))) + + +STRUCTURE_PATH = resolve_structure(PROTEIN_ID) OUTPUT_DIR = "/Volumes/cos-lab-rwu60/Longping/OGlycoTM_Final_Version/Figures/Figure6" # O-GlcNAc site data (HEK293T) - IDR REGION EXAMPLE diff --git a/data_analysis/Figure6F_HYOU1_pymol.py b/data_analysis/Figure6F_HYOU1_pymol.py index 086dc32..bc41613 100644 --- a/data_analysis/Figure6F_HYOU1_pymol.py +++ b/data_analysis/Figure6F_HYOU1_pymol.py @@ -17,7 +17,27 @@ # Configuration # ============================================ -STRUCTURE_PATH = "/Volumes/cos-lab-rwu60/Longping/OGlycoTM_Final_Version/data_source/alphafold_structures/AF-Q9Y4L1-F1-model_v6.pdb" +STRUCTURE_DIR = "/Volumes/cos-lab-rwu60/Longping/OGlycoTM_Final_Version/data_source/alphafold_structures" +PROTEIN_ID = "Q9Y4L1" + + +def resolve_structure(protein_id): + """Newest cached AlphaFold model for this accession. + + Never hand-write an AlphaFold filename or version: v4/v5 files are deleted upstream and v6 + abolished the -F2-/-F3- fragment scheme, so a pinned `model_v6` path is a latent breakage. + Take the newest cached file, or resolve one with mzml_utils.structure.fetch_structure() + (as data_analysis/pymol_site_panels.py does). + """ + import glob + import re + hits = glob.glob(os.path.join(STRUCTURE_DIR, f"AF-{protein_id}-F1-model_v*.pdb")) + if not hits: + raise SystemExit(f"no cached AlphaFold model for {protein_id} in {STRUCTURE_DIR}") + return max(hits, key=lambda p: int(re.search(r"model_v(\d+)", p).group(1))) + + +STRUCTURE_PATH = resolve_structure(PROTEIN_ID) OUTPUT_DIR = "/Volumes/cos-lab-rwu60/Longping/OGlycoTM_Final_Version/Figures/Figure6" # O-GlcNAc site data (HEK293T) - STRUCTURED REGION EXAMPLE diff --git a/data_analysis/Figure6F_pymol.py b/data_analysis/Figure6F_pymol.py index 8c95ec0..ae124b7 100644 --- a/data_analysis/Figure6F_pymol.py +++ b/data_analysis/Figure6F_pymol.py @@ -18,7 +18,27 @@ # ============================================ # File paths -STRUCTURE_PATH = "/Volumes/cos-lab-rwu60/Longping/OGlycoTM_Final_Version/data_source/alphafold_structures/AF-Q9NYJ8-F1-model_v6.pdb" +STRUCTURE_DIR = "/Volumes/cos-lab-rwu60/Longping/OGlycoTM_Final_Version/data_source/alphafold_structures" +PROTEIN_ID = "Q9NYJ8" + + +def resolve_structure(protein_id): + """Newest cached AlphaFold model for this accession. + + Never hand-write an AlphaFold filename or version: v4/v5 files are deleted upstream and v6 + abolished the -F2-/-F3- fragment scheme, so a pinned `model_v6` path is a latent breakage. + Take the newest cached file, or resolve one with mzml_utils.structure.fetch_structure() + (as data_analysis/pymol_site_panels.py does). + """ + import glob + import re + hits = glob.glob(os.path.join(STRUCTURE_DIR, f"AF-{protein_id}-F1-model_v*.pdb")) + if not hits: + raise SystemExit(f"no cached AlphaFold model for {protein_id} in {STRUCTURE_DIR}") + return max(hits, key=lambda p: int(re.search(r"model_v(\d+)", p).group(1))) + + +STRUCTURE_PATH = resolve_structure(PROTEIN_ID) OUTPUT_DIR = "/Volumes/cos-lab-rwu60/Longping/OGlycoTM_Final_Version/Figures/Figure6" # O-GlcNAc site data (HEK293T) diff --git a/data_analysis/Figure6F_ray_modes.py b/data_analysis/Figure6F_ray_modes.py index a6fdc72..bd42f23 100644 --- a/data_analysis/Figure6F_ray_modes.py +++ b/data_analysis/Figure6F_ray_modes.py @@ -10,7 +10,27 @@ # Configuration # ============================================ -STRUCTURE_PATH = "/Volumes/cos-lab-rwu60/Longping/OGlycoTM_Final_Version/data_source/alphafold_structures/AF-Q9NYJ8-F1-model_v6.pdb" +STRUCTURE_DIR = "/Volumes/cos-lab-rwu60/Longping/OGlycoTM_Final_Version/data_source/alphafold_structures" +PROTEIN_ID = "Q9NYJ8" + + +def resolve_structure(protein_id): + """Newest cached AlphaFold model for this accession. + + Never hand-write an AlphaFold filename or version: v4/v5 files are deleted upstream and v6 + abolished the -F2-/-F3- fragment scheme, so a pinned `model_v6` path is a latent breakage. + Take the newest cached file, or resolve one with mzml_utils.structure.fetch_structure() + (as data_analysis/pymol_site_panels.py does). + """ + import glob + import re + hits = glob.glob(os.path.join(STRUCTURE_DIR, f"AF-{protein_id}-F1-model_v*.pdb")) + if not hits: + raise SystemExit(f"no cached AlphaFold model for {protein_id} in {STRUCTURE_DIR}") + return max(hits, key=lambda p: int(re.search(r"model_v(\d+)", p).group(1))) + + +STRUCTURE_PATH = resolve_structure(PROTEIN_ID) OUTPUT_DIR = "/Volumes/cos-lab-rwu60/Longping/OGlycoTM_Final_Version/Figures/Figure6" # O-GlcNAc site data (HEK293T) diff --git a/data_analysis/annotate_edem1_spectra.py b/data_analysis/annotate_edem1_spectra.py index 7a2abd4..68d84d2 100644 --- a/data_analysis/annotate_edem1_spectra.py +++ b/data_analysis/annotate_edem1_spectra.py @@ -56,7 +56,7 @@ if mzml_path not in readers: print(f" Opening {os.path.basename(mzml_path)}...") - readers[mzml_path] = mzml_utils.MzMLReader(mzml_path) + readers[mzml_path] = mzml_utils.open_spectra(mzml_path) reader = readers[mzml_path] # Get spectrum diff --git a/data_analysis/annotate_ethcd_secretory_sites.py b/data_analysis/annotate_ethcd_secretory_sites.py index c079616..9651939 100644 --- a/data_analysis/annotate_ethcd_secretory_sites.py +++ b/data_analysis/annotate_ethcd_secretory_sites.py @@ -202,7 +202,7 @@ def main(): if not os.path.exists(cal_mzml): print(f' SKIP: no mzML for {raw_file}') continue - readers[raw_file] = mzml_utils.MzMLReader(cal_mzml) + readers[raw_file] = mzml_utils.open_spectra(cal_mzml) reader = readers[raw_file] # 1. MS1 isolation diff --git a/data_analysis/annotate_figure6d_all_psms.py b/data_analysis/annotate_figure6d_all_psms.py index 16afe9d..72a76cc 100644 --- a/data_analysis/annotate_figure6d_all_psms.py +++ b/data_analysis/annotate_figure6d_all_psms.py @@ -203,7 +203,7 @@ def main(): if not os.path.exists(cal_mzml): print(f' SKIP: no mzML for {raw_file}') continue - readers[raw_file] = mzml_utils.MzMLReader(cal_mzml) + readers[raw_file] = mzml_utils.open_spectra(cal_mzml) reader = readers[raw_file] # 1. MS1 isolation diff --git a/data_analysis/annotate_figure6d_spectra.py b/data_analysis/annotate_figure6d_spectra.py index 40efac6..8f61fb4 100644 --- a/data_analysis/annotate_figure6d_spectra.py +++ b/data_analysis/annotate_figure6d_spectra.py @@ -258,7 +258,7 @@ def main(): if not os.path.exists(cal_mzml): print(f' SKIP: no mzML for {raw_file}') continue - readers[raw_file] = mzml_utils.MzMLReader(cal_mzml) + readers[raw_file] = mzml_utils.open_spectra(cal_mzml) reader = readers[raw_file] # 1. MS1 isolation check @@ -316,7 +316,7 @@ def main(): if alt_raw not in readers: cal_mzml = os.path.join(MZML_DIR, f'{alt_raw}_calibrated.mzML') if os.path.exists(cal_mzml): - readers[alt_raw] = mzml_utils.MzMLReader(cal_mzml) + readers[alt_raw] = mzml_utils.open_spectra(cal_mzml) if alt_raw in readers: alt_path = os.path.join(folder, f'{site_id}_s{alt_scan}_{alt_psm["Confidence.Level"]}_{alt_act}.pdf') if not os.path.exists(alt_path): diff --git a/data_analysis/annotate_low_prob_hcd_spectra.py b/data_analysis/annotate_low_prob_hcd_spectra.py index fb8821e..99ac867 100644 --- a/data_analysis/annotate_low_prob_hcd_spectra.py +++ b/data_analysis/annotate_low_prob_hcd_spectra.py @@ -17,18 +17,14 @@ """ import os -import sys import json import numpy as np import pandas as pd -from pyteomics import mzml from collections import defaultdict import matplotlib.pyplot as plt from matplotlib.backends.backend_pdf import PdfPages -# Add GlycoSpectrumAnnotator to path -sys.path.insert(0, '/Users/longpingfu/Downloads/GlycoSpectrumAnnotator') - +from mzml_utils import open_spectra from spectrum_annotator_ddzby import ( SpectrumAnnotator, FragmentCalculator, @@ -107,36 +103,25 @@ def find_calibrated_mzml(file_name, mzml_dir): return None -def extract_spectrum_data(mzml_reader, scan_number): +def extract_spectrum_data(reader, scan_number): """ Extract complete spectrum data for a specific scan using indexed access. """ - scan_id = f"controllerType=0 controllerNumber=1 scan={scan_number}" - - try: - spectrum = mzml_reader.get_by_id(scan_id) - except KeyError: - # Try to find by scan number in ID - for spec_id in mzml_reader.index.keys(): - if f"scan={scan_number}" in spec_id: - spectrum = mzml_reader.get_by_id(spec_id) - break - else: - return None + spectrum = reader.get_spectrum(int(scan_number)) + if spectrum is None: + return None result = { - 'mz_array': spectrum.get('m/z array', np.array([])), - 'intensity_array': spectrum.get('intensity array', np.array([])), - 'ms_level': spectrum.get('ms level'), + 'mz_array': spectrum.mz, + 'intensity_array': spectrum.intensity, + 'ms_level': spectrum.ms_level, } - # Get precursor info - if 'precursorList' in spectrum: - precursor = spectrum['precursorList']['precursor'][0] - if 'selectedIonList' in precursor: - sel_ion = precursor['selectedIonList']['selectedIon'][0] - result['precursor_mz'] = sel_ion.get('selected ion m/z') - result['precursor_charge'] = sel_ion.get('charge state') + # MS1 reports precursor_mz as 0.0; leave both keys absent there, as the old + # pyteomics path did, since callers fall back with .get(..., 0) + if spectrum.precursor_mz: + result['precursor_mz'] = spectrum.precursor_mz + result['precursor_charge'] = spectrum.precursor_charge return result @@ -284,7 +269,7 @@ def main(): print(f" Processing: {os.path.basename(mzml_path)} ({len(scans)} scans)...", end=" ", flush=True) try: - with mzml.MzML(mzml_path, use_index=True) as reader: + with open_spectra(mzml_path) as reader: scans_processed = 0 for idx, scan_number, row in scans: diff --git a/data_analysis/annotate_nonsecretory_ogalnac.py b/data_analysis/annotate_nonsecretory_ogalnac.py index 037c1cf..e0bdd8f 100644 --- a/data_analysis/annotate_nonsecretory_ogalnac.py +++ b/data_analysis/annotate_nonsecretory_ogalnac.py @@ -15,8 +15,7 @@ import os import traceback -from pyteomics import mzml as pyteomics_mzml -from mzml_utils import PROTON, NEUTRON_MASS +from mzml_utils import PROTON, NEUTRON_MASS, open_spectra from spectrum_annotator_ddzby import SpectrumAnnotator DATA_BASE = '/Volumes/cos-lab-rwu60/Longping/OGlycoTM_Final_Version' @@ -181,45 +180,20 @@ def is_secretory(acc): for offset in range(-20, 5): # MS1 up to 20 before, HCD up to 4 before need_scans.add(s + offset) - # Stream through mzML once, collect needed spectra - print(f'Streaming {os.path.basename(mzml_path)} ({len(group)} PSMs)...', flush=True) + # Indexed lookup of just the scans we need -- no full-file stream + print(f'Reading {os.path.basename(mzml_path)} ({len(group)} PSMs)...', flush=True) spectra = {} - reader = pyteomics_mzml.MzML(mzml_path) - for spec in reader: - sid = spec.get('id', '') - if 'scan=' not in sid: - continue - scan = int(sid.split('scan=')[-1]) - if scan in need_scans: - mz_array = spec.get('m/z array', np.array([])) - int_array = spec.get('intensity array', np.array([])) - ms_level = spec.get('ms level', 0) - fs = spec.get('filter string', spec.get('scanList', {}).get('scan', [{}])[0].get('filter string', '')) - - precursor_mz = None - if 'precursorList' in spec: - precs = spec['precursorList'].get('precursor', []) - if precs: - ions = precs[0].get('selectedIonList', {}).get('selectedIon', []) - if ions: - precursor_mz = ions[0].get('selected ion m/z') - elif 'selected precursors' in spec: - precs = spec['selected precursors'] - if precs: - precursor_mz = precs[0].get('selected ion m/z') - - rt = spec.get('scanList', {}).get('scan', [{}])[0].get('scan start time', 0) - + with open_spectra(mzml_path) as reader: + for scan in sorted(need_scans): + spec = reader.get_spectrum(int(scan)) + if spec is None: + continue spectra[scan] = { - 'mz': mz_array, 'intensity': int_array, - 'ms_level': ms_level, 'filter_string': str(fs), - 'precursor_mz': precursor_mz, 'rt': rt, + 'mz': spec.mz, 'intensity': spec.intensity, + 'ms_level': spec.ms_level, 'filter_string': str(spec.filter_string or ''), + 'precursor_mz': spec.precursor_mz, 'rt': spec.rt, } - # Stop early if we have all needed scans - if scan > max(need_scans) + 100: - break - print(f' Collected {len(spectra)} spectra', flush=True) for _, psm in group.iterrows(): diff --git a/data_analysis/batch_annotate_spectra.py b/data_analysis/batch_annotate_spectra.py index de6c259..15acbed 100644 --- a/data_analysis/batch_annotate_spectra.py +++ b/data_analysis/batch_annotate_spectra.py @@ -6,7 +6,6 @@ Groups PSMs by raw file to read each calibrated mzML only once. """ -import numpy as np import pandas as pd import matplotlib matplotlib.use('Agg') @@ -17,7 +16,8 @@ import time import traceback from pathlib import Path -from pyteomics import mzml + +import mzml_utils from spectrum_annotator_ddzby import SpectrumAnnotator @@ -101,27 +101,20 @@ def annotate_file(raw_file, cell_type, psm_group, output_base): scan_nums = set(psm_group['scan_num'].tolist()) - # Read spectrum data from calibrated mzML + # Read spectrum data from calibrated mzML -- indexed lookup, not a full stream spectra = {} - reader = mzml.MzML(cal_mzml) - for spec in reader: - sid = spec.get('id', '') - if 'scan=' not in sid: - continue - scan = int(sid.split('scan=')[-1]) - if scan in scan_nums: - # Look up activation from pre-built cache - act_type = ACTIVATION_CACHE.get((raw_file, scan), 'HCD') - + with mzml_utils.open_spectra(cal_mzml) as reader: + for scan in scan_nums: + spec = reader.get_spectrum(int(scan)) + if spec is None: + continue spectra[scan] = { - 'mz': spec['m/z array'], - 'intensity': spec['intensity array'], - 'ms_level': spec.get('ms level', 2), - 'activation': act_type, + 'mz': spec.mz, + 'intensity': spec.intensity, + 'ms_level': spec.ms_level, + # activation from the pre-built cache, not re-derived per scan + 'activation': ACTIVATION_CACHE.get((raw_file, scan), 'HCD'), } - if len(spectra) == len(scan_nums): - break - reader.close() # Annotate each PSM n_annotated = 0 diff --git a/data_analysis/create_spring_spectrum_candidates_v2.py b/data_analysis/create_spring_spectrum_candidates_v2.py index 795cf0f..049158e 100644 --- a/data_analysis/create_spring_spectrum_candidates_v2.py +++ b/data_analysis/create_spring_spectrum_candidates_v2.py @@ -11,17 +11,13 @@ import os import shutil -import sys from pathlib import Path import matplotlib.pyplot as plt import pandas as pd -sys.path.insert(0, "/Users/longpingfu/Downloads/GlycoSpectrumAnnotator") -sys.path.insert(0, "/Users/longpingfu/Downloads/mzml-utils/src") - -from mzml_utils import MzMLReader # noqa: E402 -from spectrum_annotator_ddzby import ( # noqa: E402 +from mzml_utils import open_spectra +from spectrum_annotator_ddzby import ( SpectrumAnnotator, parse_modifications_from_string, ) @@ -65,7 +61,7 @@ def raw_file_for_ethcd(cache: pd.DataFrame, cell_type: str, ethcd_scan: int) -> return str(hit.iloc[0]["raw_file"]) -def find_precursor_matched_hcd(reader: MzMLReader, ethcd_scan: int, old_hcd_scan: int) -> tuple[int, str]: +def find_precursor_matched_hcd(reader, ethcd_scan: int, old_hcd_scan: int) -> tuple[int, str]: ethcd = reader.get_spectrum(int(ethcd_scan)) old = reader.get_spectrum(int(old_hcd_scan)) if ( @@ -138,7 +134,7 @@ def bond_coverage(matched, peptide_len: int, ion_types: set[str]) -> tuple[int, return len(bonds), f"{len(bonds)}/{peptide_len - 1}", sorted(labels) -def annotate_hcd(row: pd.Series, reader: MzMLReader, hcd_scan: int, mzml_name: str, out_path: Path) -> dict: +def annotate_hcd(row: pd.Series, reader, hcd_scan: int, mzml_name: str, out_path: Path) -> dict: spec = reader.get_spectrum(int(hcd_scan)) mods = parse_modifications_from_string(str(row["assigned_modifications"])) ann = SpectrumAnnotator( @@ -197,7 +193,7 @@ def main() -> None: "", ] - readers: dict[Path, MzMLReader] = {} + readers: dict[Path, object] = {} # open_spectra returns MzMLReader or SpectrumCache try: for _, row in metrics.iterrows(): old_dir = Path(row["candidate_dir"]) @@ -208,7 +204,7 @@ def main() -> None: raw_file = raw_file_for_ethcd(cache, str(row["cell_type"]), int(row["ethcd_scan"])) mzml_path = calibrated_mzml_path(str(row["cell_type"]), raw_file) if mzml_path not in readers: - readers[mzml_path] = MzMLReader(str(mzml_path)) + readers[mzml_path] = open_spectra(str(mzml_path)) reader = readers[mzml_path] checked_hcd, status = find_precursor_matched_hcd( diff --git a/data_analysis/fix_ogalnac_site_quant.py b/data_analysis/fix_ogalnac_site_quant.py index 8a24ac0..970a347 100644 --- a/data_analysis/fix_ogalnac_site_quant.py +++ b/data_analysis/fix_ogalnac_site_quant.py @@ -139,7 +139,7 @@ def main(): for suffix in ['_calibrated.mzML', '_mz_calibrated.mzML']: path = os.path.join(CELL_MZML_DIRS[cell], f'{raw}{suffix}') if os.path.exists(path): - readers[raw] = mzml_utils.MzMLReader(path) + readers[raw] = mzml_utils.open_spectra(path) break if raw in readers: diff --git a/data_analysis/generate_ogalnac_spectra.py b/data_analysis/generate_ogalnac_spectra.py index 941a2bb..003db3b 100644 --- a/data_analysis/generate_ogalnac_spectra.py +++ b/data_analysis/generate_ogalnac_spectra.py @@ -140,7 +140,7 @@ def main(): n_psms = len(psm_group) print(f'[{i+1}/{n_files}] {raw_file} ({cell_type}): {n_psms} PSMs...', end=' ', flush=True) - reader = mzml_utils.MzMLReader(cal_mzml) + reader = mzml_utils.open_spectra(cal_mzml) n_ok = 0 for _, psm in psm_group.iterrows(): diff --git a/data_analysis/generate_selected_spectra.py b/data_analysis/generate_selected_spectra.py index 4b8c28e..aa25253 100644 --- a/data_analysis/generate_selected_spectra.py +++ b/data_analysis/generate_selected_spectra.py @@ -15,7 +15,6 @@ """ import os -import sys import json import subprocess import numpy as np @@ -23,9 +22,6 @@ import matplotlib.pyplot as plt from matplotlib.backends.backend_pdf import PdfPages -# Add GlycoSpectrumAnnotator to path for updated annotator -sys.path.insert(0, '/Users/longpingfu/Downloads/GlycoSpectrumAnnotator') - from spectrum_annotator_ddzby import ( SpectrumAnnotator, FragmentCalculator, diff --git a/data_analysis/regenerate_ogalnac_spectra.py b/data_analysis/regenerate_ogalnac_spectra.py index 4e353ba..cd26f22 100644 --- a/data_analysis/regenerate_ogalnac_spectra.py +++ b/data_analysis/regenerate_ogalnac_spectra.py @@ -243,7 +243,7 @@ def main(): mzml_path = os.path.join(MZML_DIR, f'{raw_file}_calibrated.mzML') if mzml_path not in readers: print(f' Opening {os.path.basename(mzml_path)}...') - readers[mzml_path] = mzml_utils.MzMLReader(mzml_path) + readers[mzml_path] = mzml_utils.open_spectra(mzml_path) reader = readers[mzml_path] row = { diff --git a/data_analysis/regenerate_tyr_spectra_v5.py b/data_analysis/regenerate_tyr_spectra_v5.py index f6a7658..6cc8758 100644 --- a/data_analysis/regenerate_tyr_spectra_v5.py +++ b/data_analysis/regenerate_tyr_spectra_v5.py @@ -283,7 +283,7 @@ def main(): if mzml_path not in readers: print(f' Opening {os.path.basename(mzml_path)}...') - readers[mzml_path] = mzml_utils.MzMLReader(mzml_path) + readers[mzml_path] = mzml_utils.open_spectra(mzml_path) reader = readers[mzml_path] row = { diff --git a/data_analysis/regenerate_with_localization.py b/data_analysis/regenerate_with_localization.py index 1f4750c..f09a93e 100644 --- a/data_analysis/regenerate_with_localization.py +++ b/data_analysis/regenerate_with_localization.py @@ -9,18 +9,13 @@ """ import os -import sys import json -import numpy as np import pandas as pd -from pyteomics import mzml from collections import defaultdict import matplotlib.pyplot as plt from PyPDF2 import PdfMerger -# Add GlycoSpectrumAnnotator to path -sys.path.insert(0, '/Users/longpingfu/Downloads/GlycoSpectrumAnnotator') - +from mzml_utils import open_spectra from spectrum_annotator_ddzby import ( SpectrumAnnotator, FragmentCalculator, @@ -73,23 +68,15 @@ def find_calibrated_mzml(file_name, mzml_dir): return None -def extract_spectrum_data(mzml_reader, scan_number): +def extract_spectrum_data(reader, scan_number): """Extract spectrum data for a specific scan.""" - scan_id = f"controllerType=0 controllerNumber=1 scan={scan_number}" - - try: - spectrum = mzml_reader.get_by_id(scan_id) - except KeyError: - for spec_id in mzml_reader.index.keys(): - if f"scan={scan_number}" in spec_id: - spectrum = mzml_reader.get_by_id(spec_id) - break - else: - return None + spectrum = reader.get_spectrum(int(scan_number)) + if spectrum is None: + return None return { - 'mz_array': spectrum.get('m/z array', np.array([])), - 'intensity_array': spectrum.get('intensity array', np.array([])), + 'mz_array': spectrum.mz, + 'intensity_array': spectrum.intensity, } @@ -187,7 +174,7 @@ def main(): print(f" Processing: {os.path.basename(mzml_path)} ({len(scans)} scans)...", end=" ", flush=True) try: - with mzml.MzML(mzml_path, use_index=True) as reader: + with open_spectra(mzml_path) as reader: scans_processed = 0 for idx, scan_number, row in scans: From 1dfe1c061c90f1f38e1f7e9989058c17c10be6a0 Mon Sep 17 00:00:00 2001 From: Longping Fu Date: Tue, 4 Aug 2026 20:36:18 -0400 Subject: [PATCH 2/2] Point CLAUDE.md at open_spectra, not MzMLReader The file still told the next session to use MzMLReader directly, which now contradicts both the global rule and every script in data_analysis/. Co-Authored-By: Claude Opus 5 (1M context) --- CLAUDE.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index 421f83c..f645a2f 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -182,12 +182,12 @@ PyMOL ray trace modes: ## MS/MS Spectrum Annotation (Python) -**CRITICAL: ALWAYS use `mzml_utils.MzMLReader` for spectrum access — NEVER `pyteomics.mzml` (streaming, extremely slow). Always check for cached data (pickles, filter_string_cache.csv) before opening mzML files.** +**CRITICAL: ALWAYS read spectra via `mzml_utils.open_spectra(path)` — NEVER `pyteomics.mzml` (streaming, extremely slow), and do not construct `MzMLReader` directly. `open_spectra` is a drop-in that returns a fast `SpectrumCache` when `/spectra_cache/.spectra.db` exists and an `MzMLReader` otherwise, so it is never slower and gets faster the moment a cache is built. Every active script here was converted 2026-08-04. Always check for cached data (pickles, filter_string_cache.csv) before opening mzML files.** **Use GlycoSpectrumAnnotator** (`spectrum_annotator_ddzby`, installed as editable package at `/Users/longpingfu/Downloads/GlycoSpectrumAnnotator/`) for all annotation. The local `spectrum_annotator.py` / `fragment_calculator.py` forks were **deleted 2026-08-04** — they were strict older subsets (646 vs 1490 and 1275 vs 1855 lines), missing the glycan library, N-glycan support, precursor-envelope filtering and the isotope-consistency flags. Every script behind a published spectrum already used the installed package. Import from `spectrum_annotator_ddzby` / `spectrum_annotator_ddzby.fragment_calculator`, never from a local module. Python modules for spectrum annotation: -- **mzml_utils** (`import mzml_utils`) - Indexed mzML reader (`MzMLReader`), ion search, fragment calculator, deisotoping, spectral similarity, protease digestion +- **mzml_utils** (`import mzml_utils`) - Cache-aware spectrum reader (`open_spectra`), ion search, fragment calculator, deisotoping, spectral similarity, protease digestion - **GlycoSpectrumAnnotator** (`spectrum_annotator_ddzby`) - Publication-quality annotated spectra with correct butterfly diagram, glycan labels, deisotoping, S/N filtering, charge-reduced exclusion - **OGlyco_DBA tools** (`/Users/longpingfu/Downloads/OGlyco_DBA/data_analysis/`) - `opair_utils.py`, `oglyco_validation.py`, `mass_degeneracy.py` for validation workflows - **extract_ethcd_spectra.py** - Extracts EThcD spectra from calibrated mzML files