corems.mass_spectrum.calc.PeakPicking

@author: Yuri E. Corilo @date: Jun 27, 2019

  1"""
  2@author: Yuri E. Corilo
  3@date: Jun 27, 2019
  4"""
  5
  6import warnings
  7from numpy import (
  8    hstack,
  9    inf,
 10    isnan,
 11    where,
 12    array,
 13    polyfit,
 14    nan,
 15    pad,
 16    zeros,
 17    searchsorted,
 18)
 19from corems.encapsulation.constant import Labels
 20from corems.mass_spectra.calc import SignalProcessing as sp
 21
 22
 23class PeakPicking:
 24    """Class for peak picking.
 25
 26    Parameters
 27    ----------
 28    None
 29
 30    Attributes
 31    ----------
 32    None
 33
 34    Methods
 35    -------
 36    * prepare_peak_picking_data().
 37        Prepare the mz, abundance, and frequence data for peak picking.
 38    * cut_mz_domain_peak_picking().
 39        Cut the m/z domain for peak picking.
 40    * extrapolate_axes_for_pp(mz=None, abund=None, freq=None).
 41        Extrapolate the m/z axis and fill the abundance axis with 0s.
 42    * do_peak_picking().
 43        Perform peak picking.
 44    * find_minima(apex_index, abundance, len_abundance, right=True).
 45        Find the minima of a peak.
 46    * linear_fit_calc(intes, massa, index_term, index_sign).
 47        Algebraic solution to a linear fit.
 48    * calculate_resolving_power(intes, massa, current_index).
 49        Calculate the resolving power of a peak.
 50    * cal_minima(mass, abun).
 51        Calculate the minima of a peak.
 52    * calc_centroid(mass, abund, freq).
 53        Calculate the centroid of a peak.
 54    * get_threshold(intes).
 55        Get the intensity threshold for peak picking.
 56    * algebraic_quadratic(list_mass, list_y).
 57        Find the apex of a peak - algebraically.
 58    * find_apex_fit_quadratic(mass, abund, freq, current_index).
 59        Find the apex of a peak.
 60    * check_prominence(abun, current_index, len_abundance, peak_height_diff).
 61        Check the prominence of a peak.
 62    * use_the_max(mass, abund, current_index, len_abundance, peak_height_diff).
 63        Use the max peak height as the centroid.
 64    * calc_centroid_legacy(mass, abund, freq).
 65        Legacy centroid calculation. Deprecated - for deletion.
 66
 67    """
 68
 69    def prepare_peak_picking_data(self):
 70        """Prepare the data for peak picking.
 71
 72        This function will prepare the m/z, abundance, and frequency data for peak picking according to the settings.
 73
 74        Returns
 75        -------
 76        mz : ndarray
 77            The m/z axis.
 78        abundance : ndarray
 79            The abundance axis.
 80        freq : ndarray or None
 81            The frequency axis, if available.
 82        """
 83        # First apply cut_mz_domain_peak_picking
 84        mz, abundance, freq = self.cut_mz_domain_peak_picking()
 85
 86        # Then extrapolate the axes for peak picking
 87        if self.settings.picking_point_extrapolate > 0:
 88            mz, abundance, freq = self.extrapolate_axes_for_pp(mz, abundance, freq)
 89        return mz, abundance, freq
 90
 91    def cut_mz_domain_peak_picking(self):
 92        """
 93        Cut the m/z domain for peak picking.
 94
 95        Simplified function
 96
 97        Returns
 98        -------
 99        mz_domain_X_low_cutoff : ndarray
100            The m/z values within the specified range.
101        mz_domain_low_Y_cutoff : ndarray
102            The abundance values within the specified range.
103        freq_domain_low_Y_cutoff : ndarray or None
104            The frequency values within the specified range, if available.
105
106        """
107        max_picking_mz = self.settings.max_picking_mz
108        min_picking_mz = self.settings.min_picking_mz
109
110        # min_start =  where(self.mz_exp_profile  > min_picking_mz)[0][0]
111        # max_final =  where(self.mz_exp_profile < max_picking_mz)[-1][-1]
112        min_start = searchsorted(a=self.mz_exp_profile, v=min_picking_mz)
113        max_final = searchsorted(a=self.mz_exp_profile, v=max_picking_mz)
114
115        if self.has_frequency:
116            if self.freq_exp_profile.any():
117                return (
118                    self.mz_exp_profile[min_start:max_final],
119                    self.abundance_profile[min_start:max_final],
120                    self.freq_exp_profile[min_start:max_final],
121                )
122
123        else:
124            return (
125                self.mz_exp_profile[min_start:max_final],
126                self.abundance_profile[min_start:max_final],
127                None,
128            )
129
130    def legacy_cut_mz_domain_peak_picking(self):
131        """
132        Cut the m/z domain for peak picking.
133        DEPRECATED
134        Returns
135        -------
136        mz_domain_X_low_cutoff : ndarray
137            The m/z values within the specified range.
138        mz_domain_low_Y_cutoff : ndarray
139            The abundance values within the specified range.
140        freq_domain_low_Y_cutoff : ndarray or None
141            The frequency values within the specified range, if available.
142
143        """
144        max_picking_mz = self.settings.max_picking_mz
145        min_picking_mz = self.settings.min_picking_mz
146
147        min_final = where(self.mz_exp_profile > min_picking_mz)[-1][-1]
148        min_start = where(self.mz_exp_profile > min_picking_mz)[0][0]
149
150        (
151            mz_domain_X_low_cutoff,
152            mz_domain_low_Y_cutoff,
153        ) = (
154            self.mz_exp_profile[min_start:min_final],
155            self.abundance_profile[min_start:min_final],
156        )
157
158        max_final = where(self.mz_exp_profile < max_picking_mz)[-1][-1]
159        max_start = where(self.mz_exp_profile < max_picking_mz)[0][0]
160
161        if self.has_frequency:
162            if self.freq_exp_profile.any():
163                freq_domain_low_Y_cutoff = self.freq_exp_profile[min_start:min_final]
164
165                return (
166                    mz_domain_X_low_cutoff[max_start:max_final],
167                    mz_domain_low_Y_cutoff[max_start:max_final],
168                    freq_domain_low_Y_cutoff[max_start:max_final],
169                )
170
171        else:
172            return (
173                mz_domain_X_low_cutoff[max_start:max_final],
174                mz_domain_low_Y_cutoff[max_start:max_final],
175                None,
176            )
177
178    @staticmethod
179    def extrapolate_axis(initial_array, pts):
180        """
181        This function will extrapolate an input array in both directions by N pts.
182
183        Parameters
184        ----------
185        initial_array : ndarray
186            The input array.
187        pts : int
188            The number of points to extrapolate.
189
190        Returns
191        -------
192        ndarray
193            The extrapolated array.
194
195        Notes
196        --------
197        This is a static method.
198        """
199        initial_array_len = len(initial_array)
200        right_delta = initial_array[-1] - initial_array[-2]
201        left_delta = initial_array[1] - initial_array[0]
202
203        # Create an array with extra space for extrapolation
204        pad_array = zeros(initial_array_len + 2 * pts)
205
206        # Copy original array into the middle of the padded array
207        pad_array[pts : pts + initial_array_len] = initial_array
208
209        # Extrapolate the right side
210        for pt in range(pts):
211            final_value = initial_array[-1]
212            value_to_add = right_delta * (pt + 1)
213            new_value = final_value + value_to_add
214            pad_array[initial_array_len + pts + pt] = new_value
215
216        # Extrapolate the left side
217        for pt in range(pts):
218            first_value = initial_array[0]
219            value_to_subtract = left_delta * (pt + 1)
220            new_value = first_value - value_to_subtract
221            pad_array[pts - pt - 1] = new_value
222
223        return pad_array
224
225    def extrapolate_axes_for_pp(self, mz=None, abund=None, freq=None):
226        """Extrapolate the m/z axis and fill the abundance axis with 0s.
227
228        Parameters
229        ----------
230        mz : ndarray or None
231            The m/z axis, if available. If None, the experimental m/z axis is used.
232        abund : ndarray or None
233            The abundance axis, if available. If None, the experimental abundance axis is used.
234        freq : ndarray or None
235            The frequency axis, if available. If None, the experimental frequency axis is used.
236
237        Returns
238        -------
239        mz : ndarray
240            The extrapolated m/z axis.
241        abund : ndarray
242            The abundance axis with 0s filled.
243        freq : ndarray or None
244            The extrapolated frequency axis, if available.
245
246        Notes
247        --------
248        This function will extrapolate the mz axis by the number of datapoints specified in the settings,
249        and fill the abundance axis with 0s.
250        This should prevent peak picking issues at the spectrum edge.
251
252        """
253        # Check if the input arrays are provided
254        if mz is None or abund is None:
255            mz, abund = self.mz_exp_profile, self.abundance_profile
256            if self.has_frequency:
257                freq = self.freq_exp_profile
258            else:
259                freq = None
260        pts = self.settings.picking_point_extrapolate
261        if pts == 0:
262            return mz, abund, freq
263
264        mz = self.extrapolate_axis(mz, pts)
265        abund = pad(abund, (pts, pts), mode="constant", constant_values=(0, 0))
266        if freq is not None:
267            freq = self.extrapolate_axis(freq, pts)
268        return mz, abund, freq
269
270    def do_peak_picking(self):
271        """Perform peak picking."""
272        mz, abundance, freq = self.prepare_peak_picking_data()
273
274        if (
275            self.label == Labels.bruker_frequency
276            or self.label == Labels.midas_frequency
277        ):
278            self.calc_centroid(mz, abundance, freq)
279
280        elif self.label == Labels.thermo_profile:
281            self.calc_centroid(mz, abundance, freq)
282
283        elif self.label == Labels.bruker_profile:
284            self.calc_centroid(mz, abundance, freq)
285
286        elif self.label == Labels.booster_profile:
287            self.calc_centroid(mz, abundance, freq)
288
289        elif self.label == Labels.simulated_profile:
290            self.calc_centroid(mz, abundance, freq)
291
292        else:
293            raise Exception("Unknow mass spectrum type", self.label)
294
295    def find_minima(self, apex_index, abundance, len_abundance, right=True):
296        """Find the minima of a peak.
297
298        Parameters
299        ----------
300        apex_index : int
301            The index of the peak apex.
302        abundance : ndarray
303            The abundance values.
304        len_abundance : int
305            The length of the abundance array.
306        right : bool, optional
307            Flag indicating whether to search for minima to the right of the apex (default is True).
308
309        Returns
310        -------
311        int
312            The index of the minima.
313
314        """
315        j = apex_index
316
317        if right:
318            minima = abundance[j] > abundance[j + 1]
319        else:
320            minima = abundance[j] > abundance[j - 1]
321
322        while minima:
323            if j == 1 or j == len_abundance - 2:
324                break
325
326            if right:
327                j += 1
328
329                minima = abundance[j] >= abundance[j + 1]
330
331            else:
332                j -= 1
333                minima = abundance[j] >= abundance[j - 1]
334
335        if right:
336            return j
337        else:
338            return j
339
340    @staticmethod
341    def linear_fit_calc(intes, massa, index_term, index_sign):
342        """
343        Algebraic solution to a linear fit - roughly 25-50x faster than numpy polyfit when passing only two vals and doing a 1st order fit
344
345        Parameters
346        ----------
347        intes : ndarray
348            The intensity values.
349        massa : ndarray
350            The mass values.
351        index_term : int
352            The index of the current term.
353        index_sign : str
354            The index sign
355
356        Returns
357        -------
358        ndarray
359            The coefficients of the linear fit.
360
361        Notes
362        --------
363        This is a static method.
364        """
365        if index_sign == "+":
366            x1, x2 = massa[index_term], massa[index_term + 1]
367            y1, y2 = intes[index_term], intes[index_term + 1]
368        elif index_sign == "-":
369            x1, x2 = massa[index_term], massa[index_term - 1]
370            y1, y2 = intes[index_term], intes[index_term - 1]
371        else:
372            warnings.warn("error in linear fit calc, unknown index sign")
373
374        # Calculate the slope (m)
375        slope = (y2 - y1) / (x2 - x1)
376
377        # Calculate the intercept (b)
378        intercept = y1 - slope * x1
379
380        # The coefficients array would be [slope, intercept]
381        coefficients = array([slope, intercept])
382        return coefficients
383
384    def calculate_resolving_power(self, intes, massa, current_index):
385        """Calculate the resolving power of a peak.
386
387        Parameters
388        ----------
389        intes : ndarray
390            The intensity values.
391        massa : ndarray
392            The mass values.
393        current_index : int
394            The index of the current peak.
395
396        Returns
397        -------
398        float
399            The resolving power of the peak, or ``nan`` if a half-maximum
400            crossing does not exist on both sides (spectrum edge truncates
401            the peak).
402
403        Notes
404        --------
405        This is a conservative calculation of resolving power,
406        the peak need to be resolved at least at the half-maximum magnitude,
407        otherwise, the combined full width at half maximum is used to calculate resolving power.
408
409        """
410
411        peak_height = intes[current_index]
412        target_peak_height = peak_height / 2
413        n_points = len(intes)
414
415        index_minus = current_index
416        while index_minus > 0 and intes[index_minus] >= target_peak_height:
417            index_minus -= 1
418        if intes[index_minus] >= target_peak_height:
419            return nan
420
421        index_plus = current_index
422        while index_plus < n_points - 1 and intes[index_plus] >= target_peak_height:
423            index_plus += 1
424        if intes[index_plus] >= target_peak_height:
425            return nan
426
427        if self.mspeaks_settings.legacy_centroid_polyfit:
428            x = [massa[index_minus], massa[index_minus + 1]]
429            y = [intes[index_minus], intes[index_minus + 1]]
430            coefficients = polyfit(x, y, 1)
431        else:
432            coefficients = self.linear_fit_calc(
433                intes, massa, index_minus, index_sign="+"
434            )
435
436        a = coefficients[0]
437        b = coefficients[1]
438        if self.mspeaks_settings.legacy_resolving_power:
439            y_intercept = intes[index_minus] + (
440                (intes[index_minus + 1] - intes[index_minus]) / 2
441            )
442        else:
443            y_intercept = target_peak_height
444        massa1 = (y_intercept - b) / a
445
446        if self.mspeaks_settings.legacy_centroid_polyfit:
447            x = [massa[index_plus], massa[index_plus - 1]]
448            y = [intes[index_plus], intes[index_plus - 1]]
449            coefficients = polyfit(x, y, 1)
450        else:
451            coefficients = self.linear_fit_calc(
452                intes, massa, index_plus, index_sign="-"
453            )
454
455        a = coefficients[0]
456        b = coefficients[1]
457
458        if self.mspeaks_settings.legacy_resolving_power:
459            y_intercept = intes[index_plus - 1] + (
460                (intes[index_plus] - intes[index_plus - 1]) / 2
461            )
462        else:
463            y_intercept = target_peak_height
464
465        massa2 = (y_intercept - b) / a
466
467        if massa1 > massa2:
468            resolvingpower = massa[current_index] / (massa1 - massa2)
469
470        else:
471            resolvingpower = massa[current_index] / (massa2 - massa1)
472
473        return resolvingpower
474
475    def cal_minima(self, mass, abun):
476        """Calculate the minima of a peak.
477
478        Parameters
479        ----------
480        mass : ndarray
481            The mass values.
482        abun : ndarray
483            The abundance values.
484
485        Returns
486        -------
487        ndarray or None
488            The mass values at the minima, if found.
489
490        """
491        abun = -abun
492
493        dy = abun[1:] - abun[:-1]
494
495        # replaces nan for infinity
496        indices_nan = where(isnan(abun))[0]
497
498        if indices_nan.size:
499            abun[indices_nan] = inf
500            dy[where(isnan(dy))[0]] = inf
501
502        indexes = where((hstack((dy, 0)) < 0) & (hstack((0, dy)) > 0))[0]
503
504        if indexes.size:
505            return mass[indexes], abun[indexes]
506
507    def calc_centroid(self, mass, abund, freq):
508        """Calculate the centroid of a peak.
509
510        Parameters
511        ----------
512        mass : ndarray
513            The mass values.
514        abund : ndarray
515            The abundance values.
516        freq : ndarray or None
517            The frequency values, if available.
518
519        Returns
520        -------
521        None
522
523        """
524
525        max_height = self.mspeaks_settings.peak_height_max_percent
526        max_prominence = self.mspeaks_settings.peak_max_prominence_percent
527        min_peak_datapoints = self.mspeaks_settings.min_peak_datapoints
528        peak_derivative_threshold = self.mspeaks_settings.peak_derivative_threshold
529        max_abun = max(abund)
530        peak_height_diff = lambda hi, li: ((abund[hi] - abund[li]) / max_abun) * 100
531
532        domain = mass
533        signal = abund
534        len_signal = len(signal)
535
536        signal_threshold, factor = self.get_threshold(abund)
537        max_signal = factor
538
539        correct_baseline = False
540
541        include_indexes = sp.peak_picking_first_derivative(
542            domain,
543            signal,
544            max_height,
545            max_prominence,
546            max_signal,
547            min_peak_datapoints,
548            peak_derivative_threshold,
549            signal_threshold=signal_threshold,
550            correct_baseline=correct_baseline,
551            abun_norm=1,
552            plot_res=False,
553        )
554
555        for indexes_tuple in include_indexes:
556            apex_index = indexes_tuple[1]
557
558            peak_indexes = self.check_prominence(
559                abund, apex_index, len_signal, peak_height_diff
560            )
561
562            if peak_indexes:
563                mz_exp_centroid, freq_centr, intes_centr = self.find_apex_fit_quadratic(
564                    mass, abund, freq, apex_index
565                )
566
567                if mz_exp_centroid:
568                    peak_resolving_power = self.calculate_resolving_power(
569                        abund, mass, apex_index
570                    )
571                    s2n = intes_centr / self.baseline_noise_std
572                    self.add_mspeak(
573                        self.polarity,
574                        mz_exp_centroid,
575                        abund[apex_index],
576                        peak_resolving_power,
577                        s2n,
578                        indexes_tuple,
579                        exp_freq=freq_centr,
580                        ms_parent=self,
581                    )
582                # pyplot.plot(domain[start_index: final_index + 1], signal[start_index:final_index + 1], c='black')
583                # pyplot.show()
584
585    def get_threshold(self, intes):
586        """Get the intensity threshold for peak picking.
587
588        Parameters
589        ----------
590        intes : ndarray
591            The intensity values.
592
593        Returns
594        -------
595        float
596            The intensity threshold.
597        float
598            The factor to multiply the intensity threshold by.
599        """
600
601        intes = array(intes).astype(float)
602
603        self._resolve_noise_threshold_method()
604        noise_threshold_method = self.settings.noise_threshold_method
605
606        if noise_threshold_method == "minima":
607            if self.is_centroid:
608                warnings.warn(
609                    "Auto threshould is disabled for centroid data, returning 0"
610                )
611                factor = 1
612                abundance_threshold = 1e-20
613            # print(self.settings.noise_threshold_min_std)
614            else:
615                abundance_threshold = self.baseline_noise + (
616                    self.settings.noise_threshold_min_std * self.baseline_noise_std
617                )
618                factor = 1
619
620        elif noise_threshold_method == "signal_noise":
621            abundance_threshold = self.settings.noise_threshold_min_s2n
622            if self.is_centroid:
623                factor = 1
624            else:
625                factor = self.baseline_noise_std
626
627        elif noise_threshold_method == "relative_abundance":
628            abundance_threshold = self.settings.noise_threshold_min_relative_abundance
629            factor = intes.max() / 100
630
631        elif noise_threshold_method == "absolute_abundance":
632            abundance_threshold = self.settings.noise_threshold_absolute_abundance
633            factor = 1
634
635        elif noise_threshold_method == "log":
636            if self.is_centroid:
637                raise Exception("log noise Not tested for centroid data")
638            abundance_threshold = self.settings.noise_threshold_log_nsigma
639            factor = self.baseline_noise_std
640
641        else:
642            raise Exception(
643                "%s method was not implemented, please refer to corems.mass_spectrum.calc.NoiseCalc Class"
644                % noise_threshold_method
645            )
646
647        return abundance_threshold, factor
648
649    @staticmethod
650    def algebraic_quadratic(list_mass, list_y):
651        """
652        Find the apex of a peak - algebraically.
653        Faster than using numpy polyfit by ~28x per fit.
654
655        Parameters
656        ----------
657        list_mass : ndarray
658            list of m/z values (3 points)
659        list_y : ndarray
660            list of abundance values (3 points)
661
662        Returns
663        -------
664        a, b, c: float
665            coefficients of the quadratic equation.
666
667        Notes
668        --------
669        This is a static method.
670        """
671        x_1, x_2, x_3 = list_mass
672        y_1, y_2, y_3 = list_y
673
674        a = (
675            y_1 / ((x_1 - x_2) * (x_1 - x_3))
676            + y_2 / ((x_2 - x_1) * (x_2 - x_3))
677            + y_3 / ((x_3 - x_1) * (x_3 - x_2))
678        )
679
680        b = (
681            -y_1 * (x_2 + x_3) / ((x_1 - x_2) * (x_1 - x_3))
682            - y_2 * (x_1 + x_3) / ((x_2 - x_1) * (x_2 - x_3))
683            - y_3 * (x_1 + x_2) / ((x_3 - x_1) * (x_3 - x_2))
684        )
685
686        c = (
687            y_1 * x_2 * x_3 / ((x_1 - x_2) * (x_1 - x_3))
688            + y_2 * x_1 * x_3 / ((x_2 - x_1) * (x_2 - x_3))
689            + y_3 * x_1 * x_2 / ((x_3 - x_1) * (x_3 - x_2))
690        )
691        return a, b, c
692
693    def find_apex_fit_quadratic(self, mass, abund, freq, current_index):
694        """
695        Find the apex of a peak.
696
697        Parameters
698        ----------
699        mass : ndarray
700            The mass values.
701        abund : ndarray
702            The abundance values.
703        freq : ndarray or None
704            The frequency values, if available.
705        current_index : int
706            The index of the current peak.
707
708
709        Returns
710        -------
711        float
712            The m/z value of the peak apex.
713        float
714            The frequency value of the peak apex, if available.
715        float
716            The abundance value of the peak apex.
717
718        """
719        # calc prominence
720        # peak_indexes = self.check_prominence(abund, current_index, len_abundance, peak_height_diff )
721
722        # if not peak_indexes:
723
724        #    return None, None, None, None
725
726        # else:
727
728        # fit parabola to three most abundant datapoints
729        list_mass = [
730            mass[current_index - 1],
731            mass[current_index],
732            mass[current_index + 1],
733        ]
734        list_y = [
735            abund[current_index - 1],
736            abund[current_index],
737            abund[current_index + 1],
738        ]
739
740        if self.mspeaks_settings.legacy_centroid_polyfit:
741            z = polyfit(list_mass, list_y, 2)
742            a = z[0]
743            b = z[1]
744        else:
745            a, b, c = self.algebraic_quadratic(list_mass, list_y)
746
747        calculated = -b / (2 * a)
748
749        if calculated < 1 or int(calculated) != int(list_mass[1]):
750            mz_exp_centroid = list_mass[1]
751
752        else:
753            mz_exp_centroid = calculated
754
755        if (
756            self.label == Labels.bruker_frequency
757            or self.label == Labels.midas_frequency
758        ):
759            # fit parabola to three most abundant frequency datapoints
760            list_freq = [
761                freq[current_index - 1],
762                freq[current_index],
763                freq[current_index + 1],
764            ]
765            if self.mspeaks_settings.legacy_centroid_polyfit:
766                z = polyfit(list_mass, list_y, 2)
767                a = z[0]
768                b = z[1]
769            else:
770                a, b, c = self.algebraic_quadratic(list_mass, list_y)
771
772            calculated_freq = -b / (2 * a)
773
774            if calculated_freq < 1 or int(calculated_freq) != freq[current_index]:
775                freq_centr = list_freq[1]
776
777            else:
778                freq_centr = calculated_freq
779
780        else:
781            freq_centr = None
782
783        if self.mspeaks_settings.legacy_centroid_polyfit:
784            abundance_centroid = abund[current_index]
785        else:
786            abundance_centroid = a * mz_exp_centroid**2 + b * mz_exp_centroid + c
787
788        return mz_exp_centroid, freq_centr, abundance_centroid
789
790    def check_prominence(
791        self, abun, current_index, len_abundance, peak_height_diff
792    ) -> tuple or False:
793        """Check the prominence of a peak.
794
795        Parameters
796        ----------
797        abun : ndarray
798            The abundance values.
799        current_index : int
800            The index of the current peak.
801        len_abundance : int
802            The length of the abundance array.
803        peak_height_diff : function
804            The function to calculate the peak height difference.
805
806        Returns
807        -------
808        tuple or False
809            A tuple containing the indexes of the peak, if the prominence is above the threshold.
810            Otherwise, False.
811
812        """
813
814        final_index = self.find_minima(current_index, abun, len_abundance, right=True)
815
816        start_index = self.find_minima(current_index, abun, len_abundance, right=False)
817
818        peak_indexes = (current_index - 1, current_index, current_index + 1)
819
820        if (
821            min(
822                peak_height_diff(current_index, start_index),
823                peak_height_diff(current_index, final_index),
824            )
825            > self.mspeaks_settings.peak_min_prominence_percent
826        ):
827            return peak_indexes
828
829        else:
830            return False
831
832    def use_the_max(self, mass, abund, current_index, len_abundance, peak_height_diff):
833        """Use the max peak height as the centroid
834
835        Parameters
836        ----------
837        mass : ndarray
838            The mass values.
839        abund : ndarray
840            The abundance values.
841        current_index : int
842            The index of the current peak.
843        len_abundance : int
844            The length of the abundance array.
845        peak_height_diff : function
846            The function to calculate the peak height difference.
847
848        Returns
849        -------
850        float
851            The m/z value of the peak apex.
852        float
853            The abundance value of the peak apex.
854        tuple or None
855            A tuple containing the indexes of the peak, if the prominence is above the threshold.
856            Otherwise, None.
857        """
858
859        peak_indexes = self.check_prominence(
860            abund, current_index, len_abundance, peak_height_diff
861        )
862
863        if not peak_indexes:
864            return None, None, None
865
866        else:
867            return mass[current_index], abund[current_index], peak_indexes
868
869    def calc_centroid_legacy(self, mass, abund, freq):
870        """Legacy centroid calculation
871        Deprecated - for deletion.
872
873        """
874        warnings.warn(
875            "Legacy centroid calculation is deprecated. Please use the new centroid calculation method."
876        )
877        pass
878        if False:
879            len_abundance = len(abund)
880
881            max_abundance = max(abund)
882
883            peak_height_diff = lambda hi, li: (
884                ((abund[hi] - abund[li]) / max_abundance) * 100
885            )
886
887            abundance_threshold, factor = self.get_threshold(abund)
888            # print(abundance_threshold, factor)
889            # find indices of all peaks
890            dy = abund[1:] - abund[:-1]
891
892            # replaces nan for infi nity
893            indices_nan = where(isnan(abund))[0]
894
895            if indices_nan.size:
896                abund[indices_nan] = inf
897                dy[where(isnan(dy))[0]] = inf
898
899            indexes = where((hstack((dy, 0)) < 0) & (hstack((0, dy)) > 0))[0]
900
901            # noise threshold
902            if indexes.size and abundance_threshold is not None:
903                indexes = indexes[abund[indexes] / factor >= abundance_threshold]
904            # filter out 'peaks' within 3 points of the spectrum limits
905            # remove entries within 3 points of upper limit
906            indexes = [x for x in indexes if (len_abundance - x) > 3]
907            # remove entries within 3 points of zero
908            indexes = [x for x in indexes if x > 3]
909
910            for current_index in indexes:
911                if self.label == Labels.simulated_profile:
912                    mz_exp_centroid, intes_centr, peak_indexes = self.use_the_max(
913                        mass, abund, current_index, len_abundance, peak_height_diff
914                    )
915                    if mz_exp_centroid:
916                        peak_resolving_power = self.calculate_resolving_power(
917                            abund, mass, current_index
918                        )
919                        s2n = intes_centr / self.baseline_noise_std
920                        freq_centr = None
921                        self.add_mspeak(
922                            self.polarity,
923                            mz_exp_centroid,
924                            abund[current_index],
925                            peak_resolving_power,
926                            s2n,
927                            peak_indexes,
928                            exp_freq=freq_centr,
929                            ms_parent=self,
930                        )
931
932                else:
933                    mz_exp_centroid, freq_centr, intes_centr, peak_indexes = (
934                        self.find_apex_fit_quadratic(
935                            mass,
936                            abund,
937                            freq,
938                            current_index,
939                            len_abundance,
940                            peak_height_diff,
941                        )
942                    )
943                    if mz_exp_centroid:
944                        try:
945                            peak_resolving_power = self.calculate_resolving_power(
946                                abund, mass, current_index
947                            )
948                        except IndexError:
949                            print("index error, skipping peak")
950                            continue
951
952                        s2n = intes_centr / self.baseline_noise_std
953                        self.add_mspeak(
954                            self.polarity,
955                            mz_exp_centroid,
956                            abund[current_index],
957                            peak_resolving_power,
958                            s2n,
959                            peak_indexes,
960                            exp_freq=freq_centr,
961                            ms_parent=self,
962                        )
class PeakPicking: 
 24class PeakPicking:
 25    """Class for peak picking.
 26
 27    Parameters
 28    ----------
 29    None
 30
 31    Attributes
 32    ----------
 33    None
 34
 35    Methods
 36    -------
 37    * prepare_peak_picking_data().
 38        Prepare the mz, abundance, and frequence data for peak picking.
 39    * cut_mz_domain_peak_picking().
 40        Cut the m/z domain for peak picking.
 41    * extrapolate_axes_for_pp(mz=None, abund=None, freq=None).
 42        Extrapolate the m/z axis and fill the abundance axis with 0s.
 43    * do_peak_picking().
 44        Perform peak picking.
 45    * find_minima(apex_index, abundance, len_abundance, right=True).
 46        Find the minima of a peak.
 47    * linear_fit_calc(intes, massa, index_term, index_sign).
 48        Algebraic solution to a linear fit.
 49    * calculate_resolving_power(intes, massa, current_index).
 50        Calculate the resolving power of a peak.
 51    * cal_minima(mass, abun).
 52        Calculate the minima of a peak.
 53    * calc_centroid(mass, abund, freq).
 54        Calculate the centroid of a peak.
 55    * get_threshold(intes).
 56        Get the intensity threshold for peak picking.
 57    * algebraic_quadratic(list_mass, list_y).
 58        Find the apex of a peak - algebraically.
 59    * find_apex_fit_quadratic(mass, abund, freq, current_index).
 60        Find the apex of a peak.
 61    * check_prominence(abun, current_index, len_abundance, peak_height_diff).
 62        Check the prominence of a peak.
 63    * use_the_max(mass, abund, current_index, len_abundance, peak_height_diff).
 64        Use the max peak height as the centroid.
 65    * calc_centroid_legacy(mass, abund, freq).
 66        Legacy centroid calculation. Deprecated - for deletion.
 67
 68    """
 69
 70    def prepare_peak_picking_data(self):
 71        """Prepare the data for peak picking.
 72
 73        This function will prepare the m/z, abundance, and frequency data for peak picking according to the settings.
 74
 75        Returns
 76        -------
 77        mz : ndarray
 78            The m/z axis.
 79        abundance : ndarray
 80            The abundance axis.
 81        freq : ndarray or None
 82            The frequency axis, if available.
 83        """
 84        # First apply cut_mz_domain_peak_picking
 85        mz, abundance, freq = self.cut_mz_domain_peak_picking()
 86
 87        # Then extrapolate the axes for peak picking
 88        if self.settings.picking_point_extrapolate > 0:
 89            mz, abundance, freq = self.extrapolate_axes_for_pp(mz, abundance, freq)
 90        return mz, abundance, freq
 91
 92    def cut_mz_domain_peak_picking(self):
 93        """
 94        Cut the m/z domain for peak picking.
 95
 96        Simplified function
 97
 98        Returns
 99        -------
100        mz_domain_X_low_cutoff : ndarray
101            The m/z values within the specified range.
102        mz_domain_low_Y_cutoff : ndarray
103            The abundance values within the specified range.
104        freq_domain_low_Y_cutoff : ndarray or None
105            The frequency values within the specified range, if available.
106
107        """
108        max_picking_mz = self.settings.max_picking_mz
109        min_picking_mz = self.settings.min_picking_mz
110
111        # min_start =  where(self.mz_exp_profile  > min_picking_mz)[0][0]
112        # max_final =  where(self.mz_exp_profile < max_picking_mz)[-1][-1]
113        min_start = searchsorted(a=self.mz_exp_profile, v=min_picking_mz)
114        max_final = searchsorted(a=self.mz_exp_profile, v=max_picking_mz)
115
116        if self.has_frequency:
117            if self.freq_exp_profile.any():
118                return (
119                    self.mz_exp_profile[min_start:max_final],
120                    self.abundance_profile[min_start:max_final],
121                    self.freq_exp_profile[min_start:max_final],
122                )
123
124        else:
125            return (
126                self.mz_exp_profile[min_start:max_final],
127                self.abundance_profile[min_start:max_final],
128                None,
129            )
130
131    def legacy_cut_mz_domain_peak_picking(self):
132        """
133        Cut the m/z domain for peak picking.
134        DEPRECATED
135        Returns
136        -------
137        mz_domain_X_low_cutoff : ndarray
138            The m/z values within the specified range.
139        mz_domain_low_Y_cutoff : ndarray
140            The abundance values within the specified range.
141        freq_domain_low_Y_cutoff : ndarray or None
142            The frequency values within the specified range, if available.
143
144        """
145        max_picking_mz = self.settings.max_picking_mz
146        min_picking_mz = self.settings.min_picking_mz
147
148        min_final = where(self.mz_exp_profile > min_picking_mz)[-1][-1]
149        min_start = where(self.mz_exp_profile > min_picking_mz)[0][0]
150
151        (
152            mz_domain_X_low_cutoff,
153            mz_domain_low_Y_cutoff,
154        ) = (
155            self.mz_exp_profile[min_start:min_final],
156            self.abundance_profile[min_start:min_final],
157        )
158
159        max_final = where(self.mz_exp_profile < max_picking_mz)[-1][-1]
160        max_start = where(self.mz_exp_profile < max_picking_mz)[0][0]
161
162        if self.has_frequency:
163            if self.freq_exp_profile.any():
164                freq_domain_low_Y_cutoff = self.freq_exp_profile[min_start:min_final]
165
166                return (
167                    mz_domain_X_low_cutoff[max_start:max_final],
168                    mz_domain_low_Y_cutoff[max_start:max_final],
169                    freq_domain_low_Y_cutoff[max_start:max_final],
170                )
171
172        else:
173            return (
174                mz_domain_X_low_cutoff[max_start:max_final],
175                mz_domain_low_Y_cutoff[max_start:max_final],
176                None,
177            )
178
179    @staticmethod
180    def extrapolate_axis(initial_array, pts):
181        """
182        This function will extrapolate an input array in both directions by N pts.
183
184        Parameters
185        ----------
186        initial_array : ndarray
187            The input array.
188        pts : int
189            The number of points to extrapolate.
190
191        Returns
192        -------
193        ndarray
194            The extrapolated array.
195
196        Notes
197        --------
198        This is a static method.
199        """
200        initial_array_len = len(initial_array)
201        right_delta = initial_array[-1] - initial_array[-2]
202        left_delta = initial_array[1] - initial_array[0]
203
204        # Create an array with extra space for extrapolation
205        pad_array = zeros(initial_array_len + 2 * pts)
206
207        # Copy original array into the middle of the padded array
208        pad_array[pts : pts + initial_array_len] = initial_array
209
210        # Extrapolate the right side
211        for pt in range(pts):
212            final_value = initial_array[-1]
213            value_to_add = right_delta * (pt + 1)
214            new_value = final_value + value_to_add
215            pad_array[initial_array_len + pts + pt] = new_value
216
217        # Extrapolate the left side
218        for pt in range(pts):
219            first_value = initial_array[0]
220            value_to_subtract = left_delta * (pt + 1)
221            new_value = first_value - value_to_subtract
222            pad_array[pts - pt - 1] = new_value
223
224        return pad_array
225
226    def extrapolate_axes_for_pp(self, mz=None, abund=None, freq=None):
227        """Extrapolate the m/z axis and fill the abundance axis with 0s.
228
229        Parameters
230        ----------
231        mz : ndarray or None
232            The m/z axis, if available. If None, the experimental m/z axis is used.
233        abund : ndarray or None
234            The abundance axis, if available. If None, the experimental abundance axis is used.
235        freq : ndarray or None
236            The frequency axis, if available. If None, the experimental frequency axis is used.
237
238        Returns
239        -------
240        mz : ndarray
241            The extrapolated m/z axis.
242        abund : ndarray
243            The abundance axis with 0s filled.
244        freq : ndarray or None
245            The extrapolated frequency axis, if available.
246
247        Notes
248        --------
249        This function will extrapolate the mz axis by the number of datapoints specified in the settings,
250        and fill the abundance axis with 0s.
251        This should prevent peak picking issues at the spectrum edge.
252
253        """
254        # Check if the input arrays are provided
255        if mz is None or abund is None:
256            mz, abund = self.mz_exp_profile, self.abundance_profile
257            if self.has_frequency:
258                freq = self.freq_exp_profile
259            else:
260                freq = None
261        pts = self.settings.picking_point_extrapolate
262        if pts == 0:
263            return mz, abund, freq
264
265        mz = self.extrapolate_axis(mz, pts)
266        abund = pad(abund, (pts, pts), mode="constant", constant_values=(0, 0))
267        if freq is not None:
268            freq = self.extrapolate_axis(freq, pts)
269        return mz, abund, freq
270
271    def do_peak_picking(self):
272        """Perform peak picking."""
273        mz, abundance, freq = self.prepare_peak_picking_data()
274
275        if (
276            self.label == Labels.bruker_frequency
277            or self.label == Labels.midas_frequency
278        ):
279            self.calc_centroid(mz, abundance, freq)
280
281        elif self.label == Labels.thermo_profile:
282            self.calc_centroid(mz, abundance, freq)
283
284        elif self.label == Labels.bruker_profile:
285            self.calc_centroid(mz, abundance, freq)
286
287        elif self.label == Labels.booster_profile:
288            self.calc_centroid(mz, abundance, freq)
289
290        elif self.label == Labels.simulated_profile:
291            self.calc_centroid(mz, abundance, freq)
292
293        else:
294            raise Exception("Unknow mass spectrum type", self.label)
295
296    def find_minima(self, apex_index, abundance, len_abundance, right=True):
297        """Find the minima of a peak.
298
299        Parameters
300        ----------
301        apex_index : int
302            The index of the peak apex.
303        abundance : ndarray
304            The abundance values.
305        len_abundance : int
306            The length of the abundance array.
307        right : bool, optional
308            Flag indicating whether to search for minima to the right of the apex (default is True).
309
310        Returns
311        -------
312        int
313            The index of the minima.
314
315        """
316        j = apex_index
317
318        if right:
319            minima = abundance[j] > abundance[j + 1]
320        else:
321            minima = abundance[j] > abundance[j - 1]
322
323        while minima:
324            if j == 1 or j == len_abundance - 2:
325                break
326
327            if right:
328                j += 1
329
330                minima = abundance[j] >= abundance[j + 1]
331
332            else:
333                j -= 1
334                minima = abundance[j] >= abundance[j - 1]
335
336        if right:
337            return j
338        else:
339            return j
340
341    @staticmethod
342    def linear_fit_calc(intes, massa, index_term, index_sign):
343        """
344        Algebraic solution to a linear fit - roughly 25-50x faster than numpy polyfit when passing only two vals and doing a 1st order fit
345
346        Parameters
347        ----------
348        intes : ndarray
349            The intensity values.
350        massa : ndarray
351            The mass values.
352        index_term : int
353            The index of the current term.
354        index_sign : str
355            The index sign
356
357        Returns
358        -------
359        ndarray
360            The coefficients of the linear fit.
361
362        Notes
363        --------
364        This is a static method.
365        """
366        if index_sign == "+":
367            x1, x2 = massa[index_term], massa[index_term + 1]
368            y1, y2 = intes[index_term], intes[index_term + 1]
369        elif index_sign == "-":
370            x1, x2 = massa[index_term], massa[index_term - 1]
371            y1, y2 = intes[index_term], intes[index_term - 1]
372        else:
373            warnings.warn("error in linear fit calc, unknown index sign")
374
375        # Calculate the slope (m)
376        slope = (y2 - y1) / (x2 - x1)
377
378        # Calculate the intercept (b)
379        intercept = y1 - slope * x1
380
381        # The coefficients array would be [slope, intercept]
382        coefficients = array([slope, intercept])
383        return coefficients
384
385    def calculate_resolving_power(self, intes, massa, current_index):
386        """Calculate the resolving power of a peak.
387
388        Parameters
389        ----------
390        intes : ndarray
391            The intensity values.
392        massa : ndarray
393            The mass values.
394        current_index : int
395            The index of the current peak.
396
397        Returns
398        -------
399        float
400            The resolving power of the peak, or ``nan`` if a half-maximum
401            crossing does not exist on both sides (spectrum edge truncates
402            the peak).
403
404        Notes
405        --------
406        This is a conservative calculation of resolving power,
407        the peak need to be resolved at least at the half-maximum magnitude,
408        otherwise, the combined full width at half maximum is used to calculate resolving power.
409
410        """
411
412        peak_height = intes[current_index]
413        target_peak_height = peak_height / 2
414        n_points = len(intes)
415
416        index_minus = current_index
417        while index_minus > 0 and intes[index_minus] >= target_peak_height:
418            index_minus -= 1
419        if intes[index_minus] >= target_peak_height:
420            return nan
421
422        index_plus = current_index
423        while index_plus < n_points - 1 and intes[index_plus] >= target_peak_height:
424            index_plus += 1
425        if intes[index_plus] >= target_peak_height:
426            return nan
427
428        if self.mspeaks_settings.legacy_centroid_polyfit:
429            x = [massa[index_minus], massa[index_minus + 1]]
430            y = [intes[index_minus], intes[index_minus + 1]]
431            coefficients = polyfit(x, y, 1)
432        else:
433            coefficients = self.linear_fit_calc(
434                intes, massa, index_minus, index_sign="+"
435            )
436
437        a = coefficients[0]
438        b = coefficients[1]
439        if self.mspeaks_settings.legacy_resolving_power:
440            y_intercept = intes[index_minus] + (
441                (intes[index_minus + 1] - intes[index_minus]) / 2
442            )
443        else:
444            y_intercept = target_peak_height
445        massa1 = (y_intercept - b) / a
446
447        if self.mspeaks_settings.legacy_centroid_polyfit:
448            x = [massa[index_plus], massa[index_plus - 1]]
449            y = [intes[index_plus], intes[index_plus - 1]]
450            coefficients = polyfit(x, y, 1)
451        else:
452            coefficients = self.linear_fit_calc(
453                intes, massa, index_plus, index_sign="-"
454            )
455
456        a = coefficients[0]
457        b = coefficients[1]
458
459        if self.mspeaks_settings.legacy_resolving_power:
460            y_intercept = intes[index_plus - 1] + (
461                (intes[index_plus] - intes[index_plus - 1]) / 2
462            )
463        else:
464            y_intercept = target_peak_height
465
466        massa2 = (y_intercept - b) / a
467
468        if massa1 > massa2:
469            resolvingpower = massa[current_index] / (massa1 - massa2)
470
471        else:
472            resolvingpower = massa[current_index] / (massa2 - massa1)
473
474        return resolvingpower
475
476    def cal_minima(self, mass, abun):
477        """Calculate the minima of a peak.
478
479        Parameters
480        ----------
481        mass : ndarray
482            The mass values.
483        abun : ndarray
484            The abundance values.
485
486        Returns
487        -------
488        ndarray or None
489            The mass values at the minima, if found.
490
491        """
492        abun = -abun
493
494        dy = abun[1:] - abun[:-1]
495
496        # replaces nan for infinity
497        indices_nan = where(isnan(abun))[0]
498
499        if indices_nan.size:
500            abun[indices_nan] = inf
501            dy[where(isnan(dy))[0]] = inf
502
503        indexes = where((hstack((dy, 0)) < 0) & (hstack((0, dy)) > 0))[0]
504
505        if indexes.size:
506            return mass[indexes], abun[indexes]
507
508    def calc_centroid(self, mass, abund, freq):
509        """Calculate the centroid of a peak.
510
511        Parameters
512        ----------
513        mass : ndarray
514            The mass values.
515        abund : ndarray
516            The abundance values.
517        freq : ndarray or None
518            The frequency values, if available.
519
520        Returns
521        -------
522        None
523
524        """
525
526        max_height = self.mspeaks_settings.peak_height_max_percent
527        max_prominence = self.mspeaks_settings.peak_max_prominence_percent
528        min_peak_datapoints = self.mspeaks_settings.min_peak_datapoints
529        peak_derivative_threshold = self.mspeaks_settings.peak_derivative_threshold
530        max_abun = max(abund)
531        peak_height_diff = lambda hi, li: ((abund[hi] - abund[li]) / max_abun) * 100
532
533        domain = mass
534        signal = abund
535        len_signal = len(signal)
536
537        signal_threshold, factor = self.get_threshold(abund)
538        max_signal = factor
539
540        correct_baseline = False
541
542        include_indexes = sp.peak_picking_first_derivative(
543            domain,
544            signal,
545            max_height,
546            max_prominence,
547            max_signal,
548            min_peak_datapoints,
549            peak_derivative_threshold,
550            signal_threshold=signal_threshold,
551            correct_baseline=correct_baseline,
552            abun_norm=1,
553            plot_res=False,
554        )
555
556        for indexes_tuple in include_indexes:
557            apex_index = indexes_tuple[1]
558
559            peak_indexes = self.check_prominence(
560                abund, apex_index, len_signal, peak_height_diff
561            )
562
563            if peak_indexes:
564                mz_exp_centroid, freq_centr, intes_centr = self.find_apex_fit_quadratic(
565                    mass, abund, freq, apex_index
566                )
567
568                if mz_exp_centroid:
569                    peak_resolving_power = self.calculate_resolving_power(
570                        abund, mass, apex_index
571                    )
572                    s2n = intes_centr / self.baseline_noise_std
573                    self.add_mspeak(
574                        self.polarity,
575                        mz_exp_centroid,
576                        abund[apex_index],
577                        peak_resolving_power,
578                        s2n,
579                        indexes_tuple,
580                        exp_freq=freq_centr,
581                        ms_parent=self,
582                    )
583                # pyplot.plot(domain[start_index: final_index + 1], signal[start_index:final_index + 1], c='black')
584                # pyplot.show()
585
586    def get_threshold(self, intes):
587        """Get the intensity threshold for peak picking.
588
589        Parameters
590        ----------
591        intes : ndarray
592            The intensity values.
593
594        Returns
595        -------
596        float
597            The intensity threshold.
598        float
599            The factor to multiply the intensity threshold by.
600        """
601
602        intes = array(intes).astype(float)
603
604        self._resolve_noise_threshold_method()
605        noise_threshold_method = self.settings.noise_threshold_method
606
607        if noise_threshold_method == "minima":
608            if self.is_centroid:
609                warnings.warn(
610                    "Auto threshould is disabled for centroid data, returning 0"
611                )
612                factor = 1
613                abundance_threshold = 1e-20
614            # print(self.settings.noise_threshold_min_std)
615            else:
616                abundance_threshold = self.baseline_noise + (
617                    self.settings.noise_threshold_min_std * self.baseline_noise_std
618                )
619                factor = 1
620
621        elif noise_threshold_method == "signal_noise":
622            abundance_threshold = self.settings.noise_threshold_min_s2n
623            if self.is_centroid:
624                factor = 1
625            else:
626                factor = self.baseline_noise_std
627
628        elif noise_threshold_method == "relative_abundance":
629            abundance_threshold = self.settings.noise_threshold_min_relative_abundance
630            factor = intes.max() / 100
631
632        elif noise_threshold_method == "absolute_abundance":
633            abundance_threshold = self.settings.noise_threshold_absolute_abundance
634            factor = 1
635
636        elif noise_threshold_method == "log":
637            if self.is_centroid:
638                raise Exception("log noise Not tested for centroid data")
639            abundance_threshold = self.settings.noise_threshold_log_nsigma
640            factor = self.baseline_noise_std
641
642        else:
643            raise Exception(
644                "%s method was not implemented, please refer to corems.mass_spectrum.calc.NoiseCalc Class"
645                % noise_threshold_method
646            )
647
648        return abundance_threshold, factor
649
650    @staticmethod
651    def algebraic_quadratic(list_mass, list_y):
652        """
653        Find the apex of a peak - algebraically.
654        Faster than using numpy polyfit by ~28x per fit.
655
656        Parameters
657        ----------
658        list_mass : ndarray
659            list of m/z values (3 points)
660        list_y : ndarray
661            list of abundance values (3 points)
662
663        Returns
664        -------
665        a, b, c: float
666            coefficients of the quadratic equation.
667
668        Notes
669        --------
670        This is a static method.
671        """
672        x_1, x_2, x_3 = list_mass
673        y_1, y_2, y_3 = list_y
674
675        a = (
676            y_1 / ((x_1 - x_2) * (x_1 - x_3))
677            + y_2 / ((x_2 - x_1) * (x_2 - x_3))
678            + y_3 / ((x_3 - x_1) * (x_3 - x_2))
679        )
680
681        b = (
682            -y_1 * (x_2 + x_3) / ((x_1 - x_2) * (x_1 - x_3))
683            - y_2 * (x_1 + x_3) / ((x_2 - x_1) * (x_2 - x_3))
684            - y_3 * (x_1 + x_2) / ((x_3 - x_1) * (x_3 - x_2))
685        )
686
687        c = (
688            y_1 * x_2 * x_3 / ((x_1 - x_2) * (x_1 - x_3))
689            + y_2 * x_1 * x_3 / ((x_2 - x_1) * (x_2 - x_3))
690            + y_3 * x_1 * x_2 / ((x_3 - x_1) * (x_3 - x_2))
691        )
692        return a, b, c
693
694    def find_apex_fit_quadratic(self, mass, abund, freq, current_index):
695        """
696        Find the apex of a peak.
697
698        Parameters
699        ----------
700        mass : ndarray
701            The mass values.
702        abund : ndarray
703            The abundance values.
704        freq : ndarray or None
705            The frequency values, if available.
706        current_index : int
707            The index of the current peak.
708
709
710        Returns
711        -------
712        float
713            The m/z value of the peak apex.
714        float
715            The frequency value of the peak apex, if available.
716        float
717            The abundance value of the peak apex.
718
719        """
720        # calc prominence
721        # peak_indexes = self.check_prominence(abund, current_index, len_abundance, peak_height_diff )
722
723        # if not peak_indexes:
724
725        #    return None, None, None, None
726
727        # else:
728
729        # fit parabola to three most abundant datapoints
730        list_mass = [
731            mass[current_index - 1],
732            mass[current_index],
733            mass[current_index + 1],
734        ]
735        list_y = [
736            abund[current_index - 1],
737            abund[current_index],
738            abund[current_index + 1],
739        ]
740
741        if self.mspeaks_settings.legacy_centroid_polyfit:
742            z = polyfit(list_mass, list_y, 2)
743            a = z[0]
744            b = z[1]
745        else:
746            a, b, c = self.algebraic_quadratic(list_mass, list_y)
747
748        calculated = -b / (2 * a)
749
750        if calculated < 1 or int(calculated) != int(list_mass[1]):
751            mz_exp_centroid = list_mass[1]
752
753        else:
754            mz_exp_centroid = calculated
755
756        if (
757            self.label == Labels.bruker_frequency
758            or self.label == Labels.midas_frequency
759        ):
760            # fit parabola to three most abundant frequency datapoints
761            list_freq = [
762                freq[current_index - 1],
763                freq[current_index],
764                freq[current_index + 1],
765            ]
766            if self.mspeaks_settings.legacy_centroid_polyfit:
767                z = polyfit(list_mass, list_y, 2)
768                a = z[0]
769                b = z[1]
770            else:
771                a, b, c = self.algebraic_quadratic(list_mass, list_y)
772
773            calculated_freq = -b / (2 * a)
774
775            if calculated_freq < 1 or int(calculated_freq) != freq[current_index]:
776                freq_centr = list_freq[1]
777
778            else:
779                freq_centr = calculated_freq
780
781        else:
782            freq_centr = None
783
784        if self.mspeaks_settings.legacy_centroid_polyfit:
785            abundance_centroid = abund[current_index]
786        else:
787            abundance_centroid = a * mz_exp_centroid**2 + b * mz_exp_centroid + c
788
789        return mz_exp_centroid, freq_centr, abundance_centroid
790
791    def check_prominence(
792        self, abun, current_index, len_abundance, peak_height_diff
793    ) -> tuple or False:
794        """Check the prominence of a peak.
795
796        Parameters
797        ----------
798        abun : ndarray
799            The abundance values.
800        current_index : int
801            The index of the current peak.
802        len_abundance : int
803            The length of the abundance array.
804        peak_height_diff : function
805            The function to calculate the peak height difference.
806
807        Returns
808        -------
809        tuple or False
810            A tuple containing the indexes of the peak, if the prominence is above the threshold.
811            Otherwise, False.
812
813        """
814
815        final_index = self.find_minima(current_index, abun, len_abundance, right=True)
816
817        start_index = self.find_minima(current_index, abun, len_abundance, right=False)
818
819        peak_indexes = (current_index - 1, current_index, current_index + 1)
820
821        if (
822            min(
823                peak_height_diff(current_index, start_index),
824                peak_height_diff(current_index, final_index),
825            )
826            > self.mspeaks_settings.peak_min_prominence_percent
827        ):
828            return peak_indexes
829
830        else:
831            return False
832
833    def use_the_max(self, mass, abund, current_index, len_abundance, peak_height_diff):
834        """Use the max peak height as the centroid
835
836        Parameters
837        ----------
838        mass : ndarray
839            The mass values.
840        abund : ndarray
841            The abundance values.
842        current_index : int
843            The index of the current peak.
844        len_abundance : int
845            The length of the abundance array.
846        peak_height_diff : function
847            The function to calculate the peak height difference.
848
849        Returns
850        -------
851        float
852            The m/z value of the peak apex.
853        float
854            The abundance value of the peak apex.
855        tuple or None
856            A tuple containing the indexes of the peak, if the prominence is above the threshold.
857            Otherwise, None.
858        """
859
860        peak_indexes = self.check_prominence(
861            abund, current_index, len_abundance, peak_height_diff
862        )
863
864        if not peak_indexes:
865            return None, None, None
866
867        else:
868            return mass[current_index], abund[current_index], peak_indexes
869
870    def calc_centroid_legacy(self, mass, abund, freq):
871        """Legacy centroid calculation
872        Deprecated - for deletion.
873
874        """
875        warnings.warn(
876            "Legacy centroid calculation is deprecated. Please use the new centroid calculation method."
877        )
878        pass
879        if False:
880            len_abundance = len(abund)
881
882            max_abundance = max(abund)
883
884            peak_height_diff = lambda hi, li: (
885                ((abund[hi] - abund[li]) / max_abundance) * 100
886            )
887
888            abundance_threshold, factor = self.get_threshold(abund)
889            # print(abundance_threshold, factor)
890            # find indices of all peaks
891            dy = abund[1:] - abund[:-1]
892
893            # replaces nan for infi nity
894            indices_nan = where(isnan(abund))[0]
895
896            if indices_nan.size:
897                abund[indices_nan] = inf
898                dy[where(isnan(dy))[0]] = inf
899
900            indexes = where((hstack((dy, 0)) < 0) & (hstack((0, dy)) > 0))[0]
901
902            # noise threshold
903            if indexes.size and abundance_threshold is not None:
904                indexes = indexes[abund[indexes] / factor >= abundance_threshold]
905            # filter out 'peaks' within 3 points of the spectrum limits
906            # remove entries within 3 points of upper limit
907            indexes = [x for x in indexes if (len_abundance - x) > 3]
908            # remove entries within 3 points of zero
909            indexes = [x for x in indexes if x > 3]
910
911            for current_index in indexes:
912                if self.label == Labels.simulated_profile:
913                    mz_exp_centroid, intes_centr, peak_indexes = self.use_the_max(
914                        mass, abund, current_index, len_abundance, peak_height_diff
915                    )
916                    if mz_exp_centroid:
917                        peak_resolving_power = self.calculate_resolving_power(
918                            abund, mass, current_index
919                        )
920                        s2n = intes_centr / self.baseline_noise_std
921                        freq_centr = None
922                        self.add_mspeak(
923                            self.polarity,
924                            mz_exp_centroid,
925                            abund[current_index],
926                            peak_resolving_power,
927                            s2n,
928                            peak_indexes,
929                            exp_freq=freq_centr,
930                            ms_parent=self,
931                        )
932
933                else:
934                    mz_exp_centroid, freq_centr, intes_centr, peak_indexes = (
935                        self.find_apex_fit_quadratic(
936                            mass,
937                            abund,
938                            freq,
939                            current_index,
940                            len_abundance,
941                            peak_height_diff,
942                        )
943                    )
944                    if mz_exp_centroid:
945                        try:
946                            peak_resolving_power = self.calculate_resolving_power(
947                                abund, mass, current_index
948                            )
949                        except IndexError:
950                            print("index error, skipping peak")
951                            continue
952
953                        s2n = intes_centr / self.baseline_noise_std
954                        self.add_mspeak(
955                            self.polarity,
956                            mz_exp_centroid,
957                            abund[current_index],
958                            peak_resolving_power,
959                            s2n,
960                            peak_indexes,
961                            exp_freq=freq_centr,
962                            ms_parent=self,
963                        )

