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