corems.mass_spectrum.calc.MassSpectrumCalc

  1__author__ = "Yuri E. Corilo"
  2__date__ = "Jun 27, 2019"
  3
  4from numpy import power, multiply, sqrt, array, mean
  5from corems.mass_spectrum.calc.NoiseCalc import NoiseThresholdCalc
  6from corems.mass_spectrum.calc.PeakPicking import PeakPicking
  7
  8
  9class MassSpecCalc(PeakPicking, NoiseThresholdCalc):
 10    """Class for Mass Spectrum Calculations
 11
 12    Class including numerical calculations related to mass spectrum class
 13    Inherited PeakPicking and NoiseThresholdCalc ensuring its methods are
 14    available to the instantiated mass spectrum class object
 15
 16    Parameters
 17    -------
 18    mass_spectrum : MassSpectrum
 19        CoreMS mass spectrum object
 20
 21    Attributes
 22    --------
 23    All Attributes are derivative from the MassSpecBase Class
 24
 25    Methods
 26    --------
 27    * check_mspeaks().
 28        Check if the mspeaks attribute is populated
 29    * sort_by_abundance().
 30        Sort the mspeaks by abundance
 31    * percentile_assigned(report_error=False).
 32        Calculate the percentage of assigned peaks
 33    * resolving_power_calc(B, T).
 34        Calculate the resolving power
 35    * number_average_molecular_weight(profile=False).
 36        Calculate the number average molecular weight
 37    * weight_average_molecular_weight(profile=False).
 38        Calculate the weight average molecular weight
 39    """
 40
 41    def percentile_assigned(self, report_error: bool = False, mute_output: bool = False):
 42        """Percentage of peaks which are assigned
 43
 44        Parameters
 45        -----------
 46        report_error: bool, optional
 47            Report the error of the assigned peaks. Default is False.
 48        mute_output: bool, optional
 49            Override the verbose setting. Default is False.
 50            If True, the function will silence results
 51        """
 52        verbose = self.parameters.mass_spectrum.verbose_processing
 53        assign_abun = 0
 54        not_assign_abun = 0
 55        i = 0
 56        j = 0
 57        if report_error:
 58            error = []
 59        for mspeak in self.sort_by_abundance():
 60            if mspeak.is_assigned:
 61                i += 1
 62                assign_abun += mspeak.abundance
 63                if report_error:
 64                    error.append(mspeak.best_molecular_formula_candidate.mz_error)
 65
 66            else:
 67                j += 1
 68                not_assign_abun += mspeak.abundance
 69
 70        total_percent = (i / (i + j)) * 100
 71        total_relative_abundance = (assign_abun / (not_assign_abun + assign_abun)) * 100
 72        if report_error:
 73            rms_error = sqrt(mean(array(error) ** 2))
 74            if verbose and not mute_output:
 75                print(
 76                    "%i assigned peaks and %i unassigned peaks, total  = %.2f %%, relative abundance = %.2f %%, RMS error (best candidate) (ppm) = %.3f"
 77                    % (i, j, total_percent, total_relative_abundance, rms_error)
 78                )
 79            return i, j, total_percent, total_relative_abundance, rms_error
 80
 81        else:
 82            if verbose and not mute_output:
 83                print(
 84                    "%i assigned peaks and %i unassigned peaks , total  = %.2f %%, relative abundance = %.2f %%"
 85                    % (
 86                        i,
 87                        j,
 88                        total_percent,
 89                        total_relative_abundance,
 90                    )
 91                )
 92            return i, j, total_percent, total_relative_abundance
 93
 94    def resolving_power_calc(self, B: float, T: float):
 95        """Calculate the theoretical resolving power
 96
 97        Calls on the MSPeak object function to calculate the resolving power of a peak, this calcs for all peaks in a spectrum.
 98
 99        Parameters
100        -----------
101        T : float
102            transient time
103        B : float
104            Magnetic Filed Strength (Tesla)
105
106        References
107        ----------
108        1. Marshall et al. (Mass Spectrom Rev. 1998 Jan-Feb;17(1):1-35.)
109                DOI: 10.1002/(SICI)1098-2787(1998)17:1<1::AID-MAS1>3.0.CO;2-K
110
111        """
112
113        self.check_mspeaks()
114        return array([mspeak.resolving_power_calc(B, T) for mspeak in self.mspeaks])
115
116    def _f_to_mz(self):
117        """Ledford equation for converting frequency(Hz) to m/z
118
119        Returns
120        ----------
121        mz_domain : numpy array
122            m/z domain after conversion from frequency
123        """
124        Aterm, Bterm, Cterm = self.Aterm, self.Bterm, self.Cterm
125        # Check if the Bterm of Ledford equation scales with the ICR trap voltage or not then Bterm = Bterm*trap_voltage
126
127        if Cterm == 0:
128            if Bterm == 0:
129                # uncalibrated data
130                mz_domain = Aterm / self.freq_exp_profile
131
132            else:
133                mz_domain = (Aterm / (self.freq_exp_profile)) + (
134                    Bterm / power((self.freq_exp_profile), 2)
135                )
136
137        # @will I need you insight here, not sure what is the inverted ledford equation that Bruker refers to
138        else:
139            mz_domain = (
140                (Aterm / self.freq_exp_profile)
141                + (Bterm / power(self.freq_exp_profile, 2))
142                + Cterm
143            )
144
145        return mz_domain
146
147    def _f_to_mz_bruker(self):
148        """Frequency to m/z conversion (Bruker)
149        Bruker equations for converting frequency (Hz) to m/z,
150        nOmega acquisition is not yet implemented here.
151        However, nOmega should work for commerical Bruker 2xR systems as A Term is automatically defined for 2X or 1X by the instrument
152
153
154        Returns
155        ----------
156        numpy.array(float)
157            m/z domain after conversion from frequency
158        """
159        Aterm, Bterm, Cterm = self.Aterm, self.Bterm, self.Cterm
160        # Check if the Bterm of Ledford equation scales with the ICR trap voltage or not then Bterm = Bterm*trap_voltage
161        if Cterm == 0:
162            if Bterm == 0:
163                # uncalibrated data
164                return Aterm / self.freq_exp_profile
165
166            else:
167                # calc2
168                return Aterm / (self.freq_exp_profile + Bterm)
169
170        # @will I need you insight here, not sure what is the inverted ledford equation that Bruker refers to
171        else:
172            diff = Aterm * Aterm
173
174            # this sign(diff + 4) changes on older aquistion software
175            diff = diff + 4 * Cterm * (self.freq_exp_profile - Bterm)
176            diff = sqrt(diff)
177            diff = -Aterm + diff
178            # calc3
179            return (2 * Cterm) / diff
180            return diff / 2 * (self.freq_exp_profile - Bterm)
181
182    def number_average_molecular_weight(self, profile: bool = False):
183        """Average molecular weight calculation
184
185        Parameters
186        ----------
187        profile : bool, optional
188            is data profile or centroid mode. The default is False (e.g. Centroid data)
189
190        Returns
191        -------
192        float
193            The average molecular weight.
194
195        """
196        # mode is profile or centroid data
197        if profile:
198            a = multiply(self.mz_exp_profile, self.abundance_profile)
199            b = self.abundance_profile
200            return a.sum() / b.sum()
201
202        else:
203            return sum(self.mz_exp * self.abundance) / sum(self.abundance)
204
205    def weight_average_molecular_weight(self, profile: bool = False):
206        """
207        Weighted Average molecular weight calculation
208
209
210        Returns
211        -------
212        float
213            The weight average molecular weight.
214
215        """
216
217        # implement from MassSpectralPeaks objs
218
219        if profile:
220            a = multiply(power(self.mz_exp_profile, 2), self.abundance_profile)
221            b = self.mz_exp_profile * self.abundance_profile
222            return a.sum() / b.sum()
223
224        else:
225            return sum(power(self.mz_exp, 2) * self.abundance) / sum(
226                self.mz_exp * self.abundance
227            )