Class for peak picking.

Parameters
  • None
Attributes
  • None
Methods
  • prepare_peak_picking_data(). Prepare the mz, abundance, and frequence data for peak picking.
  • cut_mz_domain_peak_picking(). Cut the m/z domain for peak picking.
  • extrapolate_axes_for_pp(mz=None, abund=None, freq=None). Extrapolate the m/z axis and fill the abundance axis with 0s.
  • do_peak_picking(). Perform peak picking.
  • find_minima(apex_index, abundance, len_abundance, right=True). Find the minima of a peak.
  • linear_fit_calc(intes, massa, index_term, index_sign). Algebraic solution to a linear fit.
  • calculate_resolving_power(intes, massa, current_index). Calculate the resolving power of a peak.
  • cal_minima(mass, abun). Calculate the minima of a peak.
  • calc_centroid(mass, abund, freq). Calculate the centroid of a peak.
  • get_threshold(intes). Get the intensity threshold for peak picking.
  • algebraic_quadratic(list_mass, list_y). Find the apex of a peak - algebraically.
  • find_apex_fit_quadratic(mass, abund, freq, current_index). Find the apex of a peak.
  • check_prominence(abun, current_index, len_abundance, peak_height_diff). Check the prominence of a peak.
  • use_the_max(mass, abund, current_index, len_abundance, peak_height_diff). Use the max peak height as the centroid.
  • calc_centroid_legacy(mass, abund, freq). Legacy centroid calculation. Deprecated - for deletion.
