corems.ms_peak.factory.MSPeakClasses

  1__author__ = "Yuri E. Corilo"
  2__date__ = "Jun 12, 2019"
  3
  4import math
  5
  6from numpy import nan
  7
  8from corems.mass_spectra.calc import SignalProcessing as sp
  9from corems.molecular_formula.factory.MolecularFormulaFactory import MolecularFormula
 10from corems.ms_peak.calc.MSPeakCalc import MSPeakCalculation
 11
 12
 13class _MSPeak(MSPeakCalculation):
 14    """A class representing a peak in a mass spectrum.
 15
 16    Parameters:
 17    ----------
 18    polarity : int
 19        Acquisition polarity of the peak (``+1`` or ``-1``). Historically
 20        passed and stored as ``ion_charge``; that name is kept as a
 21        deprecated alias. True assignment charge lives on
 22        ``MolecularFormula.ion_charge`` after formula matching.
 23    mz_exp : float
 24        The experimental m/z value of the peak.
 25    abundance : float
 26        The abundance of the peak.
 27    resolving_power : float
 28        The resolving power of the peak.
 29    signal_to_noise : float
 30        The signal-to-noise ratio of the peak.
 31    indexes : list[int]
 32        The profile indexes of the peak.
 33    index : int
 34        The index of the peak.
 35    ms_parent : MSParent, optional
 36        The parent mass spectrum object.
 37    exp_freq : float, optional
 38        The experimental frequency of the peak.
 39
 40    Methods:
 41    -------
 42    * __len__().
 43        Returns the number of molecular formulas associated with the peak.
 44    * __setitem__(position, molecular_formula_obj).
 45        Sets the molecular formula at the specified position.
 46    * __getitem__(position) -> MolecularFormula.
 47        Returns the molecular formula at the specified position.
 48    * change_kendrick_base(kendrick_dict_base).
 49        Changes the kendrick base for the peak.
 50    * add_molecular_formula(molecular_formula_obj).
 51        Adds a molecular formula to the peak.
 52    * remove_molecular_formula(mf_obj).
 53        Removes a molecular formula from the peak.
 54    * clear_molecular_formulas().
 55        Clears all molecular formulas associated with the peak.
 56    * plot_simulation(sim_type="lorentz", ax=None, color="green", oversample_multiplier=1, delta_rp=0, mz_overlay=1).
 57        Plots the simulated peak.
 58    * plot(ax=None, color="black", derivative=True, deriv_color='red').
 59        Plots the peak.
 60    * best_molecular_formula_candidate().
 61        Returns the best molecular formula candidate for the peak.
 62    """
 63
 64    def __init__(
 65        self,
 66        polarity,
 67        mz_exp,
 68        abundance,
 69        resolving_power,
 70        signal_to_noise,
 71        indexes,
 72        index,
 73        ms_parent=None,
 74        exp_freq=None,
 75    ):
 76        self._ms_parent = ms_parent
 77        # Polarity only (±1). Assignment charge is MolecularFormula.ion_charge.
 78        self.polarity = int(polarity)
 79        self._mz_exp = float(mz_exp)
 80        self.mass = float(mz_exp) / float(self.polarity)
 81        self.abundance = float(abundance)
 82        self.resolving_power = float(resolving_power)
 83        self.signal_to_noise = float(signal_to_noise)
 84        # profile indexes
 85        self.peak_left_index = int(indexes[0])
 86        self.peak_apex_index = int(indexes[1])
 87        self.peak_right_index = int(indexes[2])
 88
 89        # mass spec obj index
 90        self.index = int(index)
 91        # parent mass spectrum obj instance
 92
 93        # updated after mass error prediction'
 94        self.predicted_std = None
 95        # updated after calibration'
 96        self.mz_cal = None
 97        # updated individual calculation'
 98        self.baseline_noise = None
 99
