diff --git a/alpharaw/mzml.py b/alpharaw/mzml.py index e3d5b99..77d52fb 100644 --- a/alpharaw/mzml.py +++ b/alpharaw/mzml.py @@ -8,6 +8,10 @@ ms_reader_provider, ) +SAFE_PRECURSOR_MZ = -1.0 +SAFE_ISOLATION_MZ = -1.0 +DEFAULT_ISOLATION_OFFSET = 1.5 + class MzMLReader(MSData_Base): """ @@ -112,10 +116,20 @@ def _import( peak_indices = np.empty(len(spec_indices) + 1, np.int64) peak_indices[0] = 0 peak_indices[1:] = np.cumsum(_peak_indices) + peak_mz = ( + np.concatenate(mzs_list) + if mzs_list + else np.empty(0, dtype=PEAK_MZ_DTYPE) + ) + peak_intensity = ( + np.concatenate(intens_list) + if intens_list + else np.empty(0, dtype=PEAK_INTENSITY_DTYPE) + ) ret_dict = { "peak_indices": peak_indices, - "peak_mz": np.concatenate(mzs_list), - "peak_intensity": np.concatenate(intens_list), + "peak_mz": peak_mz, + "peak_intensity": peak_intensity, "rt": np.array(rt_list), "precursor_mz": np.array(prec_mz_list), "precursor_charge": np.array(charge_list, dtype=np.int8), @@ -130,6 +144,106 @@ def _import( return ret_dict +def _parse_nce_from_filter_string(filter_string) -> float: + """Parse NCE from Thermo-like filter strings.""" + if not isinstance(filter_string, str) or not filter_string: + return np.nan + + try: + if "@hcd" in filter_string: + return float(filter_string.split("@hcd")[1].split(" ")[0]) + + if "@cid" in filter_string: + return float(filter_string.split("@cid")[1].split(" ")[0]) + + return np.nan + except (ValueError, IndexError): + return np.nan + + +def _parse_charge_state(selected_ion: dict | None) -> int: + if selected_ion is None: + return 0 + + charge_state = selected_ion.get("charge state") + + if charge_state is None: + return 0 + + try: + return int(charge_state) + except (TypeError, ValueError): + return 0 + + +def _get_first_precursor(item_dict: dict) -> dict | None: + precursor_list = item_dict.get("precursorList") + if not isinstance(precursor_list, dict): + return None + + precursors = precursor_list.get("precursor") + if not isinstance(precursors, list) or not precursors: + return None + + precursor = precursors[0] + if not isinstance(precursor, dict): + return None + + return precursor + + +def _get_first_selected_ion(precursor: dict | None) -> dict | None: + if precursor is None: + return None + + selected_ion_list = precursor.get("selectedIonList") + if not isinstance(selected_ion_list, dict): + return None + + selected_ions = selected_ion_list.get("selectedIon") + if not isinstance(selected_ions, list) or not selected_ions: + return None + + selected_ion = selected_ions[0] + if not isinstance(selected_ion, dict): + return None + + return selected_ion + + +def _get_isolation_window(precursor: dict | None) -> dict | None: + if precursor is None: + return None + + isolation_window = precursor.get("isolationWindow") + if not isinstance(isolation_window, dict): + return None + + return isolation_window + + +def _get_peak_array(item_dict: dict, key: str) -> np.ndarray: + peak_array = item_dict.get(key) + if peak_array is None: + raise KeyError(f"Missing '{key}' in mzML scan payload.") + return np.asarray(peak_array) + + +def _get_first_scan_entry(item_dict: dict) -> dict: + scan_list = item_dict.get("scanList") + if not isinstance(scan_list, dict): + raise KeyError("Missing 'scanList' in mzML scan payload.") + + scans = scan_list.get("scan") + if not isinstance(scans, list) or not scans: + raise KeyError("Missing 'scan' entries in mzML scan payload.") + + scan_entry = scans[0] + if not isinstance(scan_entry, dict): + raise TypeError("Expected first 'scan' entry to be a dict.") + return scan_entry + + def parse_mzml_entry(item_dict: dict) -> tuple: """ Parse mzml entries from pyteomics extracted items. @@ -144,71 +258,48 @@ def parse_mzml_entry(item_dict: dict) -> tuple: tuple items in tuple format. """ - rt = float(item_dict.get("scanList").get("scan")[0].get("scan start time")) - masses = item_dict.get("m/z array") - intensities = item_dict.get("intensity array") + scan_entry = _get_first_scan_entry(item_dict) + rt = float(scan_entry.get("scan start time")) + masses = _get_peak_array(item_dict, "m/z array") + intensities = _get_peak_array(item_dict, "intensity array") ms_level = item_dict.get("ms level") - prec_mz = -1.0 - isolation_lower_mz = -1.0 - isolation_upper_mz = -1.0 + prec_mz = SAFE_PRECURSOR_MZ + isolation_lower_mz = SAFE_ISOLATION_MZ + isolation_upper_mz = SAFE_ISOLATION_MZ charge = 0 nce = 0.0 if ms_level == 2: - try: - charge = int( - item_dict.get("precursorList") - .get("precursor")[0] - .get("selectedIonList") - .get("selectedIon")[0] - .get("charge state") - ) - except TypeError: - charge = 0 - try: - charge = int( - item_dict.get("precursorList") - .get("precursor")[0] - .get("selectedIonList") - .get("selectedIon")[0] - .get("charge state") - ) - except TypeError: - charge = 0 - - prec_mz = ( - item_dict.get("precursorList") - .get("precursor")[0] - .get("selectedIonList") - .get("selectedIon")[0] - .get("selected ion m/z") + precursor = _get_first_precursor(item_dict) + selected_ion = _get_first_selected_ion(precursor) + charge = _parse_charge_state(selected_ion) + + precursor_mz_value = ( + None if selected_ion is None else selected_ion.get("selected ion m/z") ) try: - iso_window = ( - item_dict.get("precursorList") - .get("precursor")[0] - .get("isolationWindow") + prec_mz = float(precursor_mz_value) + except (TypeError, ValueError): + prec_mz = SAFE_PRECURSOR_MZ + + if prec_mz != SAFE_PRECURSOR_MZ: + iso_window = _get_isolation_window(precursor) + iso_lower = None if iso_window is None else iso_window.get( + "isolation window lower offset" ) - iso_lower = float(iso_window.get("isolation window lower offset")) - iso_upper = float(iso_window.get("isolation window upper offset")) - isolation_upper_mz = prec_mz + iso_upper - isolation_lower_mz = prec_mz - iso_lower - except TypeError: - isolation_upper_mz = prec_mz + 1.5 - isolation_lower_mz = prec_mz - 1.5 - - nce = np.nan - try: - filter_string = ( - item_dict.get("scanList").get("scan")[0].get["filter string"] + iso_upper = None if iso_window is None else iso_window.get( + "isolation window upper offset" ) - if "@hcd" in filter_string: - nce = float(filter_string.split("@hcd")[1].split(" ")[0]) - elif "@cid" in filter_string: - nce = float(filter_string.split("@cid")[1].split(" ")[0]) - else: - nce = np.nan - except Exception: - nce = np.nan + try: + iso_lower = float(iso_lower) + iso_upper = float(iso_upper) + isolation_upper_mz = prec_mz + iso_upper + isolation_lower_mz = prec_mz - iso_lower + except (TypeError, ValueError): + isolation_upper_mz = prec_mz + DEFAULT_ISOLATION_OFFSET + isolation_lower_mz = prec_mz - DEFAULT_ISOLATION_OFFSET + + filter_string = scan_entry.get("filter string") + nce = _parse_nce_from_filter_string(filter_string) return ( rt, prec_mz, diff --git a/tests/unit/test_mzml_reader.py b/tests/unit/test_mzml_reader.py new file mode 100644 index 0000000..d9f3ffc --- /dev/null +++ b/tests/unit/test_mzml_reader.py @@ -0,0 +1,398 @@ +import numpy as np + +from alpharaw.mzml import MzMLReader, parse_mzml_entry +from alpharaw.ms_data_base import PEAK_INTENSITY_DTYPE, PEAK_MZ_DTYPE + + +def make_selected_ion(selected_ion_mz: float = 400.0, charge_state: int = 2) -> dict: + """Create a minimal selected-ion payload for mzML-like MS2 entries.""" + return { + "selected ion m/z": float(selected_ion_mz), + "charge state": int(charge_state), + } + + +def make_isolation_window( + target_mz: float = 400.0, + lower_offset: float = 1.5, + upper_offset: float = 1.5, +) -> dict: + """Create a minimal isolation-window payload for mzML-like MS2 entries.""" + return { + "isolation window target m/z": float(target_mz), + "isolation window lower offset": float(lower_offset), + "isolation window upper offset": float(upper_offset), + } + + +def make_ms1_entry( + rt: float = 1.0, + mz_array: np.ndarray | None = None, + intensity_array: np.ndarray | None = None, +) -> dict: + """Create a small mzML-like MS1 spectrum dictionary.""" + mz_array = ( + np.array([100.0, 200.0, 300.0], dtype=np.float32) + if mz_array is None + else mz_array + ) + intensity_array = ( + np.array([10.0, 20.0, 30.0], dtype=np.float32) + if intensity_array is None + else intensity_array + ) + + return { + "scanList": {"scan": [{"scan start time": float(rt)}]}, + "m/z array": mz_array, + "intensity array": intensity_array, + "ms level": 1, + } + + +def make_ms2_entry( + rt: float = 2.0, + precursor_mz: float = 400.0, + charge_state: int = 2, + isolation_lower_offset: float = 1.5, + isolation_upper_offset: float = 1.5, + filter_string: str | None = "FTMS + c NSI Full ms2 400.00@hcd27.00", + mz_array: np.ndarray | None = None, + intensity_array: np.ndarray | None = None, + include_isolation_window: bool = True, +) -> dict: + """Create a small mzML-like MS2 spectrum dictionary.""" + mz_array = ( + np.array([300.0, 400.0, 500.0], dtype=np.float32) + if mz_array is None + else mz_array + ) + intensity_array = ( + np.array([15.0, 25.0, 35.0], dtype=np.float32) + if intensity_array is None + else intensity_array + ) + + precursor = { + "selectedIonList": { + "selectedIon": [ + make_selected_ion( + selected_ion_mz=precursor_mz, + charge_state=charge_state, + ) + ] + } + } + if include_isolation_window: + precursor["isolationWindow"] = make_isolation_window( + target_mz=precursor_mz, + lower_offset=isolation_lower_offset, + upper_offset=isolation_upper_offset, + ) + + scan_entry = {"scan start time": float(rt)} + if filter_string is not None: + scan_entry["filter string"] = filter_string + + return { + "scanList": {"scan": [scan_entry]}, + "m/z array": mz_array, + "intensity array": intensity_array, + "ms level": 2, + "precursorList": {"precursor": [precursor]}, + } + + +class FakeMzMLReader: + """Simple iterable that mimics a pyteomics mzML reader object.""" + + def __init__(self, entries: list[dict]): + self._entries = list(entries) + self._index = 0 + self.closed = False + + def __len__(self) -> int: + return len(self._entries) + + def __iter__(self): + return self + + def __next__(self) -> dict: + if self._index >= len(self._entries): + raise StopIteration + entry = self._entries[self._index] + self._index += 1 + return entry + + def close(self) -> None: + self.closed = True + + +def parse_entry(entry: dict) -> tuple: + """Thin helper around parse_mzml_entry for upcoming behaviour tests.""" + return parse_mzml_entry(entry) + + +def test_parse_ms1_entry_without_precursor_fields(): + mz_array = np.array([111.1, 222.2, 333.3], dtype=np.float32) + intensity_array = np.array([10.0, 0.0, 55.5], dtype=np.float32) + entry = make_ms1_entry(rt=3.25, mz_array=mz_array, intensity_array=intensity_array) + + ( + rt, + precursor_mz, + isolation_lower_mz, + isolation_upper_mz, + ms_level, + nce, + precursor_charge, + parsed_mz_array, + parsed_intensity_array, + ) = parse_entry(entry) + + assert ms_level == 1 + assert precursor_mz == -1.0 + assert precursor_charge == 0 + assert isolation_lower_mz == -1.0 + assert isolation_upper_mz == -1.0 + assert rt == 3.25 + assert nce == 0.0 + np.testing.assert_array_equal(parsed_mz_array, mz_array) + np.testing.assert_array_equal(parsed_intensity_array, intensity_array) + + +def test_parse_ms2_entry_with_complete_precursor_fields(): + mz_array = np.array([350.0, 450.0, 550.0], dtype=np.float32) + intensity_array = np.array([100.0, 200.0, 300.0], dtype=np.float32) + entry = make_ms2_entry( + rt=5.5, + precursor_mz=523.27, + charge_state=3, + isolation_lower_offset=0.8, + isolation_upper_offset=1.2, + filter_string="FTMS + c NSI Full ms2 523.27@hcd29.00", + mz_array=mz_array, + intensity_array=intensity_array, + ) + + ( + rt, + precursor_mz, + isolation_lower_mz, + isolation_upper_mz, + ms_level, + nce, + precursor_charge, + parsed_mz_array, + parsed_intensity_array, + ) = parse_entry(entry) + + assert ms_level == 2 + assert rt == 5.5 + assert precursor_mz == 523.27 + assert precursor_charge == 3 + assert isolation_lower_mz == 522.47 + assert isolation_upper_mz == 524.47 + assert nce == 29.0 + np.testing.assert_array_equal(parsed_mz_array, mz_array) + np.testing.assert_array_equal(parsed_intensity_array, intensity_array) + + +def test_precursor_mz_uses_selected_ion_mz_when_isolation_target_differs(): + entry = make_ms2_entry( + precursor_mz=501.2, + isolation_lower_offset=1.0, + isolation_upper_offset=2.0, + ) + entry["precursorList"]["precursor"][0]["isolationWindow"][ + "isolation window target m/z" + ] = 500.0 + + (_, precursor_mz, _, _, _, _, _, _, _) = parse_entry(entry) + + # This locks in current behavior: precursor_mz comes from "selected ion m/z". + # Switching to isolation target m/z should be an explicit API/behavior decision. + assert np.isclose(precursor_mz, 501.2) + + +def test_isolation_bounds_use_selected_ion_mz_when_isolation_target_differs(): + entry = make_ms2_entry( + precursor_mz=501.2, + isolation_lower_offset=1.0, + isolation_upper_offset=2.0, + ) + entry["precursorList"]["precursor"][0]["isolationWindow"][ + "isolation window target m/z" + ] = 500.0 + + (_, _, isolation_lower_mz, isolation_upper_mz, _, _, _, _, _) = parse_entry(entry) + + # This locks in current behavior: offsets are applied around selected-ion m/z. + # Switching to isolation target m/z should be an explicit API/behavior decision. + assert np.isclose(isolation_lower_mz, 500.2) + assert np.isclose(isolation_upper_mz, 503.2) + + +def test_parse_ms2_entry_parses_hcd_nce_from_filter_string(): + entry = make_ms2_entry( + filter_string="FTMS + p NSI Full ms2 500.0000@hcd27.00 [100.0000-1000.0000]" + ) + + (_, _, _, _, _, nce, _, _, _) = parse_entry(entry) + + assert nce == 27.0 + + +def test_parse_ms2_entry_parses_cid_nce_from_filter_string(): + entry = make_ms2_entry( + filter_string="FTMS + p NSI Full ms2 500.0000@cid35.00 [100.0000-1000.0000]" + ) + + result = parse_mzml_entry(entry) + + assert result[5] == 35.0 + + +def test_parse_ms2_entry_missing_filter_string_returns_nan_nce(): + entry = make_ms2_entry(filter_string=None) + + (_, _, _, _, _, nce, _, _, _) = parse_entry(entry) + + assert np.isnan(nce) + + +def test_parse_ms2_entry_malformed_filter_string_returns_nan_nce(): + entry = make_ms2_entry() + entry["scanList"]["scan"][0]["filter string"] = {"not": "a-string"} + + (_, _, _, _, _, nce, _, _, _) = parse_entry(entry) + + assert np.isnan(nce) + + +def test_parse_ms2_entry_missing_charge_state_defaults_to_zero(): + entry = make_ms2_entry(charge_state=2) + entry["precursorList"]["precursor"][0]["selectedIonList"]["selectedIon"][0].pop( + "charge state" + ) + + (_, _, _, _, _, _, precursor_charge, _, _) = parse_entry(entry) + + assert precursor_charge == 0 + + +def test_parse_ms2_entry_invalid_charge_state_defaults_to_zero(): + entry = make_ms2_entry(charge_state=2) + entry["precursorList"]["precursor"][0]["selectedIonList"]["selectedIon"][0][ + "charge state" + ] = "not-a-charge" + + (_, _, _, _, _, _, precursor_charge, _, _) = parse_entry(entry) + + assert precursor_charge == 0 + + +def test_parse_ms2_entry_missing_precursor_list_uses_safe_defaults(): + entry = make_ms2_entry() + entry.pop("precursorList") + + ( + _, + precursor_mz, + isolation_lower_mz, + isolation_upper_mz, + _, + _, + precursor_charge, + _, + _, + ) = parse_entry(entry) + + assert precursor_mz == -1.0 + assert precursor_charge == 0 + assert isolation_lower_mz == -1.0 + assert isolation_upper_mz == -1.0 + + +def test_parse_ms2_entry_missing_selected_ion_uses_safe_defaults(): + entry = make_ms2_entry() + entry["precursorList"]["precursor"][0]["selectedIonList"].pop("selectedIon") + + ( + _, + precursor_mz, + isolation_lower_mz, + isolation_upper_mz, + _, + _, + precursor_charge, + _, + _, + ) = parse_entry(entry) + + assert precursor_mz == -1.0 + assert precursor_charge == 0 + assert isolation_lower_mz == -1.0 + assert isolation_upper_mz == -1.0 + + +def test_parse_ms2_entry_missing_selected_ion_mz_uses_safe_defaults(): + entry = make_ms2_entry() + entry["precursorList"]["precursor"][0]["selectedIonList"]["selectedIon"][0].pop( + "selected ion m/z" + ) + + (_, precursor_mz, isolation_lower_mz, isolation_upper_mz, _, _, _, _, _) = parse_entry( + entry + ) + + assert precursor_mz == -1.0 + assert isolation_lower_mz == -1.0 + assert isolation_upper_mz == -1.0 + assert isolation_lower_mz != -2.5 + assert isolation_upper_mz != 0.5 + + +def test_parse_empty_scan_payload(): + entry = make_ms1_entry( + mz_array=np.array([]), + intensity_array=np.array([]), + ) + + ( + _, + _, + _, + _, + _, + _, + _, + parsed_mz_array, + parsed_intensity_array, + ) = parse_entry(entry) + + assert parsed_mz_array.size == 0 + assert parsed_intensity_array.size == 0 + + parsed = MzMLReader()._import(FakeMzMLReader([entry])) + + assert parsed["peak_indices"][-1] == 0 + assert parsed["peak_mz"].size == 0 + assert parsed["peak_intensity"].size == 0 + + +def test_import_empty_reader_returns_empty_arrays(): + parsed = MzMLReader()._import(FakeMzMLReader([])) + + np.testing.assert_array_equal(parsed["peak_indices"], np.array([0], dtype=np.int64)) + assert parsed["peak_mz"].size == 0 + assert parsed["peak_mz"].dtype == np.dtype(PEAK_MZ_DTYPE) + assert parsed["peak_intensity"].size == 0 + assert parsed["peak_intensity"].dtype == np.dtype(PEAK_INTENSITY_DTYPE) + assert parsed["rt"].size == 0 + assert parsed["precursor_mz"].size == 0 + assert parsed["precursor_charge"].size == 0 + assert parsed["isolation_lower_mz"].size == 0 + assert parsed["isolation_upper_mz"].size == 0 + assert parsed["ms_level"].size == 0 + assert parsed["nce"].size == 0