def prepare_peak_picking_data(self): 
70    def prepare_peak_picking_data(self):
71        """Prepare the data for peak picking.
72
73        This function will prepare the m/z, abundance, and frequency data for peak picking according to the settings.
74
75        Returns
76        -------
77        mz : ndarray
78            The m/z axis.
79        abundance : ndarray
80            The abundance axis.
81        freq : ndarray or None
82            The frequency axis, if available.
83        """
84        # First apply cut_mz_domain_peak_picking
85        mz, abundance, freq = self.cut_mz_domain_peak_picking()
86
87        # Then extrapolate the axes for peak picking
88        if self.settings.picking_point_extrapolate > 0:
89            mz, abundance, freq = self.extrapolate_axes_for_pp(mz, abundance, freq)
90        return mz, abundance, freq

Prepare the data for peak picking.

This function will prepare the m/z, abundance, and frequency data for peak picking according to the settings.

Returns
  • mz (ndarray): The m/z axis.
  • abundance (ndarray): The abundance axis.
  • freq (ndarray or None): The frequency axis, if available.
def cut_mz_domain_peak_picking(self): 
 92    def cut_mz_domain_peak_picking(self):
 93        """
 94        Cut the m/z domain for peak picking.
 95
 96        Simplified function
 97
 98        Returns
 99        -------
100        mz_domain_X_low_cutoff : ndarray
101            The m/z values within the specified range.
102        mz_domain_low_Y_cutoff : ndarray
103            The abundance values within the specified range.
104        freq_domain_low_Y_cutoff : ndarray or None
105            The frequency values within the specified range, if available.
106
107        """
108        max_picking_mz = self.settings.max_picking_mz
109        min_picking_mz = self.settings.min_picking_mz
110
111        # min_start =  where(self.mz_exp_profile  > min_picking_mz)[0][0]
112        # max_final =  where(self.mz_exp_profile < max_picking_mz)[-1][-1]
113        min_start = searchsorted(a=self.mz_exp_profile, v=min_picking_mz)
114        max_final = searchsorted(a=self.mz_exp_profile, v=max_picking_mz)
115
116        if self.has_frequency:
117            if self.freq_exp_profile.any():
118                return (
119                    self.mz_exp_profile[min_start:max_final],
120                    self.abundance_profile[min_start:max_final],
121                    self.freq_exp_profile[min_start:max_final],
122                )
123
124        else:
125            return (
126                self.mz_exp_profile[min_start:max_final],
127                self.abundance_profile[min_start:max_final],
128                None,
129            )

