corems.mass_spectrum.factory.MassSpectrumClasses

   1from pathlib import Path
   2
   3import numpy as np
   4from lmfit.models import GaussianModel
   5
   6# from matplotlib import rcParamsDefault, rcParams
   7from numpy import argmax, array, asarray, float64, histogram, searchsorted, where
   8
   9try:
  10    from numpy import trapezoid
  11except ImportError:  # numpy < 2.0
  12    from numpy import trapz as trapezoid
  13from pandas import DataFrame
  14
  15from corems.encapsulation.constant import Labels
  16from corems.encapsulation.factory.parameters import MSParameters
  17from corems.encapsulation.input.parameter_from_json import (
  18    load_and_set_parameters_ms,
  19    load_and_set_toml_parameters_ms,
  20)
  21from corems.mass_spectrum.calc.KendrickGroup import KendrickGrouping
  22from corems.mass_spectrum.calc.MassSpectrumCalc import MassSpecCalc
  23from corems.mass_spectrum.calc.MeanResolvingPowerFilter import MeanResolvingPowerFilter
  24from corems.ms_peak.factory.MSPeakClasses import ICRMassPeak as MSPeak
  25
  26__author__ = "Yuri E. Corilo"
  27__date__ = "Jun 12, 2019"
  28
  29
  30def overrides(interface_class):
  31    """Checks if the method overrides a method from an interface class."""
  32
  33    def overrider(method):
  34        assert method.__name__ in dir(interface_class)
  35        return method
  36
  37    return overrider
  38
  39
  40def _nearest_index_either_order(profile, value):
  41    """Index into ``profile`` nearest ``value``; ``profile`` may run either
  42    ascending or descending (e.g. a frequency-derived spectrum's m/z = A/f
  43    runs high to low).
  44
  45    Clipped to ``[1, len(profile) - 2]`` so the result is always safe to
  46    pass to ``find_minima()``, which looks at the neighbor on each side.
  47    """
  48    n = len(profile)
  49    if profile[0] <= profile[-1]:
  50        idx = int(searchsorted(profile, value))
  51    else:
  52        idx = n - 1 - int(searchsorted(profile[::-1], value))
  53    return min(max(idx, 1), n - 2)
  54
  55
  56class MassSpecBase(MassSpecCalc, KendrickGrouping):
  57    """A mass spectrum base class, stores the profile data and instrument settings.
  58
  59    Iteration over a list of MSPeaks classes stored at the _mspeaks attributes.
  60    _mspeaks is populated under the hood by calling process_mass_spec method.
  61    Iteration is null if _mspeaks is empty.
  62
  63    Parameters
  64    ----------
  65    mz_exp : array_like
  66        The m/z values of the mass spectrum.
  67    abundance : array_like
  68        The abundance values of the mass spectrum.
  69    d_params : dict
  70        A dictionary of parameters for the mass spectrum.
  71    **kwargs
  72        Additional keyword arguments.
  73
  74    Attributes
  75    ----------
  76
  77    mspeaks : list
  78        A list of mass peaks.
  79    is_calibrated : bool
  80        Whether the mass spectrum is calibrated.
  81    is_centroid : bool
  82        Whether the mass spectrum is centroided.
  83    has_frequency : bool
  84        Whether the mass spectrum has a frequency domain.
  85    calibration_order : None or int
  86        The order of the mass spectrum's calibration.
  87    calibration_points : None or ndarray
  88        The calibration points of the mass spectrum.
  89    calibration_ref_mzs: None or ndarray
  90        The reference m/z values of the mass spectrum's calibration.
  91    calibration_meas_mzs : None or ndarray
  92        The measured m/z values of the mass spectrum's calibration.
  93    calibration_RMS : None or float
  94        The root mean square of the mass spectrum's calibration.
  95    calibration_segment : None or CalibrationSegment
  96        The calibration segment of the mass spectrum.
  97    _abundance : ndarray
  98        The abundance values of the mass spectrum.
  99    _mz_exp : ndarray
 100        The m/z values of the mass spectrum.
 101    _mspeaks : list
 102        A list of mass peaks.
 103    _dict_nominal_masses_indexes : dict
 104        A dictionary of nominal masses and their indexes.
 105    _baseline_noise : float
 106        The baseline noise of the mass spectrum.
 107    _baseline_noise_std : float
 108        The standard deviation of the baseline noise of the mass spectrum.
 109    _dynamic_range : float or None
 110        The dynamic range of the mass spectrum.
 111    _transient_settings : None or TransientSettings
 112        The transient settings of the mass spectrum.
 113    _frequency_domain : None or FrequencyDomain
 114        The frequency domain of the mass spectrum.
 115    _mz_cal_profile : None or MzCalibrationProfile
 116        The m/z calibration profile of the mass spectrum.
 117
 118    Methods
 119    -------
 120    * process_mass_spec(). Main function to process the mass spectrum,
 121    including calculating the noise threshold, peak picking, and resetting the MSpeak indexes.
 122
 123    See also: MassSpecCentroid(), MassSpecfromFreq(), MassSpecProfile()
 124    """
 125
 126    def __init__(self, mz_exp, abundance, d_params, **kwargs):
 127        self._abundance = array(abundance, dtype=float64)
 128        self._mz_exp = array(mz_exp, dtype=float64)
 129
 130        # objects created after process_mass_spec() function
 131        self._mspeaks = list()
 132        self.mspeaks = list()
 133        self._dict_nominal_masses_indexes = dict()
 134        self._baseline_noise = 0.001
 135        self._baseline_noise_std = 0.001
 136        self._dynamic_range = None
 137        # set to None: initialization occurs inside subclass MassSpecfromFreq
 138        self._transient_settings = None
 139        self._frequency_domain = None
 140        self._mz_cal_profile = None
 141        self.is_calibrated = False
 142
 143        self._set_parameters_objects(d_params)
 144        self._init_settings()
 145
 146        self.is_centroid = False
 147        self.has_frequency = False
 148
 149        self.calibration_order = None
 150        self.calibration_points = None
 151        self.calibration_ref_mzs = None
 152        self.calibration_meas_mzs = None
 153        self.calibration_RMS = None
 154        self.calibration_segment = None
 155        self.calibration_raw_error_median = None
 156        self.calibration_raw_error_stdev = None
 157
 158    def _init_settings(self):
 159        """Initializes the settings for the mass spectrum."""
 160        self._parameters = MSParameters()
 161
 162    def __len__(self):
 163        return len(self.mspeaks)
 164
 165    def __getitem__(self, position) -> MSPeak:
 166        return self.mspeaks[position]
 167
 168    def set_indexes(self, list_indexes):
 169        """Set the mass spectrum to iterate over only the selected MSpeaks indexes.
 170
 171        Parameters
 172        ----------
 173        list_indexes : list of int
 174            A list of integers representing the indexes of the MSpeaks to iterate over.
 175
 176        """
 177        self.mspeaks = [self._mspeaks[i] for i in list_indexes]
 178
 179        for i, mspeak in enumerate(self.mspeaks):
 180            mspeak.index = i
 181
 182        self._set_nominal_masses_start_final_indexes()
 183
 184    def reset_indexes(self):
 185        """Reset the mass spectrum to iterate over all MSpeaks objects.
 186
 187        This method resets the mass spectrum to its original state, allowing iteration over all MSpeaks objects.
 188        It also sets the index of each MSpeak object to its corresponding position in the mass spectrum.
 189
 190        """
 191        self.mspeaks = self._mspeaks
 192
 193        for i, mspeak in enumerate(self.mspeaks):
 194            mspeak.index = i
 195
 196        self._set_nominal_masses_start_final_indexes()
 197
 198    def add_mspeak(
 199        self,
 200        polarity,
 201        mz_exp,
 202        abundance,
 203        resolving_power,
 204        signal_to_noise,
 205        massspec_indexes,
 206        exp_freq=None,
 207        ms_parent=None,
 208        **kwargs,
 209    ):
 210        """Add a new MSPeak object to the MassSpectrum object.
 211
 212        Parameters
 213        ----------
 214        polarity : int
 215            Acquisition polarity of the MSPeak (``+1`` or ``-1``).
 216            Do **not** pass ``ion_charge``; peaks store polarity only.
 217            Assignment charge lives on ``MolecularFormula.ion_charge``.
 218        mz_exp : float
 219            The experimental m/z value of the MSPeak.
 220        abundance : float
 221            The abundance of the MSPeak.
 222        resolving_power : float
 223            The resolving power of the MSPeak.
 224        signal_to_noise : float
 225            The signal-to-noise ratio of the MSPeak.
 226        massspec_indexes : list
 227            A list of indexes of the MSPeak in the MassSpectrum object.
 228        exp_freq : float, optional
 229            The experimental frequency of the MSPeak. Defaults to None.
 230        ms_parent : MSParent, optional
 231            The MSParent object associated with the MSPeak. Defaults to None.
 232
 233        Raises
 234        ------
 235        TypeError
 236            If ``ion_charge`` is passed (removed; use ``polarity``), or if any
 237            other unexpected keyword argument is given.
 238        """
 239        if "ion_charge" in kwargs:
 240            raise TypeError(
 241                "add_mspeak() no longer accepts ion_charge=...; pass polarity "
 242                "(±1) instead. Peaks store acquisition polarity only. Assignment "
 243                "charge belongs on MolecularFormula.ion_charge after formula matching."
 244            )
 245        if kwargs:
 246            unexpected = ", ".join(sorted(kwargs))
 247            raise TypeError(
 248                f"add_mspeak() got unexpected keyword argument(s): {unexpected}"
 249            )
 250
 251        mspeak = MSPeak(
 252            polarity,
 253            mz_exp,
 254            abundance,
 255            resolving_power,
 256            signal_to_noise,
 257            massspec_indexes,
 258            len(self._mspeaks),
 259            exp_freq=exp_freq,
 260            ms_parent=ms_parent,
 261        )
 262
 263        self._mspeaks.append(mspeak)
 264
 265    def _set_parameters_objects(self, d_params):
 266        """Set the parameters of the MassSpectrum object.
 267
 268        Parameters
 269        ----------
 270        d_params : dict
 271            A dictionary containing the parameters to set.
 272
 273        Notes
 274        -----
 275        This method sets the following parameters of the MassSpectrum object:
 276        - _calibration_terms
 277        - label
 278        - analyzer
 279        - acquisition_time
 280        - instrument_label
 281        - polarity
 282        - scan_number
 283        - retention_time
 284        - mobility_rt
 285        - mobility_scan
 286        - _filename
 287        - _dir_location
 288        - _baseline_noise
 289        - _baseline_noise_std
 290        - sample_name
 291        """
 292        self._calibration_terms = (
 293            d_params.get("Aterm"),
 294            d_params.get("Bterm"),
 295            d_params.get("Cterm"),
 296        )
 297
 298        self.label = d_params.get(Labels.label)
 299
 300        self.analyzer = d_params.get("analyzer")
 301
 302        self.acquisition_time = d_params.get("acquisition_time")
 303
 304        self.instrument_label = d_params.get("instrument_label")
 305
 306        self.polarity = int(d_params.get("polarity"))
 307
 308        self.scan_number = d_params.get("scan_number")
 309
 310        self.retention_time = d_params.get("rt")
 311
 312        self.mobility_rt = d_params.get("mobility_rt")
 313
 314        self.mobility_scan = d_params.get("mobility_scan")
 315
 316        self._filename = d_params.get("filename_path")
 317
 318        self._dir_location = d_params.get("dir_location")
 319
 320        self._baseline_noise = d_params.get("baseline_noise")
 321
 322        self._baseline_noise_std = d_params.get("baseline_noise_std")
 323
 324        if d_params.get("sample_name") != "Unknown":
 325            self.sample_name = d_params.get("sample_name")
 326            if not self.sample_name:
 327                self.sample_name = self.filename.stem
 328        else:
 329            self.sample_name = self.filename.stem
 330
 331    def reset_cal_therms(self, Aterm, Bterm, C, fas=0):
 332        """Reset calibration terms and recalculate the mass-to-charge ratio and abundance.
 333
 334        Parameters
 335        ----------
 336        Aterm : float
 337            The A-term calibration coefficient.
 338        Bterm : float
 339            The B-term calibration coefficient.
 340        C : float
 341            The C-term calibration coefficient.
 342        fas : float, optional
 343            The frequency amplitude scaling factor. Default is 0.
 344        """
 345        self._calibration_terms = (Aterm, Bterm, C)
 346
 347        self._mz_exp = self._f_to_mz()
 348        self._abundance = self._abundance
 349        self.find_peaks()
 350        self.reset_indexes()
 351
 352    def clear_molecular_formulas(self):
 353        """Clear the molecular formulas for all mspeaks in the MassSpectrum.
 354
 355        Returns
 356        -------
 357        numpy.ndarray
 358            An array of the cleared molecular formulas for each mspeak in the MassSpectrum.
 359        """
 360        self.check_mspeaks()
 361        return array([mspeak.clear_molecular_formulas() for mspeak in self.mspeaks])
 362
 363    def process_mass_spec(self, keep_profile=True):
 364        """Process the mass spectrum.
 365
 366        Parameters
 367        ----------
 368        keep_profile : bool, optional
 369            Whether to keep the profile data after processing. Defaults to True.
 370
 371        Notes
 372        -----
 373        This method does the following:
 374        - calculates the noise threshold
 375        - does peak picking (creates mspeak_objs)
 376        - resets the mspeak_obj indexes
 377        """
 378
 379        # if runned mannually make sure to rerun filter_by_noise_threshold
 380        # calculates noise threshold
 381        # do peak picking( create mspeak_objs)
 382        # reset mspeak_obj the indexes
 383
 384        self.cal_noise_threshold()
 385
 386        self.find_peaks()
 387        self.reset_indexes()
 388
 389        if self.mspeaks:
 390            self._dynamic_range = self.max_abundance / self.min_abundance
 391        else:
 392            self._dynamic_range = 0
 393        if not keep_profile:
 394            self._abundance *= 0
 395            self._mz_exp *= 0
 396
 397    def _resolve_noise_threshold_method(self):
 398        """Resolve ``by_format`` and store the concrete noise-threshold method.
 399
 400        Profile spectra use ``log``. Centroid spectra use
 401        ``relative_abundance``. An explicit method is left unchanged.
 402        ``by_format`` is replaced on this spectrum's settings so saved
 403        parameters record what was used.
 404
 405        Returns
 406        -------
 407        str
 408            Concrete noise-threshold method.
 409        """
 410        method = self.settings.noise_threshold_method
 411        if method != "by_format":
 412            return method
 413        resolved = "relative_abundance" if self.is_centroid else "log"
 414        self.settings.noise_threshold_method = resolved
 415        return resolved
 416
 417    def cal_noise_threshold(self):
 418        """Calculate the noise threshold of the mass spectrum."""
 419
 420        self._resolve_noise_threshold_method()
 421
 422        if self.label == Labels.simulated_profile:
 423            self._baseline_noise, self._baseline_noise_std = 0.1, 1
 424
 425        if self.settings.noise_threshold_method == "log":
 426            self._baseline_noise, self._baseline_noise_std = (
 427                self.run_log_noise_threshold_calc()
 428            )
 429
 430        else:
 431            self._baseline_noise, self._baseline_noise_std = (
 432                self.run_noise_threshold_calc()
 433            )
 434
 435    @property
 436    def parameters(self):
 437        """Return the parameters of the mass spectrum."""
 438        return self._parameters
 439
 440    @parameters.setter
 441    def parameters(self, instance_MSParameters):
 442        self._parameters = instance_MSParameters
 443
 444    def set_parameter_from_json(self, parameters_path):
 445        """Set the parameters of the mass spectrum from a JSON file.
 446
 447        Parameters
 448        ----------
 449        parameters_path : str
 450            The path to the JSON file containing the parameters.
 451        """
 452        load_and_set_parameters_ms(self, parameters_path=parameters_path)
 453
 454    def set_parameter_from_toml(self, parameters_path):
 455        load_and_set_toml_parameters_ms(self, parameters_path=parameters_path)
 456
 457    @property
 458    def mspeaks_settings(self):
 459        """Return the MS peak settings of the mass spectrum."""
 460        return self.parameters.ms_peak
 461
 462    @mspeaks_settings.setter
 463    def mspeaks_settings(self, instance_MassSpecPeakSetting):
 464        self.parameters.ms_peak = instance_MassSpecPeakSetting
 465
 466    @property
 467    def settings(self):
 468        """Return the settings of the mass spectrum."""
 469        return self.parameters.mass_spectrum
 470
 471    @settings.setter
 472    def settings(self, instance_MassSpectrumSetting):
 473        self.parameters.mass_spectrum = instance_MassSpectrumSetting
 474
 475    @property
 476    def molecular_search_settings(self):
 477        """Return the molecular search settings of the mass spectrum."""
 478        return self.parameters.molecular_search
 479
 480    @molecular_search_settings.setter
 481    def molecular_search_settings(self, instance_MolecularFormulaSearchSettings):
 482        self.parameters.molecular_search = instance_MolecularFormulaSearchSettings
 483
 484    @property
 485    def mz_cal_profile(self):
 486        """Return the calibrated m/z profile of the mass spectrum."""
 487        return self._mz_cal_profile
 488
 489    @mz_cal_profile.setter
 490    def mz_cal_profile(self, mz_cal_list):
 491        if len(mz_cal_list) == len(self._mz_exp):
 492            self._mz_cal_profile = mz_cal_list
 493        else:
 494            raise Exception(
 495                "calibrated array (%i) is not of the same size of the data (%i)"
 496                % (len(mz_cal_list), len(self.mz_exp_profile))
 497            )
 498
 499    @property
 500    def mz_cal(self):
 501        """Return the calibrated m/z values of the mass spectrum."""
 502        return array([mspeak.mz_cal for mspeak in self.mspeaks])
 503
 504    @mz_cal.setter
 505    def mz_cal(self, mz_cal_list):
 506        if len(mz_cal_list) == len(self.mspeaks):
 507            self.is_calibrated = True
 508            for index, mz_cal in enumerate(mz_cal_list):
 509                self.mspeaks[index].mz_cal = mz_cal
 510        else:
 511            raise Exception(
 512                "calibrated array (%i) is not of the same size of the data (%i)"
 513                % (len(mz_cal_list), len(self._mspeaks))
 514            )
 515
 516    @property
 517    def mz_exp(self):
 518        """Return the experimental m/z values of the mass spectrum."""
 519        self.check_mspeaks()
 520
 521        if self.is_calibrated:
 522            return array([mspeak.mz_cal for mspeak in self.mspeaks])
 523
 524        else:
 525            return array([mspeak.mz_exp for mspeak in self.mspeaks])
 526
 527    @property
 528    def freq_exp_profile(self):
 529        """Return the experimental frequency profile of the mass spectrum."""
 530        return self._frequency_domain
 531
 532    @freq_exp_profile.setter
 533    def freq_exp_profile(self, new_data):
 534        self._frequency_domain = array(new_data)
 535
 536    @property
 537    def freq_exp_pp(self):
 538        """Return the experimental frequency values of the mass spectrum that are used for peak picking."""
 539        _, _, freq = self.prepare_peak_picking_data()
 540        return freq
 541
 542    @property
 543    def mz_exp_profile(self):
 544        """Return the experimental m/z profile of the mass spectrum."""
 545        if self.is_calibrated:
 546            return self.mz_cal_profile
 547        else:
 548            return self._mz_exp
 549
 550    @mz_exp_profile.setter
 551    def mz_exp_profile(self, new_data):
 552        self._mz_exp = array(new_data)
 553
 554    @property
 555    def mz_exp_pp(self):
 556        """Return the experimental m/z values of the mass spectrum that are used for peak picking."""
 557        mz, _, _ = self.prepare_peak_picking_data()
 558        return mz
 559
 560    @property
 561    def abundance_profile(self):
 562        """Return the abundance profile of the mass spectrum."""
 563        return self._abundance
 564
 565    @abundance_profile.setter
 566    def abundance_profile(self, new_data):
 567        self._abundance = array(new_data)
 568
 569    @property
 570    def abundance_profile_pp(self):
 571        """Return the abundance profile of the mass spectrum that is used for peak picking."""
 572        _, abundance, _ = self.prepare_peak_picking_data()
 573        return abundance
 574
 575    @property
 576    def abundance(self):
 577        """Return the abundance values of the mass spectrum."""
 578        self.check_mspeaks()
 579        return array([mspeak.abundance for mspeak in self.mspeaks])
 580
 581    def freq_exp(self):
 582        """Return the experimental frequency values of the mass spectrum."""
 583        self.check_mspeaks()
 584        return array([mspeak.freq_exp for mspeak in self.mspeaks])
 585
 586    @property
 587    def resolving_power(self):
 588        """Return the resolving power values of the mass spectrum."""
 589        self.check_mspeaks()
 590        return array([mspeak.resolving_power for mspeak in self.mspeaks])
 591
 592    @property
 593    def signal_to_noise(self):
 594        self.check_mspeaks()
 595        return array([mspeak.signal_to_noise for mspeak in self.mspeaks])
 596
 597    @property
 598    def nominal_mz(self):
 599        """Return the nominal m/z values of the mass spectrum."""
 600        if self._dict_nominal_masses_indexes:
 601            return sorted(list(self._dict_nominal_masses_indexes.keys()))
 602        else:
 603            raise ValueError("Nominal indexes not yet set")
 604
 605    def get_mz_and_abundance_peaks_tuples(self):
 606        """Return a list of tuples containing the m/z and abundance values of the mass spectrum."""
 607        self.check_mspeaks()
 608        return [(mspeak.mz_exp, mspeak.abundance) for mspeak in self.mspeaks]
 609
 610    @property
 611    def kmd(self):
 612        """Return the Kendrick mass defect values of the mass spectrum."""
 613        self.check_mspeaks()
 614        return array([mspeak.kmd for mspeak in self.mspeaks])
 615
 616    @property
 617    def kendrick_mass(self):
 618        """Return the Kendrick mass values of the mass spectrum."""
 619        self.check_mspeaks()
 620        return array([mspeak.kendrick_mass for mspeak in self.mspeaks])
 621
 622    @property
 623    def max_mz_exp(self):
 624        """Return the maximum experimental m/z value of the mass spectrum."""
 625        return max([mspeak.mz_exp for mspeak in self.mspeaks])
 626
 627    @property
 628    def min_mz_exp(self):
 629        """Return the minimum experimental m/z value of the mass spectrum."""
 630        return min([mspeak.mz_exp for mspeak in self.mspeaks])
 631
 632    @property
 633    def max_abundance(self):
 634        """Return the maximum abundance value of the mass spectrum."""
 635        return max([mspeak.abundance for mspeak in self.mspeaks])
 636
 637    @property
 638    def max_signal_to_noise(self):
 639        """Return the maximum signal-to-noise ratio of the mass spectrum."""
 640        return max([mspeak.signal_to_noise for mspeak in self.mspeaks])
 641
 642    @property
 643    def most_abundant_mspeak(self):
 644        """Return the most abundant MSpeak object of the mass spectrum."""
 645        return max(self.mspeaks, key=lambda m: m.abundance)
 646
 647    @property
 648    def min_abundance(self):
 649        """Return the minimum abundance value of the mass spectrum."""
 650        return min([mspeak.abundance for mspeak in self.mspeaks])
 651
 652    # takes too much cpu time
 653    @property
 654    def dynamic_range(self):
 655        """Return the dynamic range of the mass spectrum."""
 656        return self._dynamic_range
 657
 658    @property
 659    def baseline_noise(self):
 660        """Return the baseline noise of the mass spectrum."""
 661        if self._baseline_noise:
 662            return self._baseline_noise
 663        else:
 664            return None
 665
 666    @property
 667    def baseline_noise_std(self):
 668        """Return the standard deviation of the baseline noise of the mass spectrum."""
 669        if self._baseline_noise_std == 0:
 670            return self._baseline_noise_std
 671        if self._baseline_noise_std:
 672            return self._baseline_noise_std
 673        else:
 674            return None
 675
 676    @property
 677    def Aterm(self):
 678        """Return the A-term calibration coefficient of the mass spectrum."""
 679        return self._calibration_terms[0]
 680
 681    @property
 682    def Bterm(self):
 683        """Return the B-term calibration coefficient of the mass spectrum."""
 684        return self._calibration_terms[1]
 685
 686    @property
 687    def Cterm(self):
 688        """Return the C-term calibration coefficient of the mass spectrum."""
 689        return self._calibration_terms[2]
 690
 691    @property
 692    def filename(self):
 693        """Return the filename of the mass spectrum."""
 694        return Path(self._filename)
 695
 696    @property
 697    def dir_location(self):
 698        """Return the directory location of the mass spectrum."""
 699        return self._dir_location
 700
 701    def sort_by_mz(self):
 702        """Sort the mass spectrum by m/z values."""
 703        return sorted(self, key=lambda m: m.mz_exp)
 704
 705    def sort_by_abundance(self, reverse=False):
 706        """Sort the mass spectrum by abundance values."""
 707        return sorted(self, key=lambda m: m.abundance, reverse=reverse)
 708
 709    @property
 710    def tic(self):
 711        """Return the total ion current of the mass spectrum."""
 712        return trapezoid(self.abundance_profile, self.mz_exp_profile)
 713
 714    def check_mspeaks_warning(self):
 715        """Check if the mass spectrum has MSpeaks objects.
 716
 717        Raises
 718        ------
 719        Warning
 720            If the mass spectrum has no MSpeaks objects.
 721        """
 722        import warnings
 723
 724        if self.mspeaks:
 725            pass
 726        else:
 727            warnings.warn("mspeaks list is empty, continuing without filtering data")
 728
 729    def check_mspeaks(self):
 730        """Check if the mass spectrum has MSpeaks objects.
 731
 732        Raises
 733        ------
 734        Exception
 735            If the mass spectrum has no MSpeaks objects.
 736        """
 737        if self.mspeaks:
 738            pass
 739        else:
 740            raise Exception(
 741                "mspeaks list is empty, please run process_mass_spec() first"
 742            )
 743
 744    def remove_assignment_by_index(self, indexes):
 745        """Remove the molecular formula assignment of the MSpeaks objects at the specified indexes.
 746
 747        Parameters
 748        ----------
 749        indexes : list of int
 750            A list of indexes of the MSpeaks objects to remove the molecular formula assignment from.
 751        """
 752        for i in indexes:
 753            self.mspeaks[i].clear_molecular_formulas()
 754
 755    def filter_by_index(self, list_indexes):
 756        """Filter the mass spectrum by the specified indexes.
 757
 758        Parameters
 759        ----------
 760        list_indexes : list of int
 761            A list of indexes of the MSpeaks objects to drop.
 762
 763        """
 764
 765        self.mspeaks = [
 766            self.mspeaks[i] for i in range(len(self.mspeaks)) if i not in list_indexes
 767        ]
 768
 769        for i, mspeak in enumerate(self.mspeaks):
 770            mspeak.index = i
 771
 772        self._set_nominal_masses_start_final_indexes()
 773
 774    def filter_by_mz(self, min_mz, max_mz):
 775        """Filter the mass spectrum by the specified m/z range.
 776
 777        Parameters
 778        ----------
 779        min_mz : float
 780            The minimum m/z value to keep.
 781        max_mz : float
 782            The maximum m/z value to keep.
 783
 784        """
 785        self.check_mspeaks_warning()
 786        indexes = [
 787            index
 788            for index, mspeak in enumerate(self.mspeaks)
 789            if not min_mz <= mspeak.mz_exp <= max_mz
 790        ]
 791        self.filter_by_index(indexes)
 792
 793    def filter_by_s2n(self, min_s2n, max_s2n=False):
 794        """Filter the mass spectrum by the specified signal-to-noise ratio range.
 795
 796        Parameters
 797        ----------
 798        min_s2n : float
 799            The minimum signal-to-noise ratio to keep.
 800        max_s2n : float, optional
 801            The maximum signal-to-noise ratio to keep. Defaults to False (no maximum).
 802
 803        """
 804        self.check_mspeaks_warning()
 805        if max_s2n:
 806            indexes = [
 807                index
 808                for index, mspeak in enumerate(self.mspeaks)
 809                if not min_s2n <= mspeak.signal_to_noise <= max_s2n
 810            ]
 811        else:
 812            indexes = [
 813                index
 814                for index, mspeak in enumerate(self.mspeaks)
 815                if mspeak.signal_to_noise <= min_s2n
 816            ]
 817        self.filter_by_index(indexes)
 818
 819    def filter_by_abundance(self, min_abund, max_abund=False):
 820        """Filter the mass spectrum by the specified abundance range.
 821
 822        Parameters
 823        ----------
 824        min_abund : float
 825            The minimum abundance to keep.
 826        max_abund : float, optional
 827            The maximum abundance to keep. Defaults to False (no maximum).
 828
 829        """
 830        self.check_mspeaks_warning()
 831        if max_abund:
 832            indexes = [
 833                index
 834                for index, mspeak in enumerate(self.mspeaks)
 835                if not min_abund <= mspeak.abundance <= max_abund
 836            ]
 837        else:
 838            indexes = [
 839                index
 840                for index, mspeak in enumerate(self.mspeaks)
 841                if mspeak.abundance <= min_abund
 842            ]
 843        self.filter_by_index(indexes)
 844
 845    def filter_by_max_resolving_power(self, B, T):
 846        """Filter the mass spectrum by the specified maximum resolving power.
 847
 848        Parameters
 849        ----------
 850        B : float
 851        T : float
 852
 853        """
 854
 855        rpe = lambda m, z: (1.274e7 * z * B * T) / (m * z)
 856
 857        self.check_mspeaks_warning()
 858
 859        indexes_to_remove = [
 860            index
 861            for index, mspeak in enumerate(self.mspeaks)
 862            if mspeak.resolving_power >= rpe(mspeak.mz_exp, mspeak.polarity)
 863        ]
 864        self.filter_by_index(indexes_to_remove)
 865
 866    def filter_by_mean_resolving_power(
 867        self, ndeviations=3, plot=False, guess_pars=False
 868    ):
 869        """Filter the mass spectrum by the specified mean resolving power.
 870
 871        Parameters
 872        ----------
 873        ndeviations : float, optional
 874            The number of standard deviations to use for filtering. Defaults to 3.
 875        plot : bool, optional
 876            Whether to plot the resolving power distribution. Defaults to False.
 877        guess_pars : bool, optional
 878            Whether to guess the parameters for the Gaussian model. Defaults to False.
 879
 880        """
 881        self.check_mspeaks_warning()
 882        indexes_to_remove = MeanResolvingPowerFilter(
 883            self, ndeviations, plot, guess_pars
 884        ).main()
 885        self.filter_by_index(indexes_to_remove)
 886
 887    def filter_by_min_resolving_power(
 888        self, B, T, apodization_method: str = None, tolerance: float = 0
 889    ):
 890        """Filter the mass spectrum by the calculated minimum theoretical resolving power.
 891
 892        This is currently designed only for FTICR data, and accounts only for magnitude mode data
 893        Accurate results require passing the apodisaion method used to calculate the resolving power.
 894        see the ICRMassPeak function `resolving_power_calc` for more details.
 895
 896        Parameters
 897        ----------
 898        B : Magnetic field strength in Tesla, float
 899        T : transient length in seconds, float
 900        apodization_method : str, optional
 901            The apodization method to use for calculating the resolving power. Defaults to None.
 902        tolerance : float, optional
 903            The tolerance for the threshold. Defaults to 0, i.e. no tolerance
 904
 905        """
 906        if self.analyzer != "ICR":
 907            raise Exception("This method is only applicable to ICR mass spectra. ")
 908
 909        self.check_mspeaks_warning()
 910
 911        indexes_to_remove = [
 912            index
 913            for index, mspeak in enumerate(self.mspeaks)
 914            if mspeak.resolving_power
 915            < (1 - tolerance)
 916            * mspeak.resolving_power_calc(B, T, apodization_method=apodization_method)
 917        ]
 918        self.filter_by_index(indexes_to_remove)
 919
 920    def filter_by_noise_threshold(self):
 921        """Filter the mass spectrum by the noise threshold."""
 922
 923        threshold = self.get_noise_threshold()[1][0]
 924
 925        self.check_mspeaks_warning()
 926
 927        indexes_to_remove = [
 928            index
 929            for index, mspeak in enumerate(self.mspeaks)
 930            if mspeak.abundance <= threshold
 931        ]
 932        self.filter_by_index(indexes_to_remove)
 933
 934    def find_peaks(self):
 935        """Find the peaks of the mass spectrum."""
 936        # needs to clear previous results from peak_picking
 937        self._mspeaks = list()
 938
 939        # then do peak picking
 940        self.do_peak_picking()
 941        # print("A total of %i peaks were found" % len(self._mspeaks))
 942
 943    def change_kendrick_base_all_mspeaks(self, kendrick_dict_base):
 944        """Change the Kendrick base of all MSpeaks objects.
 945
 946        Parameters
 947        ----------
 948        kendrick_dict_base : dict
 949            A dictionary of the Kendrick base to change to.
 950
 951        Notes
 952        -----
 953        Example of kendrick_dict_base parameter: kendrick_dict_base = {"C": 1, "H": 2} or {"C": 1, "H": 1, "O":1} etc
 954        """
 955        self.parameters.ms_peak.kendrick_base = kendrick_dict_base
 956
 957        for mspeak in self.mspeaks:
 958            mspeak.change_kendrick_base(kendrick_dict_base)
 959
 960    def get_nominal_mz_first_last_indexes(self, nominal_mass):
 961        """Return the first and last indexes of the MSpeaks objects with the specified nominal mass.
 962
 963        Parameters
 964        ----------
 965        nominal_mass : int
 966            The nominal mass to get the indexes for.
 967
 968        Returns
 969        -------
 970        tuple
 971            A tuple containing the first and last indexes of the MSpeaks objects with the specified nominal mass.
 972        """
 973        if self._dict_nominal_masses_indexes:
 974            if nominal_mass in self._dict_nominal_masses_indexes.keys():
 975                return (
 976                    self._dict_nominal_masses_indexes.get(nominal_mass)[0],
 977                    self._dict_nominal_masses_indexes.get(nominal_mass)[1] + 1,
 978                )
 979
 980            else:
 981                # import warnings
 982                # uncomment warn to distribution
 983                # warnings.warn("Nominal mass not found in _dict_nominal_masses_indexes, returning (0, 0) for nominal mass %i"%nominal_mass)
 984                return (0, 0)
 985        else:
 986            raise Exception(
 987                "run process_mass_spec() function before trying to access the data"
 988            )
 989
 990    def get_masses_count_by_nominal_mass(self):
 991        """Return a dictionary of the nominal masses and their counts."""
 992
 993        dict_nominal_masses_count = {}
 994
 995        all_nominal_masses = list(set([i.nominal_mz_exp for i in self.mspeaks]))
 996
 997        for nominal_mass in all_nominal_masses:
 998            if nominal_mass not in dict_nominal_masses_count:
 999                dict_nominal_masses_count[nominal_mass] = len(
1000                    list(self.get_nominal_mass_indexes(nominal_mass))
1001                )
1002
1003        return dict_nominal_masses_count
1004
1005    def datapoints_count_by_nominal_mz(self, mz_overlay=0.1):
1006        """Return a dictionary of the nominal masses and their counts.
1007
1008        Parameters
1009        ----------
1010        mz_overlay : float, optional
1011            The m/z overlay to use for counting. Defaults to 0.1.
1012
1013        Returns
1014        -------
1015        dict
1016            A dictionary of the nominal masses and their counts.
1017        """
1018        dict_nominal_masses_count = {}
1019
1020        all_nominal_masses = list(set([i.nominal_mz_exp for i in self.mspeaks]))
1021
1022        for nominal_mass in all_nominal_masses:
1023            if nominal_mass not in dict_nominal_masses_count:
1024                min_mz = nominal_mass - mz_overlay
1025
1026                max_mz = nominal_mass + 1 + mz_overlay
1027
1028                indexes = indexes = where(
1029                    (self.mz_exp_profile > min_mz) & (self.mz_exp_profile < max_mz)
1030                )
1031
1032                dict_nominal_masses_count[nominal_mass] = indexes[0].size
1033
1034        return dict_nominal_masses_count
1035
1036    def get_nominal_mass_indexes(self, nominal_mass, overlay=0.1):
1037        """Return the indexes of the MSpeaks objects with the specified nominal mass.
1038
1039        Parameters
1040        ----------
1041        nominal_mass : int
1042            The nominal mass to get the indexes for.
1043        overlay : float, optional
1044            The m/z overlay to use for counting. Defaults to 0.1.
1045
1046        Returns
1047        -------
1048        generator
1049            A generator of the indexes of the MSpeaks objects with the specified nominal mass.
1050        """
1051        min_mz_to_look = nominal_mass - overlay
1052        max_mz_to_look = nominal_mass + 1 + overlay
1053
1054        return (
1055            i
1056            for i in range(len(self.mspeaks))
1057            if min_mz_to_look <= self.mspeaks[i].mz_exp <= max_mz_to_look
1058        )
1059
1060        # indexes = (i for i in range(len(self.mspeaks)) if min_mz_to_look <= self.mspeaks[i].mz_exp <= max_mz_to_look)
1061        # return indexes
1062
1063    def _set_nominal_masses_start_final_indexes(self):
1064        """Set the start and final indexes of the MSpeaks objects for all nominal masses."""
1065        dict_nominal_masses_indexes = {}
1066
1067        all_nominal_masses = set(i.nominal_mz_exp for i in self.mspeaks)
1068
1069        for nominal_mass in all_nominal_masses:
1070            # indexes = self.get_nominal_mass_indexes(nominal_mass)
1071            # Convert the iterator to a list to avoid multiple calls
1072            indexes = list(self.get_nominal_mass_indexes(nominal_mass))
1073
1074            # If the list is not empty, find the first and last; otherwise, set None
1075            if indexes:
1076                first, last = indexes[0], indexes[-1]
1077            else:
1078                first = last = None
1079            # defaultvalue = None
1080            # first = last = next(indexes, defaultvalue)
1081            # for last in indexes:
1082            #    pass
1083
1084            dict_nominal_masses_indexes[nominal_mass] = (first, last)
1085
1086        self._dict_nominal_masses_indexes = dict_nominal_masses_indexes
1087
1088    def plot_centroid(self, ax=None, c="g"):
1089        """Plot the centroid data of the mass spectrum.
1090
1091        Parameters
1092        ----------
1093        ax : matplotlib.axes.Axes, optional
1094            The matplotlib axes to plot on. Defaults to None.
1095        c : str, optional
1096            The color to use for the plot. Defaults to 'g' (green).
1097
1098        Returns
1099        -------
1100        matplotlib.axes.Axes
1101            The matplotlib axes containing the plot.
1102
1103        Raises
1104        ------
1105        Exception
1106            If no centroid data is found.
1107        """
1108
1109        import matplotlib.pyplot as plt
1110
1111        if self._mspeaks:
1112            if ax is None:
1113                ax = plt.gca()
1114
1115            markerline_a, stemlines_a, baseline_a = ax.stem(
1116                self.mz_exp, self.abundance, linefmt="-", markerfmt=" "
1117            )
1118
1119            plt.setp(markerline_a, "color", c, "linewidth", 2)
1120            plt.setp(stemlines_a, "color", c, "linewidth", 2)
1121            plt.setp(baseline_a, "color", c, "linewidth", 2)
1122
1123            ax.set_xlabel("$\t{m/z}$", fontsize=12)
1124            ax.set_ylabel("Abundance", fontsize=12)
1125            ax.tick_params(axis="both", which="major", labelsize=12)
1126
1127            ax.axes.spines["top"].set_visible(False)
1128            ax.axes.spines["right"].set_visible(False)
1129
1130            ax.get_yaxis().set_visible(False)
1131            ax.spines["left"].set_visible(False)
1132
1133        else:
1134            raise Exception("No centroid data found, please run process_mass_spec")
1135
1136        return ax
1137
1138    def plot_profile_and_noise_threshold(self, ax=None, legend=False):
1139        """Plot the profile data and noise threshold of the mass spectrum.
1140
1141        Parameters
1142        ----------
1143        ax : matplotlib.axes.Axes, optional
1144            The matplotlib axes to plot on. Defaults to None.
1145        legend : bool, optional
1146            Whether to show the legend. Defaults to False.
1147
1148        Returns
1149        -------
1150        matplotlib.axes.Axes
1151            The matplotlib axes containing the plot.
1152
1153        Raises
1154        ------
1155        Exception
1156            If no noise threshold is found.
1157        """
1158        import matplotlib.pyplot as plt
1159
1160        if self.baseline_noise_std and self.baseline_noise_std:
1161            # x = (self.mz_exp_profile.min(), self.mz_exp_profile.max())
1162            baseline = (self.baseline_noise, self.baseline_noise)
1163
1164            # std = self.parameters.mass_spectrum.noise_threshold_min_std
1165            # threshold = self.baseline_noise_std + (std * self.baseline_noise_std)
1166            x, y = self.get_noise_threshold()
1167
1168            if ax is None:
1169                ax = plt.gca()
1170
1171            ax.plot(
1172                self.mz_exp_profile,
1173                self.abundance_profile,
1174                color="green",
1175                label="Spectrum",
1176            )
1177            ax.plot(x, (baseline, baseline), color="yellow", label="Baseline Noise")
1178            ax.plot(x, y, color="red", label="Noise Threshold")
1179
1180            ax.set_xlabel("$\t{m/z}$", fontsize=12)
1181            ax.set_ylabel("Abundance", fontsize=12)
1182            ax.tick_params(axis="both", which="major", labelsize=12)
1183
1184            ax.axes.spines["top"].set_visible(False)
1185            ax.axes.spines["right"].set_visible(False)
1186
1187            ax.get_yaxis().set_visible(False)
1188            ax.spines["left"].set_visible(False)
1189            if legend:
1190                ax.legend()
1191
1192        else:
1193            raise Exception("Calculate noise threshold first")
1194
1195        return ax
1196
1197    def plot_mz_domain_profile(self, color="green", ax=None):
1198        """Plot the m/z domain profile of the mass spectrum.
1199
1200        Parameters
1201        ----------
1202        color : str, optional
1203            The color to use for the plot. Defaults to 'green'.
1204        ax : matplotlib.axes.Axes, optional
1205            The matplotlib axes to plot on. Defaults to None.
1206
1207        Returns
1208        -------
1209        matplotlib.axes.Axes
1210            The matplotlib axes containing the plot.
1211        """
1212
1213        import matplotlib.pyplot as plt
1214
1215        if ax is None:
1216            ax = plt.gca()
1217        ax.plot(self.mz_exp_profile, self.abundance_profile, color=color)
1218        ax.set(xlabel="m/z", ylabel="abundance")
1219
1220        return ax
1221
1222    def to_excel(self, out_file_path, write_metadata=True):
1223        """Export the mass spectrum to an Excel file.
1224
1225        Parameters
1226        ----------
1227        out_file_path : str
1228            The path to the Excel file to export to.
1229        write_metadata : bool, optional
1230            Whether to write the metadata to the Excel file. Defaults to True.
1231
1232        Returns
1233        -------
1234        None
1235        """
1236        from corems.mass_spectrum.output.export import HighResMassSpecExport
1237
1238        exportMS = HighResMassSpecExport(out_file_path, self)
1239        exportMS.to_excel(write_metadata=write_metadata)
1240
1241    def to_hdf(self, out_file_path):
1242        """Export the mass spectrum to an HDF file.
1243
1244        Parameters
1245        ----------
1246        out_file_path : str
1247            The path to the HDF file to export to.
1248
1249        Returns
1250        -------
1251        None
1252        """
1253        from corems.mass_spectrum.output.export import HighResMassSpecExport
1254
1255        exportMS = HighResMassSpecExport(out_file_path, self)
1256        exportMS.to_hdf()
1257
1258    def to_csv(self, out_file_path, write_metadata=True):
1259        """Export the mass spectrum to a CSV file.
1260
1261        Parameters
1262        ----------
1263        out_file_path : str
1264            The path to the CSV file to export to.
1265        write_metadata : bool, optional
1266            Whether to write the metadata to the CSV file. Defaults to True.
1267
1268        """
1269        from corems.mass_spectrum.output.export import HighResMassSpecExport
1270
1271        exportMS = HighResMassSpecExport(out_file_path, self)
1272        exportMS.to_csv(write_metadata=write_metadata)
1273
1274    def to_pandas(self, out_file_path, write_metadata=True):
1275        """Export the mass spectrum to a Pandas dataframe with pkl extension.
1276
1277        Parameters
1278        ----------
1279        out_file_path : str
1280            The path to the CSV file to export to.
1281        write_metadata : bool, optional
1282            Whether to write the metadata to the CSV file. Defaults to True.
1283
1284        """
1285        from corems.mass_spectrum.output.export import HighResMassSpecExport
1286
1287        exportMS = HighResMassSpecExport(out_file_path, self)
1288        exportMS.to_pandas(write_metadata=write_metadata)
1289
1290    def to_parquet(self, out_file_path, write_metadata=True):
1291        """Export the mass spectrum to a Parquet file.
1292
1293        Parameters
1294        ----------
1295        out_file_path : str
1296            The path to the Parquet file to export to.
1297        write_metadata : bool, optional
1298            Whether to write the metadata to a JSON file. Defaults to True.
1299
1300        Notes
1301        -----
1302        Requires ``pyarrow`` (or ``fastparquet``) for pandas Parquet I/O.
1303        """
1304        from corems.mass_spectrum.output.export import HighResMassSpecExport
1305
1306        exportMS = HighResMassSpecExport(out_file_path, self)
1307        exportMS.to_parquet(write_metadata=write_metadata)
1308
1309    def to_dataframe(self, additional_columns=None):
1310        """Return the mass spectrum as a Pandas dataframe.
1311
1312        Parameters
1313        ----------
1314        additional_columns : list, optional
1315            A list of additional columns to include in the dataframe. Defaults to None.
1316            Suitable columns are: "Aromaticity Index", "Aromaticity Index (modified)", and "NOSC"
1317
1318        Returns
1319        -------
1320        pandas.DataFrame
1321            The mass spectrum as a Pandas dataframe.
1322        """
1323        from corems.mass_spectrum.output.export import HighResMassSpecExport
1324
1325        exportMS = HighResMassSpecExport(self.filename, self)
1326        return exportMS.get_pandas_df(additional_columns=additional_columns)
1327
1328    def to_json(self):
1329        """Return the mass spectrum as a JSON file."""
1330        from corems.mass_spectrum.output.export import HighResMassSpecExport
1331
1332        exportMS = HighResMassSpecExport(self.filename, self)
1333        return exportMS.to_json()
1334
1335    def parameters_json(self):
1336        """Return the parameters of the mass spectrum as a JSON string."""
1337        from corems.mass_spectrum.output.export import HighResMassSpecExport
1338
1339        exportMS = HighResMassSpecExport(self.filename, self)
1340        return exportMS.parameters_to_json()
1341
1342    def parameters_toml(self):
1343        """Return the parameters of the mass spectrum as a TOML string."""
1344        from corems.mass_spectrum.output.export import HighResMassSpecExport
1345
1346        exportMS = HighResMassSpecExport(self.filename, self)
1347        return exportMS.parameters_to_toml()
1348
1349
1350class MassSpecProfile(MassSpecBase):
1351    """A mass spectrum class when the entry point is on profile format
1352
1353    Notes
1354    -----
1355    Stores the profile data and instrument settings.
1356    Iteration over a list of MSPeaks classes stored at the _mspeaks attributes.
1357    _mspeaks is populated under the hood by calling process_mass_spec method.
1358    Iteration is null if _mspeaks is empty. Many more attributes and methods inherited from MassSpecBase().
1359
1360    Parameters
1361    ----------
1362    data_dict : dict
1363        A dictionary containing the profile data.
1364    d_params : dict{'str': float, int or str}
1365        contains the instrument settings and processing settings
1366    auto_process : bool, optional
1367        Whether to automatically process the mass spectrum. Defaults to True.
1368
1369
1370    Attributes
1371    ----------
1372    _abundance : ndarray
1373        The abundance values of the mass spectrum.
1374    _mz_exp : ndarray
1375        The m/z values of the mass spectrum.
1376    _mspeaks : list
1377        A list of mass peaks.
1378
1379    Methods
1380    ----------
1381    * process_mass_spec(). Process the mass spectrum.
1382
1383    see also: MassSpecBase(), MassSpecfromFreq(), MassSpecCentroid()
1384    """
1385
1386    def __init__(self, data_dict, d_params, auto_process=True):
1387        # print(data_dict.keys())
1388        super().__init__(
1389            data_dict.get(Labels.mz), data_dict.get(Labels.abundance), d_params
1390        )
1391
1392        if auto_process:
1393            self.process_mass_spec()
1394
1395
1396class MassSpecfromFreq(MassSpecBase):
1397    """A mass spectrum class when data entry is on frequency domain
1398
1399    Notes
1400    -----
1401    - Transform to m/z based on the settings stored at d_params
1402    - Stores the profile data and instrument settings
1403    - Iteration over a list of MSPeaks classes stored at the _mspeaks attributes
1404    - _mspeaks is populated under the hood by calling process_mass_spec method
1405    - iteration is null if _mspeaks is empty
1406
1407    Parameters
1408    ----------
1409    frequency_domain : list(float)
1410        all datapoints in frequency domain in Hz
1411    magnitude :  frequency_domain : list(float)
1412        all datapoints in for magnitude of each frequency datapoint
1413    d_params : dict{'str': float, int or str}
1414        contains the instrument settings and processing settings
1415    auto_process : bool, optional
1416        Whether to automatically process the mass spectrum. Defaults to True.
1417    keep_profile : bool, optional
1418        Whether to keep the profile data. Defaults to True.
1419
1420    Attributes
1421    ----------
1422    has_frequency : bool
1423        Whether the mass spectrum has frequency data.
1424    _frequency_domain : list(float)
1425        Frequency domain in Hz
1426    label : str
1427        store label (Bruker, Midas Transient, see Labels class ). It across distinct processing points
1428    _abundance : ndarray
1429        The abundance values of the mass spectrum.
1430    _mz_exp : ndarray
1431        The m/z values of the mass spectrum.
1432    _mspeaks : list
1433        A list of mass peaks.
1434    See Also: all the attributes of MassSpecBase class
1435
1436    Methods
1437    ----------
1438    * _set_mz_domain().
1439        calculates the m_z based on the setting of d_params
1440    * process_mass_spec().  Process the mass spectrum.
1441
1442    see also: MassSpecBase(), MassSpecProfile(), MassSpecCentroid()
1443    """
1444
1445    def __init__(
1446        self,
1447        frequency_domain,
1448        magnitude,
1449        d_params,
1450        auto_process=True,
1451        keep_profile=True,
1452    ):
1453        super().__init__(None, magnitude, d_params)
1454
1455        self._frequency_domain = frequency_domain
1456        self.has_frequency = True
1457        self._set_mz_domain()
1458        self._sort_mz_domain()
1459
1460        self.magnetron_frequency = None
1461        self.magnetron_frequency_sigma = None
1462
1463        # use this call to automatically process data as the object is created, Setting need to be changed before initiating the class to be in effect
1464
1465        if auto_process:
1466            self.process_mass_spec(keep_profile=keep_profile)
1467
1468    def _sort_mz_domain(self):
1469        """Sort the mass spectrum by m/z values."""
1470
1471        if self._mz_exp[0] > self._mz_exp[-1]:
1472            self._mz_exp = self._mz_exp[::-1]
1473            self._abundance = self._abundance[::-1]
1474            self._frequency_domain = self._frequency_domain[::-1]
1475
1476    def _set_mz_domain(self):
1477        """Set the m/z domain of the mass spectrum based on the settings of d_params."""
1478        if self.label == Labels.bruker_frequency:
1479            self._mz_exp = self._f_to_mz_bruker()
1480
1481        else:
1482            self._mz_exp = self._f_to_mz()
1483
1484    @property
1485    def transient_settings(self):
1486        """Return the transient settings of the mass spectrum."""
1487        return self.parameters.transient
1488
1489    @transient_settings.setter
1490    def transient_settings(self, instance_TransientSetting):
1491        self.parameters.transient = instance_TransientSetting
1492
1493    def calc_magnetron_freq(self, max_magnetron_freq=50, magnetron_freq_bins=300):
1494        """Calculates the magnetron frequency of the mass spectrum.
1495
1496        Parameters
1497        ----------
1498        max_magnetron_freq : float, optional
1499            The maximum magnetron frequency. Defaults to 50.
1500        magnetron_freq_bins : int, optional
1501            The number of bins to use for the histogram. Defaults to 300.
1502
1503        Returns
1504        -------
1505        None
1506
1507        Notes
1508        -----
1509        Calculates the magnetron frequency by examining all the picked peaks and the distances between them in the frequency domain.
1510        A histogram of those values below the threshold 'max_magnetron_freq' with the 'magnetron_freq_bins' number of bins is calculated.
1511        A gaussian model is fit to this histogram - the center value of this (statistically probably) the magnetron frequency.
1512        This appears to work well or nOmega datasets, but may not work well for 1x datasets or those with very low magnetron peaks.
1513        """
1514        ms_df = DataFrame(self.freq_exp(), columns=["Freq"])
1515        ms_df["FreqDelta"] = ms_df["Freq"].diff()
1516
1517        freq_hist = histogram(
1518            ms_df[ms_df["FreqDelta"] < max_magnetron_freq]["FreqDelta"],
1519            bins=magnetron_freq_bins,
1520        )
1521
1522        mod = GaussianModel()
1523        pars = mod.guess(freq_hist[0], x=freq_hist[1][:-1])
1524        out = mod.fit(freq_hist[0], pars, x=freq_hist[1][:-1])
1525        self.magnetron_frequency = out.best_values["center"]
1526        self.magnetron_frequency_sigma = out.best_values["sigma"]
1527
1528
1529class MassSpecCentroid(MassSpecBase):
1530    """A mass spectrum class when the entry point is on centroid format
1531
1532    Notes
1533    -----
1534    - Stores the centroid data and instrument settings
1535    - Simulate profile data based on Gaussian or Lorentzian peak shape
1536    - Iteration over a list of MSPeaks classes stored at the _mspeaks attributes
1537    - _mspeaks is populated under the hood by calling process_mass_spec method
1538    - iteration is null if _mspeaks is empty
1539
1540    Parameters
1541    ----------
1542    data_dict : dict {string: numpy array float64 )
1543        contains keys [m/z, Abundance, Resolving Power, S/N]
1544    d_params : dict{'str': float, int or str}
1545        contains the instrument settings and processing settings
1546    auto_process : bool, optional
1547        Whether to automatically process the mass spectrum. Defaults to True.
1548
1549    Attributes
1550    ----------
1551    label : str
1552        store label (Bruker, Midas Transient, see Labels class)
1553    _baseline_noise : float
1554        store baseline noise
1555    _baseline_noise_std : float
1556        store baseline noise std
1557    _abundance : ndarray
1558        The abundance values of the mass spectrum.
1559    _mz_exp : ndarray
1560        The m/z values of the mass spectrum.
1561    _mspeaks : list
1562        A list of mass peaks.
1563
1564
1565    Methods
1566    ----------
1567    * process_mass_spec().
1568        Process the mass spectrum. Overriden from MassSpecBase. Populates the _mspeaks list with MSpeaks class using the centroid data.
1569    * __simulate_profile__data__().
1570        Simulate profile data based on Gaussian or Lorentzian peak shape. Needs theoretical resolving power calculation and define peak shape, intended for plotting and inspection purposes only.
1571
1572    see also: MassSpecBase(), MassSpecfromFreq(), MassSpecProfile()
1573    """
1574
1575    def __init__(self, data_dict, d_params, auto_process=True):
1576        super().__init__([], [], d_params)
1577
1578        self._set_parameters_objects(d_params)
1579
1580        if self.label == Labels.thermo_centroid:
1581            self._baseline_noise = d_params.get("baseline_noise")
1582            self._baseline_noise_std = d_params.get("baseline_noise_std")
1583
1584        self.is_centroid = True
1585        self.data_dict = data_dict
1586        self._mz_exp = data_dict[Labels.mz]
1587        self._abundance = data_dict[Labels.abundance]
1588
1589        if auto_process:
1590            self.process_mass_spec()
1591
1592    def __simulate_profile__data__(self, exp_mz_centroid, magnitude_centroid):
1593        """Simulate profile data based on Gaussian or Lorentzian peak shape
1594
1595        Notes
1596        -----
1597        Needs theoretical resolving power calculation and define peak shape.
1598        This is a quick fix to trick a line plot be able to plot as sticks for plotting and inspection purposes only.
1599
1600        Parameters
1601        ----------
1602        exp_mz_centroid : list(float)
1603            list of m/z values
1604        magnitude_centroid : list(float)
1605            list of abundance values
1606
1607
1608        Returns
1609        -------
1610        x : list(float)
1611            list of m/z values
1612        y : list(float)
1613            list of abundance values
1614        """
1615
1616        x, y = [], []
1617        for i in range(len(exp_mz_centroid)):
1618            x.append(exp_mz_centroid[i] - 0.0000001)
1619            x.append(exp_mz_centroid[i])
1620            x.append(exp_mz_centroid[i] + 0.0000001)
1621            y.append(0)
1622            y.append(magnitude_centroid[i])
1623            y.append(0)
1624        return x, y
1625
1626    @property
1627    def mz_exp_profile(self):
1628        """Return the m/z profile of the mass spectrum."""
1629        mz_list = []
1630        for mz in self.mz_exp:
1631            mz_list.append(mz - 0.0000001)
1632            mz_list.append(mz)
1633            mz_list.append(mz + 0.0000001)
1634        return mz_list
1635
1636    @mz_exp_profile.setter
1637    def mz_exp_profile(self, _mz_exp):
1638        self._mz_exp = _mz_exp
1639
1640    @property
1641    def abundance_profile(self):
1642        """Return the abundance profile of the mass spectrum."""
1643        ab_list = []
1644        for ab in self.abundance:
1645            ab_list.append(0)
1646            ab_list.append(ab)
1647            ab_list.append(0)
1648        return ab_list
1649
1650    @abundance_profile.setter
1651    def abundance_profile(self, abundance):
1652        self._abundance = abundance
1653
1654    @property
1655    def tic(self):
1656        """Return the total ion current of the mass spectrum."""
1657        return sum(self.abundance)
1658
1659    def process_mass_spec(self):
1660        """Process the mass spectrum."""
1661
1662        # overwrite process_mass_spec
1663        # mspeak objs are usually added inside the PeaKPicking class
1664        # for profile and freq based data
1665        data_dict = self.data_dict
1666        # Peak polarity only (±1); formula assignment charge is separate.
1667        peak_polarity = self.polarity
1668
1669        # Check if resolving power is present
1670        rp_present = True
1671        if not data_dict.get(Labels.rp):
1672            rp_present = False
1673        if rp_present and list(data_dict.get(Labels.rp)) == [None] * len(
1674            data_dict.get(Labels.rp)
1675        ):
1676            rp_present = False
1677
1678        # Check if s2n is present
1679        s2n_present = True
1680        if not data_dict.get(Labels.s2n):
1681            s2n_present = False
1682        if s2n_present and list(data_dict.get(Labels.s2n)) == [None] * len(
1683            data_dict.get(Labels.s2n)
1684        ):
1685            s2n_present = False
1686
1687        # Warning if no s2n data but noise thresholding is set to signal_noise
1688        if (
1689            not s2n_present
1690            and self.parameters.mass_spectrum.noise_threshold_method == "signal_noise"
1691        ):
1692            raise Exception("Signal to Noise data is missing for noise thresholding")
1693
1694        # Pull out abundance data
1695        abun = array(data_dict.get(Labels.abundance)).astype(float)
1696
1697        # Get the threshold for filtering if using minima, relative, or absolute abundance thresholding
1698        abundance_threshold, factor = self.get_threshold(abun)
1699
1700        # Set rp_i and s2n_i to None which will be overwritten if present
1701        rp_i, s2n_i = np.nan, np.nan
1702        for index, mz in enumerate(data_dict.get(Labels.mz)):
1703            if rp_present:
1704                if not data_dict.get(Labels.rp)[index]:
1705                    rp_i = np.nan
1706                else:
1707                    rp_i = float(data_dict.get(Labels.rp)[index])
1708            if s2n_present:
1709                if not data_dict.get(Labels.s2n)[index]:
1710                    s2n_i = np.nan
1711                else:
1712                    s2n_i = float(data_dict.get(Labels.s2n)[index])
1713
1714            # centroid peak does not have start and end peak index pos
1715            massspec_indexes = (index, index, index)
1716
1717            # Add peaks based on the noise thresholding method
1718            if (
1719                self.parameters.mass_spectrum.noise_threshold_method
1720                in ["minima", "relative_abundance", "absolute_abundance"]
1721                and abun[index] / factor >= abundance_threshold
1722            ):
1723                self.add_mspeak(
1724                    peak_polarity,
1725                    mz,
1726                    abun[index],
1727                    rp_i,
1728                    s2n_i,
1729                    massspec_indexes,
1730                    ms_parent=self,
1731                )
1732            if (
1733                self.parameters.mass_spectrum.noise_threshold_method == "signal_noise"
1734                and s2n_i >= self.parameters.mass_spectrum.noise_threshold_min_s2n
1735            ):
1736                self.add_mspeak(
1737                    peak_polarity,
1738                    mz,
1739                    abun[index],
1740                    rp_i,
1741                    s2n_i,
1742                    massspec_indexes,
1743                    ms_parent=self,
1744                )
1745
1746        self.mspeaks = self._mspeaks
1747        self._dynamic_range = self.max_abundance / self.min_abundance
1748        self._set_nominal_masses_start_final_indexes()
1749
1750        if self.label != Labels.thermo_centroid:
1751            if self.settings.noise_threshold_method == "log":
1752                raise Exception("log noise Not tested for centroid data")
1753                # self._baseline_noise, self._baseline_noise_std = self.run_log_noise_threshold_calc()
1754
1755            else:
1756                self._baseline_noise, self._baseline_noise_std = (
1757                    self.run_noise_threshold_calc()
1758                )
1759
1760        del self.data_dict
1761
1762    def to_profile(
1763        self,
1764        mz_profile,
1765        abundance_profile,
1766        mz_cal_profile=None,
1767        peak_profile_indices=None,
1768        is_calibrated=False,
1769    ):
1770        """Return a :class:`MassSpecProfile` sharing this spectrum's peaks.
1771
1772        Reuses the already-picked peaks (and any molecular formula
1773        assignments on them) instead of re-picking from the profile data, and
1774        re-points each peak's parent reference to the returned object. This
1775        transfers ownership of the peaks: ``self.mspeaks``/``self._mspeaks``
1776        are cleared once the new object holds them, so this spectrum and the
1777        returned one never end up sharing (and independently mutating) the
1778        same peak list. Used to rehydrate a profile mass spectrum from a
1779        CoreMS HDF5 file that stored both a profile and a peak table,
1780        without re-running peak picking.
1781
1782        Parameters
1783        ----------
1784        mz_profile : array_like
1785            Experimental (uncalibrated) m/z profile.
1786        abundance_profile : array_like
1787            Abundance profile, same length as ``mz_profile``.
1788        mz_cal_profile : array_like, optional
1789            Calibrated m/z profile, same length as ``mz_profile``. Applied
1790            whenever given, independently of ``is_calibrated``. Default is
1791            None.
1792        peak_profile_indices : array_like of shape (n_peaks, 3), optional
1793            Per-peak ``(left, apex, right)`` indices into
1794            ``mz_profile``/``abundance_profile``, in the same order as
1795            ``self.mspeaks``. Ignored entirely if its length doesn't match
1796            ``self.mspeaks``. Otherwise checked per peak -- a stored triple
1797            is only trusted when ``0 <= left <= apex <= right < len(mz_profile)``
1798            -- since a file can carry indices from before a later
1799            reprocessing changed the profile length (the writer skips
1800            storing new ones when that happens, but an older file may still
1801            have a stale set from before this check existed). When omitted,
1802            or for any one peak whose stored triple fails that check, that
1803            peak's bounds are estimated instead: the apex starts at the
1804            nearest profile index to the peak's own m/z (its calibrated m/z
1805            against ``mz_cal_profile`` if the peak is calibrated, else its
1806            experimental m/z against ``mz_profile`` -- either array may run
1807            ascending or descending), the left/right bounds are then the
1808            neighboring local minima from that starting point, and the apex
1809            is finally set to the actual local maximum between those two
1810            bounds -- not the initial nearest-index estimate, which can be
1811            off by a few points from the quadratic-fit interpolation the
1812            original picking used.
1813        is_calibrated : bool, optional
1814            Whether the stored m/z is calibrated. Sets ``is_calibrated`` on
1815            the returned spectrum directly; ``mz_cal_profile`` is applied
1816            independently of this flag whenever it is given. Default is
1817            False.
1818
1819        Returns
1820        -------
1821        MassSpecProfile
1822        """
1823        d_params = {
1824            "Aterm": self.Aterm,
1825            "Bterm": self.Bterm,
1826            "Cterm": self.Cterm,
1827            Labels.label: self.label,
1828            "analyzer": self.analyzer,
1829            "acquisition_time": self.acquisition_time,
1830            "instrument_label": self.instrument_label,
1831            "polarity": self.polarity,
1832            "scan_number": self.scan_number,
1833            "rt": self.retention_time,
1834            "mobility_rt": self.mobility_rt,
1835            "mobility_scan": self.mobility_scan,
1836            "filename_path": self._filename,
1837            "dir_location": self._dir_location,
1838            "baseline_noise": self._baseline_noise,
1839            "baseline_noise_std": self._baseline_noise_std,
1840            "sample_name": self.sample_name,
1841        }
1842
1843        new_obj = MassSpecProfile(
1844            {Labels.mz: mz_profile, Labels.abundance: abundance_profile},
1845            d_params,
1846            auto_process=False,
1847        )
1848
1849        new_obj._mspeaks = self._mspeaks
1850        new_obj.mspeaks = self.mspeaks
1851        new_obj._dict_nominal_masses_indexes = self._dict_nominal_masses_indexes
1852        new_obj._dynamic_range = self._dynamic_range
1853        # MassSpecBase.__init__ gave new_obj a fresh, default MSParameters;
1854        # replace it with this spectrum's (e.g. settings restored by a
1855        # reader's load_settings()) so processing settings survive the
1856        # conversion.
1857        new_obj.parameters = self.parameters
1858
1859        # Transfer ownership: this spectrum no longer holds these peaks.
1860        self._mspeaks = []
1861        self.mspeaks = []
1862
1863        if is_calibrated:
1864            new_obj.is_calibrated = True
1865        if mz_cal_profile is not None:
1866            new_obj.mz_cal_profile = mz_cal_profile
1867
1868        abundance_arr = asarray(new_obj.abundance_profile)
1869        mz_arr = asarray(mz_profile)
1870        mz_cal_arr = asarray(mz_cal_profile) if mz_cal_profile is not None else None
1871        n = len(abundance_arr)
1872
1873        if peak_profile_indices is not None and len(peak_profile_indices) != len(
1874            new_obj.mspeaks
1875        ):
1876            peak_profile_indices = None
1877
1878        for i, mspeak in enumerate(new_obj.mspeaks):
1879            mspeak._ms_parent = new_obj
1880
1881            stored = None
1882            if peak_profile_indices is not None:
1883                left, apex, right = (int(v) for v in peak_profile_indices[i])
1884                # Stored indices can predate a later reprocessing with a
1885                # different-length profile (the writer skips storing them
1886                # when that happens, but an older file may still have a
1887                # stale set); only trust them if they're a sane bound into
1888                # *this* profile, else fall back to the search below.
1889                if 0 <= left <= apex <= right < n:
1890                    stored = (left, apex, right)
1891
1892            if stored is not None:
1893                left, apex, right = stored
1894            else:
1895                if mspeak.mz_cal is not None and mz_cal_arr is not None:
1896                    search_profile, search_value = mz_cal_arr, mspeak.mz_cal
1897                else:
1898                    search_profile, search_value = mz_arr, mspeak._mz_exp
1899
1900                apex = _nearest_index_either_order(search_profile, search_value)
1901                right = new_obj.find_minima(apex, abundance_arr, n, right=True)
1902                left = new_obj.find_minima(apex, abundance_arr, n, right=False)
1903                # Re-anchor the apex to the tallest sample within the now-known
1904                # bounds, rather than trusting the nearest-index estimate above.
1905                apex = left + int(argmax(abundance_arr[left : right + 1]))
1906            mspeak.peak_left_index = left
1907            mspeak.peak_apex_index = apex
1908            mspeak.peak_right_index = right
1909
1910        return new_obj
1911
1912
1913class MassSpecCentroidLowRes(MassSpecCentroid):
1914    """A mass spectrum class when the entry point is on low resolution centroid format
1915
1916    Notes
1917    -----
1918    Does not store MSPeak Objs, will iterate over mz, abundance pairs instead
1919
1920    Parameters
1921    ----------
1922    data_dict : dict {string: numpy array float64 )
1923        contains keys [m/z, Abundance, Resolving Power, S/N]
1924    d_params : dict{'str': float, int or str}
1925        contains the instrument settings and processing settings
1926
1927    Attributes
1928    ----------
1929    _processed_tic : float
1930        store processed total ion current
1931    _abundance : ndarray
1932        The abundance values of the mass spectrum.
1933    _mz_exp : ndarray
1934        The m/z values of the mass spectrum.
1935    """
1936
1937    def __init__(self, data_dict, d_params):
1938        self._set_parameters_objects(d_params)
1939        self._mz_exp = array(data_dict.get(Labels.mz))
1940        self._abundance = array(data_dict.get(Labels.abundance))
1941        self._processed_tic = None
1942
1943    def __len__(self):
1944        return len(self.mz_exp)
1945
1946    def __getitem__(self, position):
1947        return (self.mz_exp[position], self.abundance[position])
1948
1949    @property
1950    def mz_exp(self):
1951        """Return the m/z values of the mass spectrum."""
1952        return self._mz_exp
1953
1954    @property
1955    def abundance(self):
1956        """Return the abundance values of the mass spectrum."""
1957        return self._abundance
1958
1959    @property
1960    def processed_tic(self):
1961        """Return the processed total ion current of the mass spectrum."""
1962        return sum(self._processed_tic)
1963
1964    @property
1965    def tic(self):
1966        """Return the total ion current of the mass spectrum."""
1967        if self._processed_tic:
1968            return self._processed_tic
1969        else:
1970            return sum(self.abundance)
1971
1972    @property
1973    def mz_abun_tuples(self):
1974        """Return the m/z and abundance values of the mass spectrum as a list of tuples."""
1975        r = lambda x: int(round(x[0], 0), int(round(x[1], 0)))
1976
1977        return [r(i) for i in self]
1978
1979    @property
1980    def mz_abun_dict(self):
1981        """Return the m/z and abundance values of the mass spectrum as a dictionary."""
1982        r = lambda x: int(round(x, 0))
1983
1984        return {r(i[0]): r(i[1]) for i in self}