corems.mass_spectrum.calc.AutoRecalibration

  1# -*- coding: utf-8 -*-
  2"""
  3Created on March 23 2023
  4
  5@author: Will Kew
  6
  7Modules for automatic mass internal recalibration
  8"""
  9
 10from lmfit.models import GaussianModel
 11import seaborn as sns
 12import pandas as pd
 13import numpy as np
 14import matplotlib.pyplot as plt
 15from corems.molecular_id.search.molecularFormulaSearch import SearchMolecularFormulas
 16import copy
 17
 18
 19class HighResRecalibration:
 20    """
 21    This class is designed for high resolution (FTICR, Orbitrap) data of complex mixture, e.g. Organic matter
 22
 23    The tool first does a broad mass range search for the most commonly expected ion type (i.e. CHO, deprotonated - for negative ESI)
 24    And then the assigned data mass error distribution is searched, with a gaussian fit to the most prominent range.
 25    This tool works when the data are of sufficient quality, and not outwith the typical expected range of the mass analyzer
 26    It presumes the mean error is out by 0-several ppm, but that the spread of error values is modest (<2ppm)
 27
 28    Parameters
 29    ----------
 30    mass_spectrum : MassSpectrum
 31        CoreMS mass spectrum object
 32    plot : bool, optional
 33        Whether to plot the error distribution. The default is False.
 34    docker : bool, optional
 35        Whether to use the docker database. The default is True. If not, it uses a dynamically generated sqlite database.
 36    ppmFWHMprior : float, optional
 37        The FWHM of the prior distribution (ppm). The default is 3.
 38    ppmRangeprior : float, optional
 39        The range of the prior distribution (ppm). The default is 15.
 40
 41    Methods
 42    --------
 43    * determine_error_boundaries(). Determine the error boundaries for recalibration space.
 44
 45    Notes
 46    -----
 47    This initialisation function creates a copy of the MassSpectrum object to avoid over-writing assignments.
 48    Possible future task is to make the base class copyable.
 49
 50    """
 51
 52    def __init__(
 53        self,
 54        mass_spectrum,
 55        plot: bool = False,
 56        docker: bool = True,
 57        ppmFWHMprior: float = 3,
 58        ppmRangeprior: float = 15,
 59    ):
 60        self.mass_spectrum = copy.deepcopy(mass_spectrum)
 61        self.plot = plot
 62        self.docker = docker
 63        self.ppmFWHMprior = ppmFWHMprior
 64        self.ppmRangeprior = ppmRangeprior
 65
 66    def set_uncal_settings(self):
 67        """Set uncalibrated formula search settings
 68
 69        This function serves the uncalibrated data (hence broad error tolerance)
 70        It only allows CHO formula in deprotonated ion type- as most common for SRFA ESI negative mode
 71
 72        This will not work for positive mode data, or for other ion types, or other expected elemental searches.
 73
 74        """
 75        # TODO rework this.
 76
 77        if self.docker:
 78            self.mass_spectrum.molecular_search_settings.url_database = "postgresql+psycopg2://coremsappdb:coremsapppnnl@localhost:5432/coremsapp"
 79        else:
 80            self.mass_spectrum.molecular_search_settings.url_database = None
 81        self.mass_spectrum.molecular_search_settings.error_method = None
 82        self.mass_spectrum.molecular_search_settings.score_method = "prob_score"
 83
 84        self.mass_spectrum.molecular_search_settings.min_ppm_error = (
 85            -1 * self.ppmRangeprior / 2
 86        )  # -7.5
 87        self.mass_spectrum.molecular_search_settings.max_ppm_error = (
 88            self.ppmRangeprior / 2
 89        )  # 7.5
 90
 91        self.mass_spectrum.molecular_search_settings.min_dbe = 0
 92        self.mass_spectrum.molecular_search_settings.max_dbe = 50
 93
 94        self.mass_spectrum.molecular_search_settings.use_isotopologue_filter = False
 95        self.mass_spectrum.molecular_search_settings.min_abun_error = -30
 96        self.mass_spectrum.molecular_search_settings.max_abun_error = 70
 97
 98        self.mass_spectrum.molecular_search_settings.use_min_peaks_filter = True
 99        self.mass_spectrum.molecular_search_settings.min_peaks_per_class = (
100            10  # default is 15
101        )
102
103        self.mass_spectrum.molecular_search_settings.usedAtoms["C"] = (1, 90)
104        self.mass_spectrum.molecular_search_settings.usedAtoms["H"] = (4, 200)
105        self.mass_spectrum.molecular_search_settings.usedAtoms["O"] = (1, 23)
106        self.mass_spectrum.molecular_search_settings.usedAtoms["N"] = (0, 0)
107        self.mass_spectrum.molecular_search_settings.usedAtoms["S"] = (0, 0)
108        self.mass_spectrum.molecular_search_settings.usedAtoms["P"] = (0, 0)
109
110        self.mass_spectrum.molecular_search_settings.isProtonated = True
111        self.mass_spectrum.molecular_search_settings.isRadical = False
112        self.mass_spectrum.molecular_search_settings.isAdduct = False
113
114    def positive_search_settings(self):
115        """Set the positive mode elemental search settings"""
116        self.mass_spectrum.molecular_search_settings.isProtonated = False
117        self.mass_spectrum.molecular_search_settings.isAdduct = True
118        self.mass_spectrum.molecular_search_settings.adduct_atoms_pos = ["Na"]
119
120    @staticmethod
121    def get_error_range(
122        errors: list, ppmFWHMprior: float = 3, plot_logic: bool = False
123    ):
124        """Get the error range from the error distribution
125
126        Using lmfit and seaborn kdeplot to extract the error range from the error distribution of assigned species.
127
128        Parameters
129        ----------
130        errors : list
131            list of the errors of the assigned species (ppm)
132        ppmFWHMprior : float, optional
133            The FWHM of the prior distribution (ppm). The default is 3.
134        plot_logic : bool, optional
135            Whether to plot the error distribution. The default is False.
136
137        Returns
138        -------
139        mean_error : float
140            mean mass error of the Gaussian distribution (ppm)
141        fwhm_error : float
142            full width half max of the gaussian error distribution (ppm)
143        ppm_thresh : list
144            recommended thresholds for the recalibration parameters (ppm)
145            Consists of [mean_error-fwhm_error,mean_error+fwhm_error]
146
147        """
148        # Create an isolated figure so stale lines from other tests/plots
149        # do not pollute get_lines()[0] (matplotlib global state issue)
150        fig, ax = plt.subplots()
151        kde = sns.kdeplot(errors, ax=ax)
152
153        kde_data = ax.get_lines()[0].get_data()
154
155        tmpdf = pd.Series(index=kde_data[0], data=kde_data[1])
156        kde_apex_ppm = tmpdf.idxmax()
157        kde_apex_val = tmpdf.max()
158
159        plt.close(fig)
160        plt.close("all")
161
162        lmmodel = GaussianModel()
163        lmpars = lmmodel.guess(kde_data[1], x=kde_data[0])
164        lmpars["sigma"].value = 2.3548 / ppmFWHMprior
165        lmpars["center"].value = kde_apex_ppm
166        lmpars["amplitude"].value = kde_apex_val
167        lmout = lmmodel.fit(kde_data[1], lmpars, x=kde_data[0])
168
169        if plot_logic:
170            fig, ax = plt.subplots(figsize=(8, 4))
171            lmout.plot_fit(
172                ax=ax, data_kws={"color": "tab:blue"}, fit_kws={"color": "tab:red"}
173            )
174            ax.set_xlabel("$m/z$ Error (ppm)")
175            ax.set_ylabel("Density")
176            plt.legend(facecolor="white", framealpha=0)
177
178        mean_error = lmout.best_values["center"]
179        std_error = lmout.best_values["sigma"]
180        # FWHM from Sigma = approx. 2.355*sigma
181        # fwhm_error = 2*np.sqrt(2*np.log(2))*std_error
182        fwhm_error = std_error * np.sqrt(8 * np.log(2))
183
184        ppm_thresh = [mean_error - fwhm_error, mean_error + fwhm_error]
185        return mean_error, fwhm_error, ppm_thresh
186
187    def determine_error_boundaries(self):
188        """Determine the error boundaries for recalibration space
189
190        This is the main function in this class
191        Sets the Molecular Formulas search settings, performs the initial formula search
192        Converts the data to a dataframe, and gets the error range
193        Returns the error thresholds.
194
195        Returns
196        -------
197        mean_error : float
198            mean mass error of the Gaussian distribution (ppm)
199        fwhm_error : float
200            full width half max of the gaussian error distribution (ppm)
201        ppm_thresh : list
202            recommended thresholds for the recalibration parameters (ppm)
203            Consists of [mean_error-fwhm_error,mean_error+fwhm_error]
204        """
205
206        # Set the search settings
207        self.set_uncal_settings()
208
209        # Set the positive mode settings
210        # To do - have user defineable settings?
211        if self.mass_spectrum.polarity == 1:
212            self.positive_search_settings()
213
214        # Search MFs
215        SearchMolecularFormulas(
216            self.mass_spectrum, first_hit=True
217        ).run_worker_mass_spectrum()
218
219        # Exporting to a DF is ~30x slower than just getting the errors, so this is fast.
220        errors = []
221        for mspeak in self.mass_spectrum.mspeaks:
222            if len(mspeak.molecular_formulas) > 0:
223                errors.append(mspeak.best_molecular_formula_candidate.mz_error)
224
225        # If there are NO assignments, it'll fail on the next step. Need to check for that
226        nassign = len(errors)
227        # Here we say at least 5 features assigned are needed - it probably should be greater, but we are just trying to stop it breaking the code
228        # We want to make sure the spectrum is capture in the database though - so we return the stats entries (0 assignments) and the number of assignments
229        if nassign < 5:
230            if self.mass_spectrum.parameters.mass_spectrum.verbose_processing:
231                print("fewer than 5 peaks assigned, cannot determine error range")
232            return np.nan, np.nan, [np.nan, np.nan]
233        else:
234            mean_error, fwhm_error, ppm_thresh = self.get_error_range(
235                errors, self.ppmFWHMprior, self.plot
236            )
237            return mean_error, fwhm_error, ppm_thresh