Cut the m/z domain for peak picking.

Simplified function

Returns
  • mz_domain_X_low_cutoff (ndarray): The m/z values within the specified range.
  • mz_domain_low_Y_cutoff (ndarray): The abundance values within the specified range.
  • freq_domain_low_Y_cutoff (ndarray or None): The frequency values within the specified range, if available.
def legacy_cut_mz_domain_peak_picking(self): 
131    def legacy_cut_mz_domain_peak_picking(self):
132        """
133        Cut the m/z domain for peak picking.
134        DEPRECATED
135        Returns
136        -------
137        mz_domain_X_low_cutoff : ndarray
138            The m/z values within the specified range.
139        mz_domain_low_Y_cutoff : ndarray
140            The abundance values within the specified range.
141        freq_domain_low_Y_cutoff : ndarray or None
142            The frequency values within the specified range, if available.
143
144        """
145        max_picking_mz = self.settings.max_picking_mz
146        min_picking_mz = self.settings.min_picking_mz
147
148        min_final = where(self.mz_exp_profile > min_picking_mz)[-1][-1]
149        min_start = where(self.mz_exp_profile > min_picking_mz)[0][0]
150
151        (
152            mz_domain_X_low_cutoff,
153            mz_domain_low_Y_cutoff,
154        ) = (
155            self.mz_exp_profile[min_start:min_final],
156            self.abundance_profile[min_start:min_final],
157        )
158
159        max_final = where(self.mz_exp_profile < max_picking_mz)[-1][-1]
160        max_start = where(self.mz_exp_profile < max_picking_mz)[0][0]
161
162        if self.has_frequency:
163            if self.freq_exp_profile.any():
164                freq_domain_low_Y_cutoff = self.freq_exp_profile[min_start:min_final]
165
166                return (
167                    mz_domain_X_low_cutoff[max_start:max_final],
168                    mz_domain_low_Y_cutoff[max_start:max_final],
169                    freq_domain_low_Y_cutoff[max_start:max_final],
170                )
171
172        else:
173            return (
174                mz_domain_X_low_cutoff[max_start:max_final],
175                mz_domain_low_Y_cutoff[max_start:max_final],
176                None,
177            )

