corems.mass_spectra.input.mzml

  1from collections import defaultdict
  2from pathlib import Path
  3from typing import Optional, Union, List, Tuple
  4
  5import numpy as np
  6import pandas as pd
  7import pymzml
  8import datetime
  9
 10from corems.encapsulation.constant import Labels
 11from corems.encapsulation.factory.parameters import default_parameters
 12from corems.mass_spectra.factory.lc_class import LCMSBase, MassSpectraBase
 13from corems.mass_spectra.input.parserbase import SpectraParserInterface
 14from corems.mass_spectrum.factory.MassSpectrumClasses import (
 15    MassSpecCentroid,
 16    MassSpecProfile,
 17)
 18
 19
 20class MZMLSpectraParser(SpectraParserInterface):
 21    """A class for parsing mzml spectrometry data files into MassSpectraBase or LCMSBase objects
 22
 23    Parameters
 24    ----------
 25    file_location : str or Path
 26        The path to the RAW file to be parsed.
 27    analyzer : str, optional
 28        The type of mass analyzer used in the instrument. Default is "Unknown".
 29    instrument_label : str, optional
 30        The name of the instrument used to acquire the data. Default is "Unknown".
 31    sample_name : str, optional
 32        The name of the sample being analyzed. If not provided, the stem of the file_location path will be used.
 33
 34    Attributes
 35    ----------
 36    file_location : Path
 37        The path to the RAW file being parsed.
 38    analyzer : str
 39        The type of mass analyzer used in the instrument.
 40    instrument_label : str
 41        The name of the instrument used to acquire the data.
 42    sample_name : str
 43        The name of the sample being analyzed.
 44
 45    Methods
 46    -------
 47    * load().
 48        Load mzML file using pymzml.run.Reader and return the data as a numpy array.
 49    * run(spectra=True).
 50        Parses the mzml file and returns a dictionary of mass spectra dataframes and a scan metadata dataframe.
 51    * get_mass_spectrum_from_scan(scan_number, polarity, auto_process=True)
 52        Parses the mzml file and returns a MassSpecBase object from a single scan.
 53    * get_mass_spectra_obj().
 54        Parses the mzml file and instantiates a MassSpectraBase object.
 55    * get_lcms_obj().
 56        Parses the mzml file and instantiates an LCMSBase object.
 57    * get_instrument_info().
 58        Return instrument information from the mzML file.
 59    * get_creation_time().
 60        Return the creation time of the mzML file as a datetime object.
 61
 62    Inherits from ThermoBaseClass and SpectraParserInterface
 63    """
 64
 65    def __init__(
 66        self,
 67        file_location,
 68        analyzer="Unknown",
 69        instrument_label="Unknown",
 70        sample_name=None,
 71    ):
 72        # implementation details
 73        if isinstance(file_location, str):
 74            # if obj is a string it defaults to create a Path obj, pass the S3Path if needed
 75            file_location = Path(file_location)
 76        if not file_location.exists():
 77            raise FileExistsError("File does not exist: " + str(file_location))
 78        self.file_location = file_location
 79        self.analyzer = analyzer
 80        self.instrument_label = instrument_label
 81
 82        if sample_name:
 83            self.sample_name = sample_name
 84        else:
 85            self.sample_name = file_location.stem
 86
 87    def load(self):
 88        """
 89        Load mzML file using pymzml.run.Reader and return the data as a numpy array.
 90
 91        Returns
 92        -------
 93        numpy.ndarray
 94            The mass spectra data as a numpy array.
 95        """
 96        data = pymzml.run.Reader(self.file_location)
 97        return data
 98
 99    def get_scan_df(self, data=None, time_range: Optional[Union[Tuple[float, float], List[Tuple[float, float]]]] = None):
