corems.mass_spectrum.calc.CalibrationCalc

  1import numpy as np
  2
  3
  4class FreqDomain_Calibration:
  5    """Frequency Domain Calibration class for mass spectrum.
  6
  7    Parameters
  8    ----------
  9    mass_spectrum : MassSpectrum
 10        The mass spectrum object.
 11    selected_mass_peaks : list
 12        List of selected mass peaks.
 13    include_isotopologue : bool, optional
 14        Flag to include isotopologues, by default False.
 15
 16    Attributes
 17    ----------
 18    mz_exp : ndarray
 19        Array of experimental m/z values.
 20    mz_calc : ndarray
 21        Array of calculated m/z values.
 22    freq_exp : ndarray
 23        Array of experimental frequencies.
 24    mass_spectrum : MassSpectrum
 25        The mass spectrum object.
 26    freq_exp_ms : ndarray
 27        Array of experimental frequencies for mass spectrum.
 28
 29    Methods
 30    -------
 31    * recal_mass_spec(mz_domain, Aterm, Bterm, Cterm).
 32        Recalibrate the mass spectrum with the given parameters.
 33    * linear().
 34        Perform linear calibration.
 35    * quadratic(iteration=False).
 36        Perform quadratic calibration.
 37    * ledford_calibration(iteration=False).
 38        Perform Ledford calibration.
 39    * step_fit(steps=4).
 40        Perform step fit calibration.
 41
 42    """
 43
 44    def __init__(self, mass_spectrum, selected_mass_peaks, include_isotopologue=False):
 45        self.selected_mspeaks = selected_mass_peaks
 46        error = list()
 47        freq_exp = list()
 48        mz_calc = list()
 49        mz_exp = list()
 50
 51        for mspeak in selected_mass_peaks:
 52            if not include_isotopologue:
 53                molecular_formulas = [
 54                    formula for formula in mspeak if not formula.is_isotopologue
 55                ]
 56            else:
 57                molecular_formulas = mspeak
 58
 59            for molecular_formula in molecular_formulas:
 60                freq_exp.append(mspeak.freq_exp)
 61                error.append(molecular_formula.mz_error)
 62                mz_calc.append(molecular_formula.mz_calc)
 63                mz_exp.append(mspeak.mz_exp)
 64
 65        self.mz_exp = np.array(mz_exp)
 66        self.mz_calc = np.array(mz_calc)
 67        self.freq_exp = np.array(freq_exp)
 68        self.mass_spectrum = mass_spectrum
 69        self.freq_exp_ms = np.array([mspeak.freq_exp for mspeak in mass_spectrum])
 70
 71    def recal_mass_spec(self, mz_domain, Aterm, Bterm, Cterm):
 72        """Recalibrate the mass spectrum with the given parameters.
 73
 74        Parameters
 75        ----------
 76        mz_domain : ndarray
 77            Array of m/z values for recalibration.
 78        Aterm : float
 79            Aterm parameter for recalibration.
 80        Bterm : float
 81            Bterm parameter for recalibration.
 82        Cterm : float
 83            Cterm parameter for recalibration.
 84
 85        """
 86        self.mass_spectrum._calibration_terms = (Aterm, Bterm, 0)
 87        self.mass_spectrum.mz_cal = mz_domain
 88
 89    def linear(self):
 90        """Perform linear calibration."""
 91        matrix = np.vstack([1 / self.freq_exp, np.ones(len(self.freq_exp))]).T
 92        Aterm, Bterm = np.linalg.lstsq(matrix, self.mz_calc, rcond=None)[0]
 93        if self.mass_spectrum.parameters.mass_spectrum.verbose_processing:
 94            print("%.2f Aterm,  %.2f Bterm" % (Aterm, Bterm))
 95            print("Linear Calibration %.2f Aterm,  %.2f Bterm " % (Aterm, Bterm))
 96        mz_domain = (Aterm / self.freq_exp_ms) + Bterm
 97        self.recal_mass_spec(mz_domain, Aterm, Bterm, 0)
 98
 99    def quadratic(self, iteration: bool = False):