Cut the m/z domain for peak picking. DEPRECATED

Returns
  • mz_domain_X_low_cutoff (ndarray): The m/z values within the specified range.
  • mz_domain_low_Y_cutoff (ndarray): The abundance values within the specified range.
  • freq_domain_low_Y_cutoff (ndarray or None): The frequency values within the specified range, if available.
@staticmethod
def extrapolate_axis(initial_array, pts): 
179    @staticmethod
180    def extrapolate_axis(initial_array, pts):
181        """
182        This function will extrapolate an input array in both directions by N pts.
183
184        Parameters
185        ----------
186        initial_array : ndarray
187            The input array.
188        pts : int
189            The number of points to extrapolate.
190
191        Returns
192        -------
193        ndarray
194            The extrapolated array.
195
196        Notes
197        --------
198        This is a static method.
199        """
200        initial_array_len = len(initial_array)
201        right_delta = initial_array[-1] - initial_array[-2]
202        left_delta = initial_array[1] - initial_array[0]
203
204        # Create an array with extra space for extrapolation
205        pad_array = zeros(initial_array_len + 2 * pts)
206
207        # Copy original array into the middle of the padded array
208        pad_array[pts : pts + initial_array_len] = initial_array
209
210        # Extrapolate the right side
211        for pt in range(pts):
212            final_value = initial_array[-1]
213            value_to_add = right_delta * (pt + 1)
214            new_value = final_value + value_to_add
215            pad_array[initial_array_len + pts + pt] = new_value
216
217        # Extrapolate the left side
218        for pt in range(pts):
219            first_value = initial_array[0]
220            value_to_subtract = left_delta * (pt + 1)
221            new_value = first_value - value_to_subtract
222            pad_array[pts - pt - 1] = new_value
223
224        return pad_array

