Route all spectrum reading through open_spectra; unpin AlphaFold models

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 `<mzml_dir>/spectra_cache/<stem>.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

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

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

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)

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)

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

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

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

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):

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:

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():

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