100        """Perform quadratic calibration.
101
102        Parameters
103        ----------
104        iteration : bool, optional
105            Flag to perform iterative calibration, by default False.
106
107        """
108        mz_calc = self.mz_calc
109        freq_exp = self.freq_exp
110        mz_exp = self.mz_exp
111
112        error = ((mz_exp - mz_calc) / mz_calc) * 1000000
113        last_rms = np.sqrt(np.mean(error**2))
114        while True:
115            matrix = np.vstack(
116                [1 / freq_exp, 1 / np.power(freq_exp, 2), np.ones(len(freq_exp))]
117            ).T
118            Aterm, Bterm, Cterm = np.linalg.lstsq(matrix, self.mz_calc, rcond=None)[0]
119            mz_exp = (Aterm / (freq_exp)) + (Bterm / np.power((freq_exp), 2)) + Cterm
120            error = ((mz_exp - mz_calc) / mz_calc) * 1000000
121            rms = np.sqrt(np.mean(error**2))
122            std = np.std(error)
123            if self.mass_spectrum.parameters.mass_spectrum.verbose_processing:
124                print("%.2f Aterm,  %.2f Bterm" % (Aterm, Bterm))
125                print(
126                    "Quadratic Calibration %.2f RMS,  %.2f std,  %.2f Aterm,  %.2f Bterm "
127                    % (rms, std, Aterm, Bterm)
128                )
129            if rms < last_rms:
130                last_rms = rms
131                freq_exp = (
132                    Aterm
133                    + np.sqrt(np.power(-Aterm, 2) - (4 * Cterm * (mz_exp - Bterm)))
134                ) / (2 * mz_exp)
135
136                mz_domain = (
137                    (Aterm / (self.freq_exp_ms))
138                    + (Bterm / np.power((self.freq_exp_ms), 2))
139                    + Cterm
140                )
141                self.recal_mass_spec(mz_domain, Aterm, Bterm, Cterm)
142                if not iteration:
143                    break
144            else:
145                break
146
147    def ledford_calibration(self, iteration: bool = False):
148        """Perform Ledford calibration.
149
150        Parameters
151        ----------
152        iteration : bool, optional
153            Flag to perform iterative calibration, by default False.
154
155        """
156        mz_calc = self.mz_calc
157        freq_exp = self.freq_exp
158        mz_exp = self.mz_exp
159
160        error = ((mz_exp - self.mz_calc) / self.mz_calc) * 1000000
161        last_rms = np.sqrt(np.mean(error**2))
162        while True:
163            matrix = np.vstack([1 / freq_exp, 1 / np.power(freq_exp, 2)]).T
164            Aterm, Bterm = np.linalg.lstsq(matrix, self.mz_calc, rcond=None)[0]
165
166            mz_exp = (Aterm / (freq_exp)) + (Bterm / np.power((freq_exp), 2))
167            error = ((mz_exp - mz_calc) / mz_calc) * 1000000
168            rms = np.sqrt(np.mean(error**2))
169            std = np.std(error)
170            if self.mass_spectrum.parameters.mass_spectrum.verbose_processing:
171                print("%.2f Aterm,  %.2f Bterm" % (Aterm, Bterm))
172                print(
173                    "Ledford Calibration %.2f RMS,  %.2f std,  %.2f Aterm,  %.2f Bterm "
174                    % (rms, std, Aterm, Bterm)
175                )
176            if rms < last_rms:
177                last_rms = rms
178                freq_exp = (
179                    Aterm + np.sqrt(np.power(-Aterm, 2) - (4 * mz_exp - Bterm))
180                ) / (2 * mz_exp)
181                mz_domain = (Aterm / (self.freq_exp_ms)) + (
182                    Bterm / np.power((self.freq_exp_ms), 2)
183                )
184                self.recal_mass_spec(mz_domain, Aterm, Bterm, 0)
185                if not iteration:
186                    break
187            else:
188                break
189
190    def step_fit(self, steps: int = 4):
191        """Perform step fit calibration.
192
193        Parameters
194        ----------
195        steps : int, optional
196            Number of steps for step fit calibration, by default 4.
197
198        """
199
200        def f_to_mz(f, A, B, C, a):
201            return (A / f) + (B / np.power(f, 2)) + (C * a / np.power(f, 2))
202
203        def mz_to_f(m, A, B, C):
204            return -A - m / B
205
206        tuple_indexes = [
207            (i, i + steps) for i in range(0, len(self.selected_mspeaks) - steps, steps)
208        ]
209
210        for current_index, tuple_index in enumerate(tuple_indexes):
211            mspeak_ii, mspeak_fi = tuple_index
212            freq_exp = list()
213            mz_calc = list()
214            mz_exp = list()
215            abu = list()
216
217            for i in range(mspeak_ii, mspeak_fi + 1):
218                best_formula = self.selected_mspeaks[i].best_molecular_formula_candidate
219
220                freq_exp.append(self.selected_mspeaks[i].freq_exp)
221                mz_calc.append(best_formula.mz_calc)
222                mz_exp.append(self.selected_mspeaks[i].mz_exp)
223                abu.append(self.selected_mspeaks[i].abundance)
224
225            freq_exp = np.array(freq_exp)
226            mz_calc = np.array(mz_calc)
227            mz_exp = np.array(mz_exp)
228            abu = np.array(abu)
229
230            if current_index == len(tuple_indexes) - 1:
231                ms_peaks_indexes = (self.selected_mspeaks[mspeak_ii].index, 0)
232
233            elif current_index == 0:
234                ms_peaks_indexes = (
235                    len(self.mass_spectrum) - 1,
236                    self.selected_mspeaks[mspeak_fi].index - 1,
237                )
238            else:
239                ms_peaks_indexes = (
240                    self.selected_mspeaks[mspeak_ii].index,
241                    self.selected_mspeaks[mspeak_fi].index - 1,
242                )
243
244            final_index, start_index = ms_peaks_indexes
245
246            matrix = np.vstack([1 / freq_exp, 1 / np.power(freq_exp, 2)]).T
247            A, B = np.linalg.lstsq(matrix, mz_calc, rcond=None)[0]
248            C = 0
249
250            for mspeak in self.mass_spectrum[start_index:final_index]:
251                mspeak.mz_cal = f_to_mz(mspeak.freq_exp, A, B, C, 0)
252
253        self.mass_spectrum.is_calibrated = True