This function will extrapolate an input array in both directions by N pts.

Parameters
  • initial_array (ndarray): The input array.
  • pts (int): The number of points to extrapolate.
Returns
  • ndarray: The extrapolated array.
Notes

This is a static method.

def extrapolate_axes_for_pp(self, mz=None, abund=None, freq=None): 
226    def extrapolate_axes_for_pp(self, mz=None, abund=None, freq=None):
227        """Extrapolate the m/z axis and fill the abundance axis with 0s.
228
229        Parameters
230        ----------
231        mz : ndarray or None
232            The m/z axis, if available. If None, the experimental m/z axis is used.
233        abund : ndarray or None
234            The abundance axis, if available. If None, the experimental abundance axis is used.
235        freq : ndarray or None
236            The frequency axis, if available. If None, the experimental frequency axis is used.
237
238        Returns
239        -------
240        mz : ndarray
241            The extrapolated m/z axis.
242        abund : ndarray
243            The abundance axis with 0s filled.
244        freq : ndarray or None
245            The extrapolated frequency axis, if available.
246
247        Notes
248        --------
249        This function will extrapolate the mz axis by the number of datapoints specified in the settings,
250        and fill the abundance axis with 0s.
251        This should prevent peak picking issues at the spectrum edge.
252
253        """
254        # Check if the input arrays are provided
255        if mz is None or abund is None:
256            mz, abund = self.mz_exp_profile, self.abundance_profile
257            if self.has_frequency:
258                freq = self.freq_exp_profile
259            else:
260                freq = None
261        pts = self.settings.picking_point_extrapolate
262        if pts == 0:
263            return mz, abund, freq
264
265        mz = self.extrapolate_axis(mz, pts)
266        abund = pad(abund, (pts, pts), mode="constant", constant_values=(0, 0))
267        if freq is not None:
268            freq = self.extrapolate_axis(freq, pts)
269        return mz, abund, freq

Extrapolate the m/z axis and fill the abundance axis with 0s.

Parameters
  • mz (ndarray or None): The m/z axis, if available. If None, the experimental m/z axis is used.
  • abund (ndarray or None): The abundance axis, if available. If None, the experimental abundance axis is used.
  • freq (ndarray or None): The frequency axis, if available. If None, the experimental frequency axis is used.
Returns
  • mz (ndarray): The extrapolated m/z axis.
  • abund (ndarray): The abundance axis with 0s filled.
  • freq (ndarray or None): The extrapolated frequency axis, if available.
Notes

This function will extrapolate the mz axis by the number of datapoints specified in the settings, and fill the abundance axis with 0s. This should prevent peak picking issues at the spectrum edge.

def do_peak_picking(self): 
271    def do_peak_picking(self):
272        """Perform peak picking."""
273        mz, abundance, freq = self.prepare_peak_picking_data()
274
275        if (
276            self.label == Labels.bruker_frequency
277            or self.label == Labels.midas_frequency
278        ):
279            self.calc_centroid(mz, abundance, freq)
280
281        elif self.label == Labels.thermo_profile:
282            self.calc_centroid(mz, abundance, freq)
283
284        elif self.label == Labels.bruker_profile:
285            self.calc_centroid(mz, abundance, freq)
286
287        elif self.label == Labels.booster_profile:
288            self.calc_centroid(mz, abundance, freq)
289
290        elif self.label == Labels.simulated_profile:
291            self.calc_centroid(mz, abundance, freq)
292
293        else:
294            raise Exception("Unknow mass spectrum type", self.label)

Perform peak picking.

def find_minima(self, apex_index, abundance, len_abundance, right=True): 
296    def find_minima(self, apex_index, abundance, len_abundance, right=True):
297        """Find the minima of a peak.
298
299        Parameters
300        ----------
301        apex_index : int
302            The index of the peak apex.
303        abundance : ndarray
304            The abundance values.
305        len_abundance : int
306            The length of the abundance array.
307        right : bool, optional
308            Flag indicating whether to search for minima to the right of the apex (default is True).
309
310        Returns
311        -------
312        int
313            The index of the minima.
314
315        """
316        j = apex_index
317
318        if right:
319            minima = abundance[j] > abundance[j + 1]
320        else:
321            minima = abundance[j] > abundance[j - 1]
322
323        while minima:
324            if j == 1 or j == len_abundance - 2:
325                break
326
327            if right:
328                j += 1
329
330                minima = abundance[j] >= abundance[j + 1]
331
332            else:
333                j -= 1
334                minima = abundance[j] >= abundance[j - 1]
335
336        if right:
337            return j
338        else:
339            return j

Find the minima of a peak.

Parameters
  • apex_index (int): The index of the peak apex.
  • abundance (ndarray): The abundance values.
  • len_abundance (int): The length of the abundance array.
  • right (bool, optional): Flag indicating whether to search for minima to the right of the apex (default is True).
Returns
  • int: The index of the minima.
@staticmethod
def linear_fit_calc(intes, massa, index_term, index_sign): 
341    @staticmethod
342    def linear_fit_calc(intes, massa, index_term, index_sign):
343        """
344        Algebraic solution to a linear fit - roughly 25-50x faster than numpy polyfit when passing only two vals and doing a 1st order fit
345
346        Parameters
347        ----------
348        intes : ndarray
349            The intensity values.
350        massa : ndarray
351            The mass values.
352        index_term : int
353            The index of the current term.
354        index_sign : str
355            The index sign
356
357        Returns
358        -------
359        ndarray
360            The coefficients of the linear fit.
361
362        Notes
363        --------
364        This is a static method.
365        """
366        if index_sign == "+":
367            x1, x2 = massa[index_term], massa[index_term + 1]
368            y1, y2 = intes[index_term], intes[index_term + 1]
369        elif index_sign == "-":
370            x1, x2 = massa[index_term], massa[index_term - 1]
371            y1, y2 = intes[index_term], intes[index_term - 1]
372        else:
373            warnings.warn("error in linear fit calc, unknown index sign")
374
375        # Calculate the slope (m)
376        slope = (y2 - y1) / (x2 - x1)
377
378        # Calculate the intercept (b)
379        intercept = y1 - slope * x1
380
381        # The coefficients array would be [slope, intercept]
382        coefficients = array([slope, intercept])
383        return coefficients

Algebraic solution to a linear fit - roughly 25-50x faster than numpy polyfit when passing only two vals and doing a 1st order fit

Parameters
  • intes (ndarray): The intensity values.
  • massa (ndarray): The mass values.
  • index_term (int): The index of the current term.
  • index_sign (str): The index sign
Returns
  • ndarray: The coefficients of the linear fit.
Notes

This is a static method.

def calculate_resolving_power(self, intes, massa, current_index): 
385    def calculate_resolving_power(self, intes, massa, current_index):
386        """Calculate the resolving power of a peak.
387
388        Parameters
389        ----------
390        intes : ndarray
391            The intensity values.
392        massa : ndarray
393            The mass values.
394        current_index : int
395            The index of the current peak.
396
397        Returns
398        -------
399        float
400            The resolving power of the peak, or ``nan`` if a half-maximum
401            crossing does not exist on both sides (spectrum edge truncates
402            the peak).
403
404        Notes
405        --------
406        This is a conservative calculation of resolving power,
407        the peak need to be resolved at least at the half-maximum magnitude,
408        otherwise, the combined full width at half maximum is used to calculate resolving power.
409
410        """
411
412        peak_height = intes[current_index]
413        target_peak_height = peak_height / 2
414        n_points = len(intes)
415
416        index_minus = current_index
417        while index_minus > 0 and intes[index_minus] >= target_peak_height:
418            index_minus -= 1
419        if intes[index_minus] >= target_peak_height:
420            return nan
421
422        index_plus = current_index
423        while index_plus < n_points - 1 and intes[index_plus] >= target_peak_height:
424            index_plus += 1
425        if intes[index_plus] >= target_peak_height:
426            return nan
427
428        if self.mspeaks_settings.legacy_centroid_polyfit:
429            x = [massa[index_minus], massa[index_minus + 1]]
430            y = [intes[index_minus], intes[index_minus + 1]]
431            coefficients = polyfit(x, y, 1)
432        else:
433            coefficients = self.linear_fit_calc(
434                intes, massa, index_minus, index_sign="+"
435            )
436
437        a = coefficients[0]
438        b = coefficients[1]
439        if self.mspeaks_settings.legacy_resolving_power:
440            y_intercept = intes[index_minus] + (
441                (intes[index_minus + 1] - intes[index_minus]) / 2
442            )
443        else:
444            y_intercept = target_peak_height
445        massa1 = (y_intercept - b) / a
446
447        if self.mspeaks_settings.legacy_centroid_polyfit:
448            x = [massa[index_plus], massa[index_plus - 1]]
449            y = [intes[index_plus], intes[index_plus - 1]]
450            coefficients = polyfit(x, y, 1)
451        else:
452            coefficients = self.linear_fit_calc(
453                intes, massa, index_plus, index_sign="-"
454            )
455
456        a = coefficients[0]
457        b = coefficients[1]
458
459        if self.mspeaks_settings.legacy_resolving_power:
460            y_intercept = intes[index_plus - 1] + (
461                (intes[index_plus] - intes[index_plus - 1]) / 2
462            )
463        else:
464            y_intercept = target_peak_height
465
466        massa2 = (y_intercept - b) / a
467
468        if massa1 > massa2:
469            resolvingpower = massa[current_index] / (massa1 - massa2)
470
471        else:
472            resolvingpower = massa[current_index] / (massa2 - massa1)
473
474        return resolvingpower

Calculate the resolving power of a peak.

Parameters
  • intes (ndarray): The intensity values.
  • massa (ndarray): The mass values.
  • current_index (int): The index of the current peak.
Returns
  • float: The resolving power of the peak, or nan if a half-maximum crossing does not exist on both sides (spectrum edge truncates the peak).
Notes

This is a conservative calculation of resolving power, the peak need to be resolved at least at the half-maximum magnitude, otherwise, the combined full width at half maximum is used to calculate resolving power.

def cal_minima(self, mass, abun): 
476    def cal_minima(self, mass, abun):
477        """Calculate the minima of a peak.
478
479        Parameters
480        ----------
481        mass : ndarray
482            The mass values.
483        abun : ndarray
484            The abundance values.
485
486        Returns
487        -------
488        ndarray or None
489            The mass values at the minima, if found.
490
491        """
492        abun = -abun
493
494        dy = abun[1:] - abun[:-1]
495
496        # replaces nan for infinity
497        indices_nan = where(isnan(abun))[0]
498
499        if indices_nan.size:
500            abun[indices_nan] = inf
501            dy[where(isnan(dy))[0]] = inf
502
503        indexes = where((hstack((dy, 0)) < 0) & (hstack((0, dy)) > 0))[0]
504
505        if indexes.size:
506            return mass[indexes], abun[indexes]

Calculate the minima of a peak.

Parameters
  • mass (ndarray): The mass values.
  • abun (ndarray): The abundance values.
Returns
  • ndarray or None: The mass values at the minima, if found.
def calc_centroid(self, mass, abund, freq): 
508    def calc_centroid(self, mass, abund, freq):
509        """Calculate the centroid of a peak.
510
511        Parameters
512        ----------
513        mass : ndarray
514            The mass values.
515        abund : ndarray
516            The abundance values.
517        freq : ndarray or None
518            The frequency values, if available.
519
520        Returns
521        -------
522        None
523
524        """
525
526        max_height = self.mspeaks_settings.peak_height_max_percent
527        max_prominence = self.mspeaks_settings.peak_max_prominence_percent
528        min_peak_datapoints = self.mspeaks_settings.min_peak_datapoints
529        peak_derivative_threshold = self.mspeaks_settings.peak_derivative_threshold
530        max_abun = max(abund)
531        peak_height_diff = lambda hi, li: ((abund[hi] - abund[li]) / max_abun) * 100
532
533        domain = mass
534        signal = abund
535        len_signal = len(signal)
536
537        signal_threshold, factor = self.get_threshold(abund)
538        max_signal = factor
539
540        correct_baseline = False
541
542        include_indexes = sp.peak_picking_first_derivative(
543            domain,
544            signal,
545            max_height,
546            max_prominence,
547            max_signal,
548            min_peak_datapoints,
549            peak_derivative_threshold,
550            signal_threshold=signal_threshold,
551            correct_baseline=correct_baseline,
552            abun_norm=1,
553            plot_res=False,
554        )
555
556        for indexes_tuple in include_indexes:
557            apex_index = indexes_tuple[1]
558
559            peak_indexes = self.check_prominence(
560                abund, apex_index, len_signal, peak_height_diff
561            )
562
563            if peak_indexes:
564                mz_exp_centroid, freq_centr, intes_centr = self.find_apex_fit_quadratic(
565                    mass, abund, freq, apex_index
566                )
567
568                if mz_exp_centroid:
569                    peak_resolving_power = self.calculate_resolving_power(
570                        abund, mass, apex_index
571                    )
572                    s2n = intes_centr / self.baseline_noise_std
573                    self.add_mspeak(
574                        self.polarity,
575                        mz_exp_centroid,
576                        abund[apex_index],
577                        peak_resolving_power,
578                        s2n,
579                        indexes_tuple,
580                        exp_freq=freq_centr,
581                        ms_parent=self,
582                    )
583                # pyplot.plot(domain[start_index: final_index + 1], signal[start_index:final_index + 1], c='black')
584                # pyplot.show()

Calculate the centroid of a peak.

Parameters
  • mass (ndarray): The mass values.
  • abund (ndarray): The abundance values.
  • freq (ndarray or None): The frequency values, if available.
Returns
  • None
def get_threshold(self, intes): 
586    def get_threshold(self, intes):
587        """Get the intensity threshold for peak picking.
588
589        Parameters
590        ----------
591        intes : ndarray
592            The intensity values.
593
594        Returns
595        -------
596        float
597            The intensity threshold.
598        float
599            The factor to multiply the intensity threshold by.
600        """
601
602        intes = array(intes).astype(float)
603
604        self._resolve_noise_threshold_method()
605        noise_threshold_method = self.settings.noise_threshold_method
606
607        if noise_threshold_method == "minima":
608            if self.is_centroid:
609                warnings.warn(
610                    "Auto threshould is disabled for centroid data, returning 0"
611                )
612                factor = 1
613                abundance_threshold = 1e-20
614            # print(self.settings.noise_threshold_min_std)
615            else:
616                abundance_threshold = self.baseline_noise + (
617                    self.settings.noise_threshold_min_std * self.baseline_noise_std
618                )
619                factor = 1
620
621        elif noise_threshold_method == "signal_noise":
622            abundance_threshold = self.settings.noise_threshold_min_s2n
623            if self.is_centroid:
624                factor = 1
625            else:
626                factor = self.baseline_noise_std
627
628        elif noise_threshold_method == "relative_abundance":
629            abundance_threshold = self.settings.noise_threshold_min_relative_abundance
630            factor = intes.max() / 100
631
632        elif noise_threshold_method == "absolute_abundance":
633            abundance_threshold = self.settings.noise_threshold_absolute_abundance
634            factor = 1
635
636        elif noise_threshold_method == "log":
637            if self.is_centroid:
638                raise Exception("log noise Not tested for centroid data")
639            abundance_threshold = self.settings.noise_threshold_log_nsigma
640            factor = self.baseline_noise_std
641
642        else:
643            raise Exception(
644                "%s method was not implemented, please refer to corems.mass_spectrum.calc.NoiseCalc Class"
645                % noise_threshold_method
646            )
647
648        return abundance_threshold, factor