100        if exp_freq:
101            self.freq_exp = float(exp_freq)
102
103        if self._ms_parent is not None:
104            kendrick_dict_base = self._ms_parent.mspeaks_settings.kendrick_base
105        else:
106            kendrick_dict_base = {"C": 1, "H": 2}
107        self._kmd, self._kendrick_mass, self._nominal_km = self._calc_kmd(
108            kendrick_dict_base
109        )
110
111        "updated after molecular formula ID"
112
113        self.molecular_formulas = []
114        self._confidence_score = None
115        # placeholder for found isotopologues index
116        self.isotopologue_indexes = []
117        # placeholder for found isotopologues molecular formula obj
118        self.found_isotopologues = {}
119
120        # Label for what type of peak it is - real signal, noise, sinc wiggle, magnetron or harmonic peak, etc.
121        self.peak_type = None
122
123    @property
124    def ion_charge(self):
125        """Deprecated alias for :attr:`polarity` (±1).
126
127        Peaks do not carry assignment charge state. After molecular formula
128        matching, use ``MolecularFormula.ion_charge`` on assigned formulas.
129        """
130        return self.polarity
131
132    @ion_charge.setter
133    def ion_charge(self, value):
134        """Deprecated. Sets :attr:`polarity`."""
135        self.polarity = int(value)
136
137    def __len__(self):
138        return len(self.molecular_formulas)
139
140    def __setitem__(self, position, molecular_formula_obj):
141        self.molecular_formulas[position] = molecular_formula_obj
142
143    def __getitem__(self, position) -> MolecularFormula:
144        return self.molecular_formulas[position]
145
146    def change_kendrick_base(self, kendrick_dict_base):
147        """Changes the kendrick base for the peak.
148
149        Parameters:
150        ----------
151        kendrick_dict_base : dict
152            The kendrick base dictionary.
153            Default is {"C": 1, "H": 2}. (CH2)
154        """
155        self._kmd, self._kendrick_mass, self._nominal_km = self._calc_kmd(
156            kendrick_dict_base
157        )
158
159    def add_molecular_formula(self, molecular_formula_obj):
160        """Adds a molecular formula to the peak.
161
162        Parameters:
163        ----------
164        molecular_formula_obj : MolecularFormula
165            The molecular formula object to be added.
166
167        Returns:
168        -------
169        MolecularFormula
170            The molecular formula object added.
171
172        """
173        # freeze state
174        molecular_formula_obj._mspeak_parent = self
175
176        # new_mol_formula = deepcopy(molecular_formula_obj)
177        # add link mass spectrum obj instance
178
179        # new_mol_formula.mspeak_parent = self
180
181        self.molecular_formulas.append(molecular_formula_obj)
182
183        return molecular_formula_obj
184
185    def remove_molecular_formula(self, mf_obj):
186        """Removes a molecular formula from the peak.
187
188        Parameters:
189        ----------
190        mf_obj : MolecularFormula
191            The molecular formula object to be removed.
192        """
193        self.molecular_formulas.remove(mf_obj)
194
195    def clear_molecular_formulas(self):
196        """Clears all molecular formulas associated with the peak."""
197        self.molecular_formulas = []
198
199    @property
200    def mz_exp(self):
201        """The experimental m/z value of the peak."""
202        if self.mz_cal:
203            return self.mz_cal
204        else:
205            return self._mz_exp
206
207    @mz_exp.setter
208    def mz_exp(self, mz_exp):
209        """Sets the experimental m/z value of the peak."""
210        self._mz_exp = mz_exp
211
212    @property
213    def area(self):
214        """The area of the peak."""
215        if self._ms_parent.is_centroid:
216            return nan
217        else:
218            return self.calc_area()
219
220    @property
221    def nominal_mz_exp(self):
222        """The experimental nominal (integer) m/z value of the peak."""
223        return math.floor(self.mz_exp)
224
225    @property
226    def kmd(self):
227        """The Kendrick mass defect of the peak."""
228        return self._kmd
229
230    @property
231    def kendrick_mass(self):
232        """The Kendrick mass of the peak."""
233        return self._kendrick_mass
234
235    @property
236    def knm(self):
237        """The Kendrick nominal mass of the peak."""
238        return self._nominal_km
239
240    @property
241    def is_assigned(self) -> bool:
242        """Whether the peak is assigned or not."""
243        return bool(self.molecular_formulas)
244
245    def plot_simulation(
246        self,
247        sim_type="lorentz",
248        ax=None,
249        color="green",
250        oversample_multiplier=1,
251        delta_rp=0,
252        mz_overlay=1,
253    ):
254        """Plots the simulated peak.
255
256        Parameters:
257        ----------
258        sim_type : str, optional
259            The type of simulation to be plotted.
260            Default is "lorentz".
261        ax : matplotlib.axes, optional
262            The axes to plot the simulated peak.
263            Default is None.
264        color : str, optional
265            The color of the simulated peak.
266            Default is "green".
267        oversample_multiplier : int, optional
268            The oversample multiplier.
269            Default is 1.
270        delta_rp : int, optional
271            A delta value to the resolving power
272            Default is 0.
273        mz_overlay : int, optional
274            The mz overlay.
275            Default is 1.
276
277        Returns:
278        -------
279        matplotlib.axes
280            The axes where the simulated peak was plotted.
281
282        """
283        if self._ms_parent:
284            import matplotlib.pyplot as plt
285
286            x, y = eval(
287                "self."
288                + sim_type
289                + "(oversample_multiplier="
290                + str(oversample_multiplier)
291                + ", delta_rp="
292                + str(delta_rp)
293                + ", mz_overlay="
294                + str(mz_overlay)
295                + ")"
296            )
297
298            if ax is None:
299                ax = plt.gca()
300            ax.plot(x, y, color=color, label="Simulation")
301            ax.set(xlabel="m/z", ylabel="abundance")
302
303            plt.legend()
304            return ax
305
306    def plot(
307        self,
308        ax=None,
309        color: str = "black",
310        derivative: bool = True,
311        deriv_color: str = "red",
312    ):  # pragma: no cover
313        """Plots the peak.
314
315        Parameters:
316        ----------
317        ax : matplotlib.axes, optional
318            The axes to plot the peak.
319            Default is None.
320        color : str, optional
321            The color of the peak.
322            Default is "black".
323        derivative : bool, optional
324            Whether to plot the derivative of the peak.
325            Default is True.
326        deriv_color : str, optional
327            The color of the derivative of the peak.
328            Default is "red".
329
330        Returns:
331        -------
332        matplotlib.axes
333            The axes where the peak was plotted.
334
335        """
336        if self._ms_parent:
337            import matplotlib.pyplot as plt
338
339            if ax is None:
340                ax = plt.gca()
341            x = self._ms_parent.mz_exp_profile[
342                self.peak_left_index : self.peak_right_index
343            ]
344            y = self._ms_parent.abundance_profile[
345                self.peak_left_index : self.peak_right_index
346            ]
347
348            ax.plot(x, y, color=color, label="Data")
349            ax.set(xlabel="m/z", ylabel="abundance")
350            if derivative and not self._ms_parent.is_centroid:
351                dy = sp.derivate(
352                    self._ms_parent.abundance_profile[
353                        self.peak_left_index : self.peak_right_index + 1
354                    ]
355                )
356                ax.plot(x, dy, c=deriv_color)
357            else:
358                ax.plot(
359                    (self.mz_exp, self.mz_exp),
360                    (0, self.abundance),
361                    color=color,
362                    label="Data",
363                )
364
365            # plt.legend()
366
367            return ax
368
369        else:
370            raise AttributeError(
371                "No parent mass spectrum object found to plot the peak."
372            )
373
374    @property
375    def best_molecular_formula_candidate(self):
376        """The best molecular formula candidate for the peak.
377
378        Returns a single best formula candidate based on the user defined score method.
379        Score method is set with:
380            molecular_search_settings.score_method
381
382        Returns
383        -------
384        MolecularFormula
385            The best molecular formula candidate for the peak.
386
387        """
388        if (
389            self._ms_parent.molecular_search_settings.score_method
390            == "N_S_P_lowest_error"
391        ):
392            return self.cia_score_N_S_P_error()
393
394        elif (
395            self._ms_parent.molecular_search_settings.score_method == "S_P_lowest_error"
396        ):
397            return self.cia_score_S_P_error()
398
399        elif self._ms_parent.molecular_search_settings.score_method == "lowest_error":
400            return self.molecular_formula_lowest_error()
401
402        elif (
403            self._ms_parent.molecular_search_settings.score_method == "air_filter_error"
404        ):
405            return self.molecular_formula_air_filter()
406
407        elif (
408            self._ms_parent.molecular_search_settings.score_method
409            == "water_filter_error"
410        ):
411            return self.molecular_formula_water_filter()
412
413        elif (
414            self._ms_parent.molecular_search_settings.score_method
415            == "earth_filter_error"
416        ):
417            return self.molecular_formula_earth_filter()
418
419        elif self._ms_parent.molecular_search_settings.score_method == "prob_score":
420            return self.molecular_formula_highest_prob_score()
421        else:
422            raise TypeError(
423                "Unknown score method selected: % s, \
424                            Please check score_method at \
425                            encapsulation.settings.molecular_id.MolecularIDSettings.MolecularFormulaSearchSettings",
426                self._ms_parent.parameters.molecular_search.score_method,
427            )
428
429
430class ICRMassPeak(_MSPeak):
431    """A class representing a peak in an ICR mass spectrum."""
432
433    def __init__(self, *args, ms_parent=None, exp_freq=None):
434        super().__init__(*args, exp_freq=exp_freq, ms_parent=ms_parent)
435
436    def resolving_power_calc(
437        self, B, T, apodization_method: str = None
438    ):  # , kaiser_beta: float = None):
439        """Calculate the theoretical resolving power of the peak.
440
441        Note this presumes magntiude mode Fourier transform and fundamental frequency detection.
442        If you are doing aFT/eFT, or 2nd harmonic detection, you should adjust the B or T values accordingly.
443        Note that the window function corrections are an approximation and may not be 100% accurate.
444        Encourage you to use with caution and/or re-do the derivation independently and update....
445
446        Parameters
447        ----------
448        T: float
449            transient time
450        B: float
451            Magnetic Filed Strength (Tesla)
452        apodization_method: str, optional
453            The apodization method applied to the transient prior to FFT.
454            Default is None.
455
456        Returns
457        -------
458        float
459            Theoretical resolving power of the peak, with magnitude mode Fourier transform and fundamental frequency detection.
460
461
462        References
463        ----------
464        1. Marshall et al. (Mass Spectrom Rev. 1998 Jan-Feb;17(1):1-35.)
465            DOI: 10.1002/(SICI)1098-2787(1998)17:1<1::AID-MAS1>3.0.CO;2-K
466        """
467        # Calculate theoretical low pressure limit (undamped) ICR resolving power
468        # Peak polarity is ±1 (legacy field was named ion_charge).
469        RP_unwindowed = (1.274e7 * self.polarity * B * T) / (
470            self.mz_exp * self.polarity
471        )
472
473        # Default correction factor (no windowing)
474        CF = 1.0
475
476        # Apply correction factor based on the window function
477        if apodization_method == "Rectangular" or apodization_method is None:
478            CF = 1.0
479        elif apodization_method == "Hamming":
480            CF = 0.685
481        elif apodization_method == "Hanning":
482            CF = 0.618
483        elif apodization_method == "Blackman":
484            CF = 0.530
485        elif apodization_method == "Full-Sine":
486            CF = 0.736
487        elif (
488            apodization_method == "Half-Sine"
489            or apodization_method == "Kaiser"
490            or apodization_method == "Half-Kaiser"
491        ):
492            raise ValueError(
493                f"Apodization method {apodization_method} not supported yet"
494            )
495        else:
496            raise ValueError(f"Unknown window function: {apodization_method}")
497
498        """ # This was placeholder code, not sure these CF and equations are correct.
499        elif apodization_method == "Half-Sine":
500            CF = 0.85  # Approximate value
501        elif apodization_method == "Kaiser":
502            if kaiser_beta is not None:
503                CF = 1.3975 / (kaiser_beta + 0.5)
504            else:
505                raise ValueError("Beta parameter is required for Kaiser window")
506        elif apodization_method == "Half-Kaiser":
507            if kaiser_beta is not None:
508                k_half = 0.85  # Adjust this factor based on empirical data or literature
509                CF = (1.3975 / (kaiser_beta + 0.5)) * k_half
510            else:
511                raise ValueError("Beta parameter is required for Half-Kaiser window")
512        """
513
514        # Adjusted resolving power after applying the window function
515        RP_windowed = RP_unwindowed * CF
516        return RP_windowed
517
518    def set_calc_resolving_power(
519        self, B: float, T: float, apodization_method: str = None
520    ):
521        """Set the resolving power of the peak to the calculated one."""
522        self.resolving_power = self.resolving_power_calc(B, T, apodization_method)
523
524    def _mz_to_f_bruker(self):
525        """[Not Functional] Convert a peak m/z value to frequency
526
527        # Currently Broken - Not sure why
528        if self.mz_cal:
529            mz_val = self.mz_cal
530        else:
531            mz_val = self.mz_exp
532        Aterm, Bterm, Cterm = self._ms_parent.Aterm, self._ms_parent.Bterm, self._ms_parent.Cterm
533        # Check if the Bterm of Ledford equation scales with the ICR trap voltage or not then Bterm = Bterm*trap_voltage
534
535        if Cterm == 0:
536
537            if Bterm == 0:
538                #uncalibrated data
539                freq_domain = Aterm / mz_val
540
541            else:
542
543                freq_domain = (Aterm / (mz_val)) - Bterm
544
545        # @will I need you insight here, not sure what is the inverted ledford equation that Bruker refers to
546        else:
547
548            freq_domain = (Aterm / mz_val) + (Bterm / power(mz_val, 2)) + Cterm
549
550        return freq_domain
551        """
552        raise RuntimeError("Function not confirmed to work, disabled.")
553
554
555class TOFMassPeak(_MSPeak):
556    """A class representing a peak in a TOF mass spectrum."""
557
558    def __init__(self, *args, exp_freq=None):
559        super().__init__(*args, exp_freq=exp_freq)
560
561    def set_calc_resolving_power(self):
562        return 0
563
564
565class OrbiMassPeak(_MSPeak):
566    """A class representing a peak in an Orbitrap mass spectrum."""
567
568    def __init__(self, *args, exp_freq=None):
569        super().__init__(*args, exp_freq=exp_freq)
570
571    def set_calc_resolving_power(self):
572        return 0