100        """
101        Return scan data as a pandas DataFrame.
102
103        Parameters
104        ----------
105        data : pymzml.run.Reader, optional
106            The mass spectra data. If None, will load the data.
107        time_range : tuple or list of tuples, optional
108            Retention time range(s) to filter scans. Can be:
109            - Single range: (start_time, end_time) in minutes
110            - Multiple ranges: [(start1, end1), (start2, end2), ...] in minutes
111            If None, returns all scans.
112
113        Returns
114        -------
115        pandas.DataFrame
116            A pandas DataFrame containing metadata for each scan, including scan number, MS level, polarity, and scan time.
117        """
118        if data is None:
119            data = self.load()
120        # Scan dict
121        # instatinate scan dict, with empty lists of size of scans
122        n_scans = data.get_spectrum_count()
123        scan_dict = {
124            "scan": np.empty(n_scans, dtype=np.int32),
125            "scan_time": np.empty(n_scans, dtype=np.float32),
126            "ms_level": [None] * n_scans,
127            "polarity": [None] * n_scans,
128            "precursor_mz": [None] * n_scans,
129            "scan_text": [None] * n_scans,
130            "scan_window_lower": np.empty(n_scans, dtype=np.float32),
131            "scan_window_upper": np.empty(n_scans, dtype=np.float32),
132            "scan_precision": [None] * n_scans,
133            "tic": np.empty(n_scans, dtype=np.float32),
134            "ms_format": [None] * n_scans,
135        }
136
137        # First pass: loop through scans to get scan info
138        for i, spec in enumerate(data):
139            scan_dict["scan"][i] = spec.ID
140            scan_dict["ms_level"][i] = spec.ms_level
141            scan_dict["scan_precision"][i] = spec._measured_precision
142            scan_dict["tic"][i] = spec.TIC
143            if spec.selected_precursors:
144                scan_dict["precursor_mz"][i] = spec.selected_precursors[0].get(
145                    "mz", None
146                )
147            if spec["negative scan"] is not None:
148                scan_dict["polarity"][i] = "negative"
149            if spec["positive scan"] is not None:
150                scan_dict["polarity"][i] = "positive"
151            if spec["negative scan"] is not None and spec["positive scan"] is not None:
152                raise ValueError(
153                    "Error: scan {0} has both negative and positive polarity".format(
154                        spec.ID
155                    )
156                )
157
158            scan_dict["scan_time"][i] = spec.get("MS:1000016")
159            scan_dict["scan_text"][i] = spec.get("MS:1000512")
160            scan_dict["scan_window_lower"][i] = spec.get("MS:1000501")
161            scan_dict["scan_window_upper"][i] = spec.get("MS:1000500")
162            if spec.get("MS:1000128"):
163                scan_dict["ms_format"][i] = "profile"
164            elif spec.get("MS:1000127"):
165                scan_dict["ms_format"][i] = "centroid"
166            else:
167                scan_dict["ms_format"][i] = None
168
169        scan_df = pd.DataFrame(scan_dict)
170
171        # Remove any non-mass spectra scans (e.g., MS level 0 or None)
172        scan_df = scan_df[scan_df.ms_level.notnull() & (scan_df.ms_level > 0)].reset_index(drop=True)
173        
174        # Apply time range filtering if specified
175        if time_range is not None:
176            time_ranges = self._normalize_time_range(time_range)
177            # Create a mask for scans within any of the time ranges
178            mask = np.zeros(len(scan_df), dtype=bool)
179            for start_time, end_time in time_ranges:
180                mask |= (scan_df["scan_time"] >= start_time) & (scan_df["scan_time"] <= end_time)
181            scan_df = scan_df[mask].reset_index(drop=True)
182
183        return scan_df
184
185    def get_ms_raw(self, spectra, scan_df, data=None, time_range: Optional[Union[Tuple[float, float], List[Tuple[float, float]]]] = None):
186        """Return a dictionary of mass spectra data as a pandas DataFrame.
187
188        Parameters
189        ----------
190        spectra : str
191            Which mass spectra data to include in the output.
192            Options: None, "ms1", "ms2", "all".
193        scan_df : pandas.DataFrame
194            Scan dataframe. Output from get_scan_df().
195        data : pymzml.run.Reader, optional
196            The mass spectra data. If None, will load the data.
197        time_range : tuple or list of tuples, optional
198            Retention time range(s) to filter scans. Can be:
199            - Single range: (start_time, end_time) in minutes
200            - Multiple ranges: [(start1, end1), (start2, end2), ...] in minutes
201            If None, returns all scans. Note: filtering is typically done at scan_df level.
202
203        Returns
204        -------
205        dict
206            A dictionary containing the mass spectra data as pandas DataFrames, with keys corresponding to the MS level.
207
208        """
209        if data is None:
210            data = self.load()
211        if spectra == "all":
212            scan_df_forspec = scan_df
213        elif spectra == "ms1":
214            scan_df_forspec = scan_df[scan_df.ms_level == 1]
215        elif spectra == "ms2":
216            scan_df_forspec = scan_df[scan_df.ms_level == 2]
217        else:
218            raise ValueError("spectra must be 'all', 'ms1', or 'ms2'")
219
220        # Result container
221        res = {}
222
223        # Row count container
224        counter = {}
225
226        # Column name container
227        cols = {}
228
229        # set at float32
230        dtype = np.float32
231
232        # First pass: get nrows
233        N = defaultdict(lambda: 0)
234        for i, spec in enumerate(data):
235            if spec.ID in scan_df_forspec.scan.values:
236                # Get ms level
237                level = "ms{}".format(spec.ms_level)
238
239                # Number of rows
240                N[level] += spec.mz.shape[0]
241
242        # Second pass: parse
243        for i, spec in enumerate(data):
244            if spec.ID in scan_df_forspec.scan.values:
245                # Number of rows
246                n = spec.mz.shape[0]
247
248                # No measurements
249                if n == 0:
250                    continue
251
252                # Dimension check
253                if len(spec.mz) != len(spec.i):
254                    # raise an error if the mz and intensity arrays are not the same length
255                    raise ValueError("m/z and intensity array dimension mismatch")
256
257                # Scan/frame info
258                id_dict = spec.id_dict
259
260                # Get ms level
261                level = "ms{}".format(spec.ms_level)
262
263                # Columns
264                cols[level] = list(id_dict.keys()) + ["mz", "intensity"]
265                m = len(cols[level])
266
267                # Subarray init
268                arr = np.empty((n, m), dtype=dtype)
269                inx = 0
270
271                # Populate scan/frame info
272                for k, v in id_dict.items():
273                    arr[:, inx] = v
274                    inx += 1
275
276                # Populate m/z
277                arr[:, inx] = spec.mz
278                inx += 1
279
280                # Populate intensity
281                arr[:, inx] = spec.i
282                inx += 1
283
284                # Initialize output container
285                if level not in res:
286                    res[level] = np.empty((N[level], m), dtype=dtype)
287                    counter[level] = 0
288
289                # Insert subarray
290                res[level][counter[level] : counter[level] + n, :] = arr
291                counter[level] += n
292
293        # Construct ms1 and ms2 mz dataframes
294        for level in res.keys():
295            res[level] = pd.DataFrame(res[level], columns=cols[level]).drop(
296                columns=["controllerType", "controllerNumber"],
297            )
298
299        return res
300
301    def run(self, spectra="all", scan_df=None, time_range: Optional[Union[Tuple[float, float], List[Tuple[float, float]]]] = None):
302        """Parse the mzML file and return a dictionary of spectra dataframes and a scan metadata dataframe.
303
304        Parameters
305        ----------
306        spectra : str, optional
307            Which mass spectra data to include in the output. Default is "all".
308            Other options: None, "ms1", "ms2".
309        scan_df : pandas.DataFrame, optional
310            Scan dataframe.  If not provided, the scan dataframe is created from the mzML file.
311        time_range : tuple or list of tuples, optional
312            Retention time range(s) to load. Can be:
313            - Single range: (start_time, end_time) in minutes
314            - Multiple ranges: [(start1, end1), (start2, end2), ...] in minutes
315            If None, loads all scans.
316
317        Returns
318        -------
319        tuple
320            A tuple containing two elements:
321            - A dictionary containing the mass spectra data as numpy arrays, with keys corresponding to the MS level.
322            - A pandas DataFrame containing metadata for each scan, including scan number, MS level, polarity, and scan time.
323        """
324
325        # Open file
326        data = self.load()
327
328        if scan_df is None:
329            scan_df = self.get_scan_df(data, time_range=time_range)
330
331        if spectra != "none":
332            res = self.get_ms_raw(spectra, scan_df, data)
333
334        else:
335            res = None
336
337        return res, scan_df
338
339    def get_mass_spectrum_from_scan(
340        self, scan_number, spectrum_mode, auto_process=True
341    ):
342        """Instatiate a mass spectrum object from the mzML file.
343
344        Parameters
345        ----------
346        scan_number : int
347            The scan number to be parsed.
348        spectrum_mode : str
349            The type of spectrum to instantiate.  Must be'profile' or 'centroid'.
350        polarity : int
351            The polarity of the scan.  Must be -1 or 1.
352        auto_process : bool, optional
353            If True, process the mass spectrum. Default is True.
354
355        Returns
356        -------
357        MassSpecProfile | MassSpecCentroid
358            The MassSpecProfile or MassSpecCentroid object containing the parsed mass spectrum.
359        """
360        # Use the batch function and return the first result
361        result_list = self.get_mass_spectra_from_scan_list(
362            [scan_number], spectrum_mode, auto_process
363        )
364        return result_list[0] if result_list else None
365
366    def get_mass_spectra_from_scan_list(
367        self, scan_list, spectrum_mode, auto_process=True
368    ):
369        """Instatiate mass spectrum objects from the mzML file.
370
371        Parameters
372        ----------
373        scan_list : list of int
374            The scan numbers to be parsed.
375        spectrum_mode : str
376            The type of spectrum to instantiate.  Must be'profile' or 'centroid'.
377        auto_process : bool, optional
378            If True, process the mass spectrum. Default is True.
379
380        Returns
381        -------
382        list of MassSpecProfile | MassSpecCentroid
383            List of MassSpecProfile or MassSpecCentroid objects containing the parsed mass spectra.
384        """
385
386        def set_metadata(
387            scan_number: int,
388            polarity: int,
389            file_location: str,
390            label=Labels.thermo_profile,
391        ):
392            """
393            Set the output parameters for creating a MassSpecProfile or MassSpecCentroid object.
394
395            Parameters
396            ----------
397            scan_number : int
398                The scan number.
399            polarity : int
400                The polarity of the data.
401            file_location : str
402                The file location.
403            label : str, optional
404                The label for the mass spectrum. Default is Labels.thermo_profile.
405
406            Returns
407            -------
408            dict
409                The output parameters ready for creating a MassSpecProfile or MassSpecCentroid object.
410            """
411            d_params = default_parameters(file_location)
412            d_params["label"] = label
413            d_params["polarity"] = polarity
414            d_params["filename_path"] = file_location
415            d_params["scan_number"] = scan_number
416
417            return d_params
418
419        # Open file
420        data = self.load()
421
422        mass_spectrum_objects = []
423        scan_set = set(scan_list)
424        # Materialize peak data during iteration. Do not retain live pymzml
425        # spectrum objects, and do not overwrite the first valid MS spectrum for
426        # a scan ID.
427        #
428        # Thermo multi-controller mzML files often reuse scan numbers across
429        # controllers (e.g. controllerType=0 MS data and controllerType=3
430        # auxiliary traces). pymzml exposes only the integer scan number as
431        # spec.ID, so sequential iteration can yield the same ID twice. Keeping
432        # the last yield (or random-access data[scan]) can replace real MS
433        # spectra with non-MS controller spectra and then fail centroid/profile
434        # checks with "spectrum is not centroided".
435        #
436        # Direct random-access via data[scan_number] also uses pymzml's
437        # byte-offset index, which is unreliable on Windows (CRLF vs LF).
438        collected = {}
439
440        for spec in data:
441            scan_id = spec.ID
442            if scan_id not in scan_set or scan_id in collected:
443                continue
444
445            # Skip non-MS / auxiliary-controller spectra that share scan numbers
446            ms_level = spec.ms_level
447            if ms_level is None or ms_level <= 0:
448                continue
449
450            if spec["negative scan"] is not None:
451                polarity = -1
452            elif spec["positive scan"] is not None:
453                polarity = 1
454            else:
455                polarity = None
456
457            mz = np.asarray(spec.mz)
458            abundance = np.asarray(spec.i)
459            collected[scan_id] = {
460                "mz": mz,
461                "abundance": abundance,
462                "polarity": polarity,
463                "is_profile": bool(spec.get("MS:1000128")),
464                "is_centroid": bool(spec.get("MS:1000127")),
465            }
466
467        for scan_number in scan_list:
468            entry = collected.get(scan_number)
469            if entry is None:
470                raise ValueError(
471                    "Scan number %d not found in mzML file" % scan_number
472                )
473
474            polarity = entry["polarity"]
475            mz = entry["mz"]
476            abundance = entry["abundance"]
477
478            # Get mass spectrum
479            if spectrum_mode == "profile":
480                # Check if profile
481                if not entry["is_profile"]:
482                    raise ValueError("spectrum is not profile")
483                data_dict = {
484                    Labels.mz: mz,
485                    Labels.abundance: abundance,
486                }
487                d_params = set_metadata(
488                    scan_number,
489                    polarity,
490                    self.file_location,
491                    label=Labels.simulated_profile,
492                )
493                mass_spectrum_obj = MassSpecProfile(
494                    data_dict, d_params, auto_process=auto_process
495                )
496            elif spectrum_mode == "centroid":
497                # Check if centroided
498                if not entry["is_centroid"]:
499                    raise ValueError("spectrum is not centroided")
500                data_dict = {
501                    Labels.mz: mz,
502                    Labels.abundance: abundance,
503                    Labels.rp: [np.nan] * len(mz),
504                    Labels.s2n: [np.nan] * len(abundance),
505                }
506                d_params = set_metadata(
507                    scan_number, polarity, self.file_location, label=Labels.corems_centroid
508                )
509                mass_spectrum_obj = MassSpecCentroid(
510                    data_dict, d_params, auto_process=auto_process
511                )
512            else:
513                raise ValueError(
514                    "spectrum_mode must be 'profile' or 'centroid', got %r"
515                    % spectrum_mode
516                )
517
518            mass_spectrum_objects.append(mass_spectrum_obj)
519
520        return mass_spectrum_objects
521
522    def get_mass_spectra_obj(self, time_range: Optional[Union[Tuple[float, float], List[Tuple[float, float]]]] = None):
523        """Instatiate a MassSpectraBase object from the mzML file.
524
525        Parameters
526        ----------
527        time_range : tuple or list of tuples, optional
528            Retention time range(s) to load. Can be:
529            - Single range: (start_time, end_time) in minutes
530            - Multiple ranges: [(start1, end1), (start2, end2), ...] in minutes
531            If None, loads all scans.
532
533        Returns
534        -------
535        MassSpectraBase
536            The MassSpectra object containing the parsed mass spectra.
537            The object is instatiated with the mzML file, analyzer, instrument, sample name, and scan dataframe.
538        """
539        _, scan_df = self.run(spectra=False, time_range=time_range)
540        mass_spectra_obj = MassSpectraBase(
541            self.file_location,
542            self.analyzer,
543            self.instrument_label,
544            self.sample_name,
545            self,
546        )
547        scan_df = scan_df.set_index("scan", drop=False)
548        mass_spectra_obj.scan_df = scan_df
549
550        return mass_spectra_obj
551
552    def get_lcms_obj(self, spectra="all", time_range: Optional[Union[Tuple[float, float], List[Tuple[float, float]]]] = None):
553        """Instatiates a LCMSBase object from the mzML file.
554
555        Parameters
556        ----------
557        spectra : str, optional
558            Which mass spectra data to include in the output. Default is all.  Other options: none, ms1, ms2.
559        time_range : tuple or list of tuples, optional
560            Retention time range(s) to load. Can be:
561            - Single range: (start_time, end_time) in minutes
562            - Multiple ranges: [(start1, end1), (start2, end2), ...] in minutes
563            If None, loads all scans. Useful for targeted workflows to improve performance.
564
565        Returns
566        -------
567        LCMSBase
568            LCMS object containing mass spectra data.
569            The object is instatiated with the mzML file, analyzer, instrument, sample name, scan dataframe,
570            and mz dataframe(s), as well as lists of scan numbers, retention times, and TICs.
571        """
572        _, scan_df = self.run(spectra="none", time_range=time_range)  # first run it to just get scan info
573        if spectra != "none":
574            res, scan_df = self.run(
575                scan_df=scan_df, spectra=spectra, time_range=time_range
576            )  # second run to parse data
577        lcms_obj = LCMSBase(
578            self.file_location,
579            self.analyzer,
580            self.instrument_label,
581            self.sample_name,
582            self,
583        )
584        if spectra != "none":
585            for key in res:
586                key_int = int(key.replace("ms", ""))
587                res[key] = res[key][res[key].intensity > 0]
588                res[key] = res[key].sort_values(by=["scan", "mz"]).reset_index(drop=True)
589                lcms_obj._ms_unprocessed[key_int] = res[key]
590        lcms_obj.scan_df = scan_df.set_index("scan", drop=False)
591        # Check if polarity is mixed
592        if len(set(scan_df.polarity)) > 1:
593            raise ValueError("Mixed polarities detected in scan data")
594        lcms_obj.polarity = scan_df.polarity[0]
595        lcms_obj._scans_number_list = list(scan_df.scan)
596        lcms_obj._retention_time_list = list(scan_df.scan_time)
597        lcms_obj._tic_list = list(scan_df.tic)
598
599        return lcms_obj
600
601    def get_scans_in_time_range(
602        self, 
603        time_range: Union[Tuple[float, float], List[Tuple[float, float]]],
604        ms_level: Optional[int] = None
605    ) -> List[int]:
606        """
607        Return scan numbers within specified retention time range(s).
608        
609        This method provides efficient filtering of scans by retention time,
610        which is particularly useful for targeted workflows where only specific
611        time windows are of interest.
612        
613        Parameters
614        ----------
615        time_range : tuple or list of tuples
616            Retention time range(s) in minutes. Can be:
617            - Single range: (start_time, end_time)
618            - Multiple ranges: [(start1, end1), (start2, end2), ...]
619        ms_level : int, optional
620            If specified, only return scans of this MS level (e.g., 1 for MS1, 2 for MS2).
621            If None, returns scans of all MS levels.
622        
623        Returns
624        -------
625        list of int
626            List of scan numbers within the specified time range(s) and MS level.
627        
628        Examples
629        --------
630        Get MS1 scans between 1.0 and 2.0 minutes:
631        
632        >>> scans = parser.get_scans_in_time_range((1.0, 2.0), ms_level=1)
633        
634        Get scans in multiple time windows:
635        
636        >>> scans = parser.get_scans_in_time_range([(0.5, 1.5), (3.0, 4.0)])
637        """
638        # Get scan dataframe filtered by time range
639        scan_df = self.get_scan_df(time_range=time_range)
640        
641        # Further filter by MS level if specified
642        if ms_level is not None:
643            scan_df = scan_df[scan_df.ms_level == ms_level]
644        
645        # Return list of scan numbers
646        return scan_df.scan.tolist()
647
648    def get_instrument_info(self):
649        """
650        Return instrument information.
651
652        Returns
653        -------
654        dict
655            A dictionary with the keys 'model' and 'serial_number'.
656        """
657        # Load the pymzml data
658        data = self.load()
659        instrument_info = data.info.get('referenceable_param_group_list_element')[0]
660        cv_params = instrument_info.findall('{http://psi.hupo.org/ms/mzml}cvParam')
661
662        # Extract details from each cvParam
663        params = []
664        for param in cv_params:
665            accession = param.get('accession')  # Get 'accession' attribute
666            name = param.get('name')           # Get 'name' attribute
667            value = param.get('value')         # Get 'value' attribute
668            params.append({
669                'accession': accession,
670                'name': name,
671                'value': value
672            })
673
674        # Loop through params and try to find the relevant information
675        instrument_dict = {
676            'model': 'Unknown',
677            'serial_number': 'Unknown'
678        }
679
680        # Assuming there are only two paramters here - one is for the serial number (agnostic to the model) and the other is for the model
681        # If there are more than two, we raise an error
682        if len(params) < 2:
683            raise ValueError("Not enough parameters found in the instrument info, cannot parse.")
684        if len(params) > 2:
685            raise ValueError("Too many parameters found in the instrument info, cannot parse.")
686        for param in params:
687            if param['accession'] == 'MS:1000529':
688                instrument_dict['serial_number'] = param['value']
689            else:
690                instrument_dict['model'] = data.OT[param['accession']]     
691
692        return instrument_dict
693
694    def get_creation_time(self) -> datetime.datetime:
695        """
696        Return the creation time of the mzML file.
697        """
698        data = self.load()
699        write_time = data.info.get('start_time')
700        if write_time:
701            # Convert the write time to a datetime object
702            return datetime.datetime.strptime(write_time, "%Y-%m-%dT%H:%M:%SZ")
703        else:
704            raise ValueError("Creation time is not available in the mzML file. "
705                           "Please ensure the file contains the 'start_time' information.")