Get the intensity threshold for peak picking.

Parameters
  • intes (ndarray): The intensity values.
Returns
  • float: The intensity threshold.
  • float: The factor to multiply the intensity threshold by.
@staticmethod
def algebraic_quadratic(list_mass, list_y): 
650    @staticmethod
651    def algebraic_quadratic(list_mass, list_y):
652        """
653        Find the apex of a peak - algebraically.
654        Faster than using numpy polyfit by ~28x per fit.
655
656        Parameters
657        ----------
658        list_mass : ndarray
659            list of m/z values (3 points)
660        list_y : ndarray
661            list of abundance values (3 points)
662
663        Returns
664        -------
665        a, b, c: float
666            coefficients of the quadratic equation.
667
668        Notes
669        --------
670        This is a static method.
671        """
672        x_1, x_2, x_3 = list_mass
673        y_1, y_2, y_3 = list_y
674
675        a = (
676            y_1 / ((x_1 - x_2) * (x_1 - x_3))
677            + y_2 / ((x_2 - x_1) * (x_2 - x_3))
678            + y_3 / ((x_3 - x_1) * (x_3 - x_2))
679        )
680
681        b = (
682            -y_1 * (x_2 + x_3) / ((x_1 - x_2) * (x_1 - x_3))
683            - y_2 * (x_1 + x_3) / ((x_2 - x_1) * (x_2 - x_3))
684            - y_3 * (x_1 + x_2) / ((x_3 - x_1) * (x_3 - x_2))
685        )
686
687        c = (
688            y_1 * x_2 * x_3 / ((x_1 - x_2) * (x_1 - x_3))
689            + y_2 * x_1 * x_3 / ((x_2 - x_1) * (x_2 - x_3))
690            + y_3 * x_1 * x_2 / ((x_3 - x_1) * (x_3 - x_2))
691        )
692        return a, b, c

Find the apex of a peak - algebraically. Faster than using numpy polyfit by ~28x per fit.

Parameters
  • list_mass (ndarray): list of m/z values (3 points)
  • list_y (ndarray): list of abundance values (3 points)
Returns
  • a, b, c (float): coefficients of the quadratic equation.
Notes

This is a static method.

def find_apex_fit_quadratic(self, mass, abund, freq, current_index): 
694    def find_apex_fit_quadratic(self, mass, abund, freq, current_index):
695        """
696        Find the apex of a peak.
697
698        Parameters
699        ----------
700        mass : ndarray
701            The mass values.
702        abund : ndarray
703            The abundance values.
704        freq : ndarray or None
705            The frequency values, if available.
706        current_index : int
707            The index of the current peak.
708
709
710        Returns
711        -------
712        float
713            The m/z value of the peak apex.
714        float
715            The frequency value of the peak apex, if available.
716        float
717            The abundance value of the peak apex.
718
719        """
720        # calc prominence
721        # peak_indexes = self.check_prominence(abund, current_index, len_abundance, peak_height_diff )
722
723        # if not peak_indexes:
724
725        #    return None, None, None, None
726
727        # else:
728
729        # fit parabola to three most abundant datapoints
730        list_mass = [
731            mass[current_index - 1],
732            mass[current_index],
733            mass[current_index + 1],
734        ]
735        list_y = [
736            abund[current_index - 1],
737            abund[current_index],
738            abund[current_index + 1],
739        ]
740
741        if self.mspeaks_settings.legacy_centroid_polyfit:
742            z = polyfit(list_mass, list_y, 2)
743            a = z[0]
744            b = z[1]
745        else:
746            a, b, c = self.algebraic_quadratic(list_mass, list_y)
747
748        calculated = -b / (2 * a)
749
750        if calculated < 1 or int(calculated) != int(list_mass[1]):
751            mz_exp_centroid = list_mass[1]
752
753        else:
754            mz_exp_centroid = calculated
755
756        if (
757            self.label == Labels.bruker_frequency
758            or self.label == Labels.midas_frequency
759        ):
760            # fit parabola to three most abundant frequency datapoints
761            list_freq = [
762                freq[current_index - 1],
763                freq[current_index],
764                freq[current_index + 1],
765            ]
766            if self.mspeaks_settings.legacy_centroid_polyfit:
767                z = polyfit(list_mass, list_y, 2)
768                a = z[0]
769                b = z[1]
770            else:
771                a, b, c = self.algebraic_quadratic(list_mass, list_y)
772
773            calculated_freq = -b / (2 * a)
774
775            if calculated_freq < 1 or int(calculated_freq) != freq[current_index]:
776                freq_centr = list_freq[1]
777
778            else:
779                freq_centr = calculated_freq
780
781        else:
782            freq_centr = None
783
784        if self.mspeaks_settings.legacy_centroid_polyfit:
785            abundance_centroid = abund[current_index]
786        else:
787            abundance_centroid = a * mz_exp_centroid**2 + b * mz_exp_centroid + c
788
789        return mz_exp_centroid, freq_centr, abundance_centroid

Find the apex of a peak.

Parameters
  • mass (ndarray): The mass values.
  • abund (ndarray): The abundance values.
  • freq (ndarray or None): The frequency values, if available.
  • current_index (int): The index of the current peak.
Returns
  • float: The m/z value of the peak apex.
  • float: The frequency value of the peak apex, if available.
  • float: The abundance value of the peak apex.
def check_prominence(self, abun, current_index, len_abundance, peak_height_diff) -> tuple: 
791    def check_prominence(
792        self, abun, current_index, len_abundance, peak_height_diff
793    ) -> tuple or False:
794        """Check the prominence of a peak.
795
796        Parameters
797        ----------
798        abun : ndarray
799            The abundance values.
800        current_index : int
801            The index of the current peak.
802        len_abundance : int
803            The length of the abundance array.
804        peak_height_diff : function
805            The function to calculate the peak height difference.
806
807        Returns
808        -------
809        tuple or False
810            A tuple containing the indexes of the peak, if the prominence is above the threshold.
811            Otherwise, False.
812
813        """
814
815        final_index = self.find_minima(current_index, abun, len_abundance, right=True)
816
817        start_index = self.find_minima(current_index, abun, len_abundance, right=False)
818
819        peak_indexes = (current_index - 1, current_index, current_index + 1)
820
821        if (
822            min(
823                peak_height_diff(current_index, start_index),
824                peak_height_diff(current_index, final_index),
825            )
826            > self.mspeaks_settings.peak_min_prominence_percent
827        ):
828            return peak_indexes
829
830        else:
831            return False

Check the prominence of a peak.

Parameters
  • abun (ndarray): The abundance values.
  • current_index (int): The index of the current peak.
  • len_abundance (int): The length of the abundance array.
  • peak_height_diff (function): The function to calculate the peak height difference.
Returns
  • tuple or False: A tuple containing the indexes of the peak, if the prominence is above the threshold. Otherwise, False.
def use_the_max(self, mass, abund, current_index, len_abundance, peak_height_diff): 
833    def use_the_max(self, mass, abund, current_index, len_abundance, peak_height_diff):
834        """Use the max peak height as the centroid
835
836        Parameters
837        ----------
838        mass : ndarray
839            The mass values.
840        abund : ndarray
841            The abundance values.
842        current_index : int
843            The index of the current peak.
844        len_abundance : int
845            The length of the abundance array.
846        peak_height_diff : function
847            The function to calculate the peak height difference.
848
849        Returns
850        -------
851        float
852            The m/z value of the peak apex.
853        float
854            The abundance value of the peak apex.
855        tuple or None
856            A tuple containing the indexes of the peak, if the prominence is above the threshold.
857            Otherwise, None.
858        """
859
860        peak_indexes = self.check_prominence(
861            abund, current_index, len_abundance, peak_height_diff
862        )
863
864        if not peak_indexes:
865            return None, None, None
866
867        else:
868            return mass[current_index], abund[current_index], peak_indexes

Use the max peak height as the centroid

Parameters
  • mass (ndarray): The mass values.
  • abund (ndarray): The abundance values.
  • current_index (int): The index of the current peak.
  • len_abundance (int): The length of the abundance array.
  • peak_height_diff (function): The function to calculate the peak height difference.
Returns
  • float: The m/z value of the peak apex.
  • float: The abundance value of the peak apex.
  • tuple or None: A tuple containing the indexes of the peak, if the prominence is above the threshold. Otherwise, None.
def calc_centroid_legacy(self, mass, abund, freq): 
870    def calc_centroid_legacy(self, mass, abund, freq):
871        """Legacy centroid calculation
872        Deprecated - for deletion.
873
874        """
875        warnings.warn(
876            "Legacy centroid calculation is deprecated. Please use the new centroid calculation method."
877        )
878        pass
879        if False:
880            len_abundance = len(abund)
881
882            max_abundance = max(abund)
883
884            peak_height_diff = lambda hi, li: (
885                ((abund[hi] - abund[li]) / max_abundance) * 100
886            )
887
888            abundance_threshold, factor = self.get_threshold(abund)
889            # print(abundance_threshold, factor)
890            # find indices of all peaks
891            dy = abund[1:] - abund[:-1]
892
893            # replaces nan for infi nity
894            indices_nan = where(isnan(abund))[0]
895
896            if indices_nan.size:
897                abund[indices_nan] = inf
898                dy[where(isnan(dy))[0]] = inf
899
900            indexes = where((hstack((dy, 0)) < 0) & (hstack((0, dy)) > 0))[0]
901
902            # noise threshold
903            if indexes.size and abundance_threshold is not None:
904                indexes = indexes[abund[indexes] / factor >= abundance_threshold]
905            # filter out 'peaks' within 3 points of the spectrum limits
906            # remove entries within 3 points of upper limit
907            indexes = [x for x in indexes if (len_abundance - x) > 3]
908            # remove entries within 3 points of zero
909            indexes = [x for x in indexes if x > 3]
910
911            for current_index in indexes:
912                if self.label == Labels.simulated_profile:
913                    mz_exp_centroid, intes_centr, peak_indexes = self.use_the_max(
914                        mass, abund, current_index, len_abundance, peak_height_diff
915                    )
916                    if mz_exp_centroid:
917                        peak_resolving_power = self.calculate_resolving_power(
918                            abund, mass, current_index
919                        )
920                        s2n = intes_centr / self.baseline_noise_std
921                        freq_centr = None
922                        self.add_mspeak(
923                            self.polarity,
924                            mz_exp_centroid,
925                            abund[current_index],
926                            peak_resolving_power,
927                            s2n,
928                            peak_indexes,
929                            exp_freq=freq_centr,
930                            ms_parent=self,
931                        )
932
933                else:
934                    mz_exp_centroid, freq_centr, intes_centr, peak_indexes = (
935                        self.find_apex_fit_quadratic(
936                            mass,
937                            abund,
938                            freq,
939                            current_index,
940                            len_abundance,
941                            peak_height_diff,
942                        )
943                    )
944                    if mz_exp_centroid:
945                        try:
946                            peak_resolving_power = self.calculate_resolving_power(
947                                abund, mass, current_index
948                            )
949                        except IndexError:
950                            print("index error, skipping peak")
951                            continue
952
953                        s2n = intes_centr / self.baseline_noise_std
954                        self.add_mspeak(
955                            self.polarity,
956                            mz_exp_centroid,
957                            abund[current_index],
958                            peak_resolving_power,
959                            s2n,
960                            peak_indexes,
961                            exp_freq=freq_centr,
962                            ms_parent=self,
963                        )

Legacy centroid calculation Deprecated - for deletion.