corems.molecular_id.search.lcms_spectral_search

  1import re
  2
  3import numpy as np
  4
  5from corems.molecular_id.factory.spectrum_search_results import SpectrumSearchResults
  6
  7
  8class LCMSSpectralSearch:
  9    """
 10    Methods for searching LCMS spectra.
 11
 12    This class is designed to be a mixin class for the :obj:`~corems.mass_spectra.factory.lc_class.LCMSBase` class.
 13
 14    """
 15
 16    @staticmethod
 17    def get_more_match_quals(
 18        query_mz_arr, lib_entry, mz_tol_da=0.1, include_fragment_types=False
 19    ):
 20        """
 21        Return additional match qualities between query and library entry.
 22
 23        Parameters
 24        ----------
 25        query_mz_arr : np.array
 26            Array of query spectrum. Shape (N, 2), with m/z in the first column
 27            and abundance in the second.
 28        lib_entry : dict
 29            Library spectrum entry, with 'mz' key containing the spectrum in
 30            the format (mz, abundance),(mz, abundance), i.e. from MetabRef.
 31        mz_tol_da : float, optional
 32            Tolerance in Da for matching peaks (in MS2). Default is 0.1.
 33        include_fragment_types : bool, optional
 34            If True, include fragment type comparisons in output.
 35            Defaults to False.
 36
 37        Returns
 38        -------
 39        tuple
 40            Tuple of (query_in_lib, query_in_lib_fract, lib_in_query, lib_in_query_fract, query_frags, lib_frags, lib_precursor_mz).
 41
 42        Notes
 43        -----
 44        query_in_lib : int
 45            Number of peaks in query that are present in the library entry (within mz_tol_da).
 46        query_in_lib_fract : float
 47            Fraction of peaks in query that are present in the library entry (within mz_tol_da).
 48        lib_in_query : int
 49            Number of peaks in the library entry that are present in the query (within mz_tol_da).
 50        lib_in_query_fract : float
 51            Fraction of peaks in the library entry that are present in the query (within mz_tol_da).
 52        query_frags : list
 53            List of unique fragment types present in the query, generally 'MLF' or 'LSF' or both.
 54        lib_frags : list
 55            List of unique fragment types present in the library entry, generally 'MLF' or 'LSF' or both.
 56
 57        Raises
 58        ------
 59        ValueError
 60            If library entry does not have 'fragment_types' key and include_fragment_types is True.
 61
 62        """
 63
 64        if "mz" in lib_entry.keys():
 65            # Get the original mz values from the library entry
 66            lib_mzs = np.array(
 67                re.findall(r"\(([^,]+),([^)]+)\)", lib_entry["mz"]), dtype=float
 68            ).reshape(-1, 2)[:, 0]
 69        elif "peaks" in lib_entry.keys() and lib_entry["peaks"] is not None:
 70            lib_mzs = lib_entry["peaks"][:, 0]
 71
 72        # Get count and fraction of peaks in query that are in lib entry
 73        query_in_lib = 0
 74        for peak in query_mz_arr:
 75            if np.any(np.isclose(lib_mzs, peak, atol=mz_tol_da)):
 76                query_in_lib += 1
 77        query_in_lib_fract = query_in_lib / len(query_mz_arr)
 78
 79        # Get count and fraction of peaks in lib that are in query
 80        lib_in_query = 0
 81        for peak in lib_mzs:
 82            if np.any(np.isclose(query_mz_arr, peak, atol=mz_tol_da)):
 83                lib_in_query += 1
 84        lib_in_query_fract = lib_in_query / len(lib_mzs)
 85
 86        if include_fragment_types:
 87            # Check that fragment types are present in the library entry
 88            if "fragment_types" not in lib_entry.keys():
 89                raise ValueError(
 90                    "Flash entropy library entry must have 'fragment_types' key to include fragment types in output."
 91                )
 92
 93            # Get types of fragments in the lib entry and convert it to a list on ", "
 94            lib_frags = [x.strip() for x in lib_entry["fragment_types"].split(",")]
 95            # make list of the fragment types that are present in the query spectrum
 96            lib_in_query_ids = list(
 97                set(
 98                    [
 99                        ind
100                        for ind, x in enumerate(lib_mzs)
101                        if len(np.where(np.isclose(query_mz_arr, x, atol=mz_tol_da))[0])
102                        > 0
103                    ]
104                )
105            )
106            query_frags = list(set([lib_frags[x] for x in lib_in_query_ids]))
107            lib_frags = list(set(lib_frags))
108
109        else:
110            query_frags = None
111            lib_frags = None
112
113        return (
114            query_in_lib,
115            query_in_lib_fract,
116            lib_in_query,
117            lib_in_query_fract,
118            query_frags,
119            lib_frags,
120        )
121
122    def fe_search(
123        self,
124        scan_list,
125        fe_lib,
126        precursor_mz_list=[],
127        use_mass_features=True,
128        peak_sep_da=None,
129        get_additional_metrics=True,
130        accumulate_results=False,
131        ms_params_key: str = "ms2",
132    ):
133        """
134        Search LC-MS MS2 spectra with FlashEntropy against a pre-built library.
135
136        Annotation thresholds and default peak-separation for cleaning come from
137        :class:`~corems.encapsulation.factory.processingSetting.SpectralSimilaritySearchSettings`
138        on this object's parameter tree
139        (``parameters.mass_spectrum[ms_params_key].spectral_similarity_search``).
140        Build ``fe_lib`` with the **same** settings so library index tolerances
141        match search cleaning (see
142        :meth:`~corems.molecular_id.search.database_interfaces.MSPInterface.get_metabolomics_spectra_library`
143        and
144        :meth:`~corems.molecular_id.search.database_interfaces.LCLipidLibraryInterface.get_lipid_library`).
145
146        Parameters
147        ----------
148        scan_list : list
149            List of scan numbers to search (must already be loaded in ``self._ms``).
150        fe_lib : :obj:`~ms_entropy.FlashEntropySearch`
151            FlashEntropy search instance. Prefer constructing it with
152            ``settings=settings_from_lcms(self)`` (or the same profile key as
153            *ms_params_key*) so FE build knobs match this search.
154        precursor_mz_list : list, optional
155            List of precursor m/z values to search, by default [], which implies
156            matched with mass features; to enable this use_mass_features must be True.
157        use_mass_features : bool, optional
158            If True, use mass features to get precursor m/z values, by default True.
159            If True, will add search results to mass features' ms2_similarity_results attribute.
160        peak_sep_da : float, optional
161            Minimum separation between m/z peaks spectra in Da. This needs match the
162            approximate resolution of the search spectra and the FlashEntropySearch
163            instance. If None (default), uses
164            ``parameters.mass_spectrum[ms_params_key].spectral_similarity_search.resolved_peak_sep_da``
165            (typically ``2 * max_ms2_tolerance_in_da`` from that settings object).
166        get_additional_metrics : bool, optional
167            If True, get additional metrics from FlashEntropy search, by default True.
168        accumulate_results : bool, optional
169            If True, accumulate results with existing spectral_search_results instead of
170            replacing them. This allows searching the same scans with multiple libraries
171            without overwriting previous results, by default False.
172        ms_params_key : str, optional
173            Key in ``parameters.mass_spectrum`` whose ``spectral_similarity_search``
174            settings are used (default ``"ms2"``). Use another key (e.g.
175            ``"ms2_cid"``) when multiple MS2 parameter profiles are configured
176            for different scan classes.
177
178        Returns
179        -------
180        None
181            Adds results to ``self.spectral_search_results`` and associates them
182            with mass features in ``self.mass_features`` when applicable.
183
184        Notes
185        -----
186        From the settings profile this method reads:
187
188        - ``ms2_min_fe_score`` — minimum entropy score to keep a hit
189        - ``include_fragment_types`` — lipid-style fragment-type metrics
190        - default peak separation when *peak_sep_da* is None
191
192        FlashEntropy **library** build fields (``max_ms2_tolerance_in_da``, etc.)
193        should already have been applied when constructing *fe_lib*. Configure
194        them on the same ``spectral_similarity_search`` instance before library
195        generation.
196
197        For an :class:`~corems.mass_spectra.factory.lc_class.LCMSCollection`,
198        sample parameters are equal at construction; build one shared library
199        with ``settings_from_lcms_collection(collection)``, then call
200        ``fe_search`` on each sample (or via collection pipeline operations).
201
202        Examples
203        --------
204        Single LCMSBase object (library + search share profile ``"ms2"``)::
205
206            from corems.encapsulation.factory.parameters import settings_from_lcms
207
208            settings = settings_from_lcms(lcms_obj, profile="ms2")
209            fe_lib, _meta = msp.get_metabolomics_spectra_library(
210                polarity="positive",
211                format="flashentropy",
212                settings=settings,
213            )
214            lcms_obj.fe_search(scan_list=ms2_scans, fe_lib=fe_lib)
215            # peak_sep / score gate from mass_spectrum["ms2"].spectral_similarity_search
216
217        Alternate MS2 profile (e.g. CID bag)::
218
219            lcms_obj.fe_search(
220                scan_list=cid_scans,
221                fe_lib=fe_lib_cid,
222                ms_params_key="ms2_cid",
223            )
224        """
225        # Annotation knobs from nested SpectralSimilaritySearchSettings (default profile "ms2")
226        if ms_params_key not in self.parameters.mass_spectrum:
227            raise KeyError(
228                f"ms_params_key={ms_params_key!r} not in parameters.mass_spectrum "
229                f"(keys={list(self.parameters.mass_spectrum)})"
230            )
231        ms2_p = self.parameters.mass_spectrum[ms_params_key].spectral_similarity_search
232        include_fragment_types = ms2_p.include_fragment_types
233        min_match_score = ms2_p.ms2_min_fe_score
234        if peak_sep_da is None:
235            peak_sep_da = ms2_p.resolved_peak_sep_da
236
237        # If precursor_mz_list is empty and use_mass_features is True, get precursor m/z values from mass features for each scan in scan_list
238        if use_mass_features and len(precursor_mz_list) == 0:
239            precursor_mz_list = []
240            for scan in scan_list:
241                mf_ids = [
242                    key
243                    for key, value in self.mass_features.items()
244                    if scan in value.ms2_mass_spectra
245                ]
246                precursor_mz = [
247                    value.mz
248                    for key, value in self.mass_features.items()
249                    if key in mf_ids
250                ]
251                precursor_mz_list.append(precursor_mz)
252
253        # Check that precursor_mz_list same length as scan_list, if not, raise error
254        if len(precursor_mz_list) != len(scan_list):
255            raise ValueError("Length of precursor_mz_list is not equal to scan_list.")
256
257        # Loop through each query spectrum / precursor match and save ids of db spectrum that are decent matches
258        overall_results_dict = {}
259        for i in np.arange(len(scan_list)):
260            scan_oi = scan_list[i]
261            if len(self._ms[scan_oi].mspeaks) > 0:
262                precursor_mzs = precursor_mz_list[i]
263                overall_results_dict[scan_oi] = {}
264                for precursor_mz in precursor_mzs:
265                    query_spectrum = fe_lib.clean_spectrum_for_search(
266                        precursor_mz=precursor_mz,
267                        peaks=np.vstack(
268                            (self._ms[scan_oi].mz_exp, self._ms[scan_oi].abundance)
269                        ).T,
270                        precursor_ions_removal_da=None,
271                        noise_threshold=self._ms[
272                            scan_oi
273                        ].parameters.mass_spectrum.noise_threshold_min_relative_abundance
274                        / 100,
275                        min_ms2_difference_in_da=peak_sep_da,
276                    )
277                    search_results = fe_lib.search(
278                        precursor_mz=precursor_mz,
279                        peaks=query_spectrum,
280                        ms1_tolerance_in_da=self.parameters.mass_spectrum[
281                            "ms1"
282                        ].molecular_search.max_ppm_error
283                        * 10**-6
284                        * precursor_mz,
285                        ms2_tolerance_in_da=peak_sep_da * 0.5,
286                        method={"identity"},
287                        precursor_ions_removal_da=None,
288                        noise_threshold=self._ms[
289                            scan_oi
290                        ].parameters.mass_spectrum.noise_threshold_min_relative_abundance
291                        / 100,
292                        target="cpu",
293                    )["identity_search"]
294                    match_inds = np.where(search_results > min_match_score)[0]
295
296                    # If any decent matches are found, add them to the results dictionary
297                    if len(match_inds) > 0:
298                        match_scores = search_results[match_inds]
299                        ref_ms_ids = [fe_lib[x]["id"] for x in match_inds]
300                        ref_mol_ids = [
301                            fe_lib[x]["molecular_data_id"] for x in match_inds
302                        ]
303                        ref_precursor_mzs = [
304                            fe_lib[x]["precursor_mz"] for x in match_inds
305                        ]
306                        ion_types = [fe_lib[x]["ion_type"] for x in match_inds]
307                        overall_results_dict[scan_oi][precursor_mz] = {
308                            "ref_mol_id": ref_mol_ids,
309                            "ref_ms_id": ref_ms_ids,
310                            "ref_precursor_mz": ref_precursor_mzs,
311                            "precursor_mz_error_ppm": [
312                                (precursor_mz - x) / precursor_mz * 10**6
313                                for x in ref_precursor_mzs
314                            ],
315                            "entropy_similarity": match_scores,
316                            "ref_ion_type": ion_types,
317                        }
318                        # Add database name, if present
319                        db_name = [fe_lib[x].get("database_name") for x in match_inds]
320                        if db_name is not None:
321                            overall_results_dict[scan_oi][precursor_mz].update(
322                                {"database_name": db_name}
323                            )
324                        if get_additional_metrics:
325                            more_match_quals = [
326                                self.get_more_match_quals(
327                                    self._ms[scan_oi].mz_exp,
328                                    fe_lib[x],
329                                    mz_tol_da=peak_sep_da,
330                                    include_fragment_types=include_fragment_types,
331                                )
332                                for x in match_inds
333                            ]
334                            overall_results_dict[scan_oi][precursor_mz].update(
335                                {
336                                    "query_mz_in_ref_n": [
337                                        x[0] for x in more_match_quals
338                                    ],
339                                    "query_mz_in_ref_fract": [
340                                        x[1] for x in more_match_quals
341                                    ],
342                                    "ref_mz_in_query_n": [
343                                        x[2] for x in more_match_quals
344                                    ],
345                                    "ref_mz_in_query_fract": [
346                                        x[3] for x in more_match_quals
347                                    ],
348                                }
349                            )
350                            if include_fragment_types:
351                                overall_results_dict[scan_oi][precursor_mz].update(
352                                    {
353                                        "query_frag_types": [
354                                            x[4] for x in more_match_quals
355                                        ],
356                                        "ref_frag_types": [
357                                            x[5] for x in more_match_quals
358                                        ],
359                                    }
360                                )
361
362        # Drop scans with no results from dictionary
363        overall_results_dict = {k: v for k, v in overall_results_dict.items() if v}
364
365        # Cast each entry as a MS2SearchResults object
366        for scan_id in overall_results_dict.keys():
367            for precursor_mz in overall_results_dict[scan_id].keys():
368                ms2_spectrum = self._ms[scan_id]
369                ms2_search_results = overall_results_dict[scan_id][precursor_mz]
370                overall_results_dict[scan_id][precursor_mz] = SpectrumSearchResults(
371                    ms2_spectrum, precursor_mz, ms2_search_results
372                )
373
374        # Add MS2SearchResults to the existing spectral search results dictionary
375        if accumulate_results:
376            # Merge results with existing spectral_search_results
377            for scan_id, precursor_dict in overall_results_dict.items():
378                if scan_id in self.spectral_search_results:
379                    # Scan already has results, merge precursor_mz dictionaries
380                    self.spectral_search_results[scan_id].update(precursor_dict)
381                else:
382                    # New scan, add entire dictionary
383                    self.spectral_search_results[scan_id] = precursor_dict
384        else:
385            # Replace existing results (original behavior)
386            self.spectral_search_results.update(overall_results_dict)
387
388        # If there are mass features, associate the results with each mass feature
389        if len(self.mass_features) > 0:
390            # Determine which results to associate with mass features
391            if accumulate_results:
392                # When accumulating, only associate new results from this search
393                # to avoid duplicating previously associated results
394                results_to_associate = overall_results_dict
395            else:
396                # When not accumulating, clear existing associations and re-associate all results
397                for mass_feature_id in self.mass_features.keys():
398                    self.mass_features[mass_feature_id].ms2_similarity_results = []
399                results_to_associate = self.spectral_search_results
400
401            for mass_feature_id, mass_feature in self.mass_features.items():
402                scan_ids = mass_feature.ms2_scan_numbers
403                for ms2_scan_id in scan_ids:
404                    precursor_mz = mass_feature.mz
405                    try:
406                        results_to_associate[ms2_scan_id][precursor_mz]
407                    except KeyError:
408                        pass
409                    else:
410                        self.mass_features[
411                            mass_feature_id
412                        ].ms2_similarity_results.append(
413                            results_to_associate[ms2_scan_id][precursor_mz]
414                        )
class LCMSSpectralSearch:
  9class LCMSSpectralSearch:
 10    """
 11    Methods for searching LCMS spectra.
 12
 13    This class is designed to be a mixin class for the :obj:`~corems.mass_spectra.factory.lc_class.LCMSBase` class.
 14
 15    """
 16
 17    @staticmethod
 18    def get_more_match_quals(
 19        query_mz_arr, lib_entry, mz_tol_da=0.1, include_fragment_types=False
 20    ):
 21        """
 22        Return additional match qualities between query and library entry.
 23
 24        Parameters
 25        ----------
 26        query_mz_arr : np.array
 27            Array of query spectrum. Shape (N, 2), with m/z in the first column
 28            and abundance in the second.
 29        lib_entry : dict
 30            Library spectrum entry, with 'mz' key containing the spectrum in
 31            the format (mz, abundance),(mz, abundance), i.e. from MetabRef.
 32        mz_tol_da : float, optional
 33            Tolerance in Da for matching peaks (in MS2). Default is 0.1.
 34        include_fragment_types : bool, optional
 35            If True, include fragment type comparisons in output.
 36            Defaults to False.
 37
 38        Returns
 39        -------
 40        tuple
 41            Tuple of (query_in_lib, query_in_lib_fract, lib_in_query, lib_in_query_fract, query_frags, lib_frags, lib_precursor_mz).
 42
 43        Notes
 44        -----
 45        query_in_lib : int
 46            Number of peaks in query that are present in the library entry (within mz_tol_da).
 47        query_in_lib_fract : float
 48            Fraction of peaks in query that are present in the library entry (within mz_tol_da).
 49        lib_in_query : int
 50            Number of peaks in the library entry that are present in the query (within mz_tol_da).
 51        lib_in_query_fract : float
 52            Fraction of peaks in the library entry that are present in the query (within mz_tol_da).
 53        query_frags : list
 54            List of unique fragment types present in the query, generally 'MLF' or 'LSF' or both.
 55        lib_frags : list
 56            List of unique fragment types present in the library entry, generally 'MLF' or 'LSF' or both.
 57
 58        Raises
 59        ------
 60        ValueError
 61            If library entry does not have 'fragment_types' key and include_fragment_types is True.
 62
 63        """
 64
 65        if "mz" in lib_entry.keys():
 66            # Get the original mz values from the library entry
 67            lib_mzs = np.array(
 68                re.findall(r"\(([^,]+),([^)]+)\)", lib_entry["mz"]), dtype=float
 69            ).reshape(-1, 2)[:, 0]
 70        elif "peaks" in lib_entry.keys() and lib_entry["peaks"] is not None:
 71            lib_mzs = lib_entry["peaks"][:, 0]
 72
 73        # Get count and fraction of peaks in query that are in lib entry
 74        query_in_lib = 0
 75        for peak in query_mz_arr:
 76            if np.any(np.isclose(lib_mzs, peak, atol=mz_tol_da)):
 77                query_in_lib += 1
 78        query_in_lib_fract = query_in_lib / len(query_mz_arr)
 79
 80        # Get count and fraction of peaks in lib that are in query
 81        lib_in_query = 0
 82        for peak in lib_mzs:
 83            if np.any(np.isclose(query_mz_arr, peak, atol=mz_tol_da)):
 84                lib_in_query += 1
 85        lib_in_query_fract = lib_in_query / len(lib_mzs)
 86
 87        if include_fragment_types:
 88            # Check that fragment types are present in the library entry
 89            if "fragment_types" not in lib_entry.keys():
 90                raise ValueError(
 91                    "Flash entropy library entry must have 'fragment_types' key to include fragment types in output."
 92                )
 93
 94            # Get types of fragments in the lib entry and convert it to a list on ", "
 95            lib_frags = [x.strip() for x in lib_entry["fragment_types"].split(",")]
 96            # make list of the fragment types that are present in the query spectrum
 97            lib_in_query_ids = list(
 98                set(
 99                    [
100                        ind
101                        for ind, x in enumerate(lib_mzs)
102                        if len(np.where(np.isclose(query_mz_arr, x, atol=mz_tol_da))[0])
103                        > 0
104                    ]
105                )
106            )
107            query_frags = list(set([lib_frags[x] for x in lib_in_query_ids]))
108            lib_frags = list(set(lib_frags))
109
110        else:
111            query_frags = None
112            lib_frags = None
113
114        return (
115            query_in_lib,
116            query_in_lib_fract,
117            lib_in_query,
118            lib_in_query_fract,
119            query_frags,
120            lib_frags,
121        )
122
123    def fe_search(
124        self,
125        scan_list,
126        fe_lib,
127        precursor_mz_list=[],
128        use_mass_features=True,
129        peak_sep_da=None,
130        get_additional_metrics=True,
131        accumulate_results=False,
132        ms_params_key: str = "ms2",
133    ):
134        """
135        Search LC-MS MS2 spectra with FlashEntropy against a pre-built library.
136
137        Annotation thresholds and default peak-separation for cleaning come from
138        :class:`~corems.encapsulation.factory.processingSetting.SpectralSimilaritySearchSettings`
139        on this object's parameter tree
140        (``parameters.mass_spectrum[ms_params_key].spectral_similarity_search``).
141        Build ``fe_lib`` with the **same** settings so library index tolerances
142        match search cleaning (see
143        :meth:`~corems.molecular_id.search.database_interfaces.MSPInterface.get_metabolomics_spectra_library`
144        and
145        :meth:`~corems.molecular_id.search.database_interfaces.LCLipidLibraryInterface.get_lipid_library`).
146
147        Parameters
148        ----------
149        scan_list : list
150            List of scan numbers to search (must already be loaded in ``self._ms``).
151        fe_lib : :obj:`~ms_entropy.FlashEntropySearch`
152            FlashEntropy search instance. Prefer constructing it with
153            ``settings=settings_from_lcms(self)`` (or the same profile key as
154            *ms_params_key*) so FE build knobs match this search.
155        precursor_mz_list : list, optional
156            List of precursor m/z values to search, by default [], which implies
157            matched with mass features; to enable this use_mass_features must be True.
158        use_mass_features : bool, optional
159            If True, use mass features to get precursor m/z values, by default True.
160            If True, will add search results to mass features' ms2_similarity_results attribute.
161        peak_sep_da : float, optional
162            Minimum separation between m/z peaks spectra in Da. This needs match the
163            approximate resolution of the search spectra and the FlashEntropySearch
164            instance. If None (default), uses
165            ``parameters.mass_spectrum[ms_params_key].spectral_similarity_search.resolved_peak_sep_da``
166            (typically ``2 * max_ms2_tolerance_in_da`` from that settings object).
167        get_additional_metrics : bool, optional
168            If True, get additional metrics from FlashEntropy search, by default True.
169        accumulate_results : bool, optional
170            If True, accumulate results with existing spectral_search_results instead of
171            replacing them. This allows searching the same scans with multiple libraries
172            without overwriting previous results, by default False.
173        ms_params_key : str, optional
174            Key in ``parameters.mass_spectrum`` whose ``spectral_similarity_search``
175            settings are used (default ``"ms2"``). Use another key (e.g.
176            ``"ms2_cid"``) when multiple MS2 parameter profiles are configured
177            for different scan classes.
178
179        Returns
180        -------
181        None
182            Adds results to ``self.spectral_search_results`` and associates them
183            with mass features in ``self.mass_features`` when applicable.
184
185        Notes
186        -----
187        From the settings profile this method reads:
188
189        - ``ms2_min_fe_score`` — minimum entropy score to keep a hit
190        - ``include_fragment_types`` — lipid-style fragment-type metrics
191        - default peak separation when *peak_sep_da* is None
192
193        FlashEntropy **library** build fields (``max_ms2_tolerance_in_da``, etc.)
194        should already have been applied when constructing *fe_lib*. Configure
195        them on the same ``spectral_similarity_search`` instance before library
196        generation.
197
198        For an :class:`~corems.mass_spectra.factory.lc_class.LCMSCollection`,
199        sample parameters are equal at construction; build one shared library
200        with ``settings_from_lcms_collection(collection)``, then call
201        ``fe_search`` on each sample (or via collection pipeline operations).
202
203        Examples
204        --------
205        Single LCMSBase object (library + search share profile ``"ms2"``)::
206
207            from corems.encapsulation.factory.parameters import settings_from_lcms
208
209            settings = settings_from_lcms(lcms_obj, profile="ms2")
210            fe_lib, _meta = msp.get_metabolomics_spectra_library(
211                polarity="positive",
212                format="flashentropy",
213                settings=settings,
214            )
215            lcms_obj.fe_search(scan_list=ms2_scans, fe_lib=fe_lib)
216            # peak_sep / score gate from mass_spectrum["ms2"].spectral_similarity_search
217
218        Alternate MS2 profile (e.g. CID bag)::
219
220            lcms_obj.fe_search(
221                scan_list=cid_scans,
222                fe_lib=fe_lib_cid,
223                ms_params_key="ms2_cid",
224            )
225        """
226        # Annotation knobs from nested SpectralSimilaritySearchSettings (default profile "ms2")
227        if ms_params_key not in self.parameters.mass_spectrum:
228            raise KeyError(
229                f"ms_params_key={ms_params_key!r} not in parameters.mass_spectrum "
230                f"(keys={list(self.parameters.mass_spectrum)})"
231            )
232        ms2_p = self.parameters.mass_spectrum[ms_params_key].spectral_similarity_search
233        include_fragment_types = ms2_p.include_fragment_types
234        min_match_score = ms2_p.ms2_min_fe_score
235        if peak_sep_da is None:
236            peak_sep_da = ms2_p.resolved_peak_sep_da
237
238        # If precursor_mz_list is empty and use_mass_features is True, get precursor m/z values from mass features for each scan in scan_list
239        if use_mass_features and len(precursor_mz_list) == 0:
240            precursor_mz_list = []
241            for scan in scan_list:
242                mf_ids = [
243                    key
244                    for key, value in self.mass_features.items()
245                    if scan in value.ms2_mass_spectra
246                ]
247                precursor_mz = [
248                    value.mz
249                    for key, value in self.mass_features.items()
250                    if key in mf_ids
251                ]
252                precursor_mz_list.append(precursor_mz)
253
254        # Check that precursor_mz_list same length as scan_list, if not, raise error
255        if len(precursor_mz_list) != len(scan_list):
256            raise ValueError("Length of precursor_mz_list is not equal to scan_list.")
257
258        # Loop through each query spectrum / precursor match and save ids of db spectrum that are decent matches
259        overall_results_dict = {}
260        for i in np.arange(len(scan_list)):
261            scan_oi = scan_list[i]
262            if len(self._ms[scan_oi].mspeaks) > 0:
263                precursor_mzs = precursor_mz_list[i]
264                overall_results_dict[scan_oi] = {}
265                for precursor_mz in precursor_mzs:
266                    query_spectrum = fe_lib.clean_spectrum_for_search(
267                        precursor_mz=precursor_mz,
268                        peaks=np.vstack(
269                            (self._ms[scan_oi].mz_exp, self._ms[scan_oi].abundance)
270                        ).T,
271                        precursor_ions_removal_da=None,
272                        noise_threshold=self._ms[
273                            scan_oi
274                        ].parameters.mass_spectrum.noise_threshold_min_relative_abundance
275                        / 100,
276                        min_ms2_difference_in_da=peak_sep_da,
277                    )
278                    search_results = fe_lib.search(
279                        precursor_mz=precursor_mz,
280                        peaks=query_spectrum,
281                        ms1_tolerance_in_da=self.parameters.mass_spectrum[
282                            "ms1"
283                        ].molecular_search.max_ppm_error
284                        * 10**-6
285                        * precursor_mz,
286                        ms2_tolerance_in_da=peak_sep_da * 0.5,
287                        method={"identity"},
288                        precursor_ions_removal_da=None,
289                        noise_threshold=self._ms[
290                            scan_oi
291                        ].parameters.mass_spectrum.noise_threshold_min_relative_abundance
292                        / 100,
293                        target="cpu",
294                    )["identity_search"]
295                    match_inds = np.where(search_results > min_match_score)[0]
296
297                    # If any decent matches are found, add them to the results dictionary
298                    if len(match_inds) > 0:
299                        match_scores = search_results[match_inds]
300                        ref_ms_ids = [fe_lib[x]["id"] for x in match_inds]
301                        ref_mol_ids = [
302                            fe_lib[x]["molecular_data_id"] for x in match_inds
303                        ]
304                        ref_precursor_mzs = [
305                            fe_lib[x]["precursor_mz"] for x in match_inds
306                        ]
307                        ion_types = [fe_lib[x]["ion_type"] for x in match_inds]
308                        overall_results_dict[scan_oi][precursor_mz] = {
309                            "ref_mol_id": ref_mol_ids,
310                            "ref_ms_id": ref_ms_ids,
311                            "ref_precursor_mz": ref_precursor_mzs,
312                            "precursor_mz_error_ppm": [
313                                (precursor_mz - x) / precursor_mz * 10**6
314                                for x in ref_precursor_mzs
315                            ],
316                            "entropy_similarity": match_scores,
317                            "ref_ion_type": ion_types,
318                        }
319                        # Add database name, if present
320                        db_name = [fe_lib[x].get("database_name") for x in match_inds]
321                        if db_name is not None:
322                            overall_results_dict[scan_oi][precursor_mz].update(
323                                {"database_name": db_name}
324                            )
325                        if get_additional_metrics:
326                            more_match_quals = [
327                                self.get_more_match_quals(
328                                    self._ms[scan_oi].mz_exp,
329                                    fe_lib[x],
330                                    mz_tol_da=peak_sep_da,
331                                    include_fragment_types=include_fragment_types,
332                                )
333                                for x in match_inds
334                            ]
335                            overall_results_dict[scan_oi][precursor_mz].update(
336                                {
337                                    "query_mz_in_ref_n": [
338                                        x[0] for x in more_match_quals
339                                    ],
340                                    "query_mz_in_ref_fract": [
341                                        x[1] for x in more_match_quals
342                                    ],
343                                    "ref_mz_in_query_n": [
344                                        x[2] for x in more_match_quals
345                                    ],
346                                    "ref_mz_in_query_fract": [
347                                        x[3] for x in more_match_quals
348                                    ],
349                                }
350                            )
351                            if include_fragment_types:
352                                overall_results_dict[scan_oi][precursor_mz].update(
353                                    {
354                                        "query_frag_types": [
355                                            x[4] for x in more_match_quals
356                                        ],
357                                        "ref_frag_types": [
358                                            x[5] for x in more_match_quals
359                                        ],
360                                    }
361                                )
362
363        # Drop scans with no results from dictionary
364        overall_results_dict = {k: v for k, v in overall_results_dict.items() if v}
365
366        # Cast each entry as a MS2SearchResults object
367        for scan_id in overall_results_dict.keys():
368            for precursor_mz in overall_results_dict[scan_id].keys():
369                ms2_spectrum = self._ms[scan_id]
370                ms2_search_results = overall_results_dict[scan_id][precursor_mz]
371                overall_results_dict[scan_id][precursor_mz] = SpectrumSearchResults(
372                    ms2_spectrum, precursor_mz, ms2_search_results
373                )
374
375        # Add MS2SearchResults to the existing spectral search results dictionary
376        if accumulate_results:
377            # Merge results with existing spectral_search_results
378            for scan_id, precursor_dict in overall_results_dict.items():
379                if scan_id in self.spectral_search_results:
380                    # Scan already has results, merge precursor_mz dictionaries
381                    self.spectral_search_results[scan_id].update(precursor_dict)
382                else:
383                    # New scan, add entire dictionary
384                    self.spectral_search_results[scan_id] = precursor_dict
385        else:
386            # Replace existing results (original behavior)
387            self.spectral_search_results.update(overall_results_dict)
388
389        # If there are mass features, associate the results with each mass feature
390        if len(self.mass_features) > 0:
391            # Determine which results to associate with mass features
392            if accumulate_results:
393                # When accumulating, only associate new results from this search
394                # to avoid duplicating previously associated results
395                results_to_associate = overall_results_dict
396            else:
397                # When not accumulating, clear existing associations and re-associate all results
398                for mass_feature_id in self.mass_features.keys():
399                    self.mass_features[mass_feature_id].ms2_similarity_results = []
400                results_to_associate = self.spectral_search_results
401
402            for mass_feature_id, mass_feature in self.mass_features.items():
403                scan_ids = mass_feature.ms2_scan_numbers
404                for ms2_scan_id in scan_ids:
405                    precursor_mz = mass_feature.mz
406                    try:
407                        results_to_associate[ms2_scan_id][precursor_mz]
408                    except KeyError:
409                        pass
410                    else:
411                        self.mass_features[
412                            mass_feature_id
413                        ].ms2_similarity_results.append(
414                            results_to_associate[ms2_scan_id][precursor_mz]
415                        )

Methods for searching LCMS spectra.

This class is designed to be a mixin class for the ~corems.mass_spectra.factory.lc_class.LCMSBase class.

@staticmethod
def get_more_match_quals(query_mz_arr, lib_entry, mz_tol_da=0.1, include_fragment_types=False):
 17    @staticmethod
 18    def get_more_match_quals(
 19        query_mz_arr, lib_entry, mz_tol_da=0.1, include_fragment_types=False
 20    ):
 21        """
 22        Return additional match qualities between query and library entry.
 23
 24        Parameters
 25        ----------
 26        query_mz_arr : np.array
 27            Array of query spectrum. Shape (N, 2), with m/z in the first column
 28            and abundance in the second.
 29        lib_entry : dict
 30            Library spectrum entry, with 'mz' key containing the spectrum in
 31            the format (mz, abundance),(mz, abundance), i.e. from MetabRef.
 32        mz_tol_da : float, optional
 33            Tolerance in Da for matching peaks (in MS2). Default is 0.1.
 34        include_fragment_types : bool, optional
 35            If True, include fragment type comparisons in output.
 36            Defaults to False.
 37
 38        Returns
 39        -------
 40        tuple
 41            Tuple of (query_in_lib, query_in_lib_fract, lib_in_query, lib_in_query_fract, query_frags, lib_frags, lib_precursor_mz).
 42
 43        Notes
 44        -----
 45        query_in_lib : int
 46            Number of peaks in query that are present in the library entry (within mz_tol_da).
 47        query_in_lib_fract : float
 48            Fraction of peaks in query that are present in the library entry (within mz_tol_da).
 49        lib_in_query : int
 50            Number of peaks in the library entry that are present in the query (within mz_tol_da).
 51        lib_in_query_fract : float
 52            Fraction of peaks in the library entry that are present in the query (within mz_tol_da).
 53        query_frags : list
 54            List of unique fragment types present in the query, generally 'MLF' or 'LSF' or both.
 55        lib_frags : list
 56            List of unique fragment types present in the library entry, generally 'MLF' or 'LSF' or both.
 57
 58        Raises
 59        ------
 60        ValueError
 61            If library entry does not have 'fragment_types' key and include_fragment_types is True.
 62
 63        """
 64
 65        if "mz" in lib_entry.keys():
 66            # Get the original mz values from the library entry
 67            lib_mzs = np.array(
 68                re.findall(r"\(([^,]+),([^)]+)\)", lib_entry["mz"]), dtype=float
 69            ).reshape(-1, 2)[:, 0]
 70        elif "peaks" in lib_entry.keys() and lib_entry["peaks"] is not None:
 71            lib_mzs = lib_entry["peaks"][:, 0]
 72
 73        # Get count and fraction of peaks in query that are in lib entry
 74        query_in_lib = 0
 75        for peak in query_mz_arr:
 76            if np.any(np.isclose(lib_mzs, peak, atol=mz_tol_da)):
 77                query_in_lib += 1
 78        query_in_lib_fract = query_in_lib / len(query_mz_arr)
 79
 80        # Get count and fraction of peaks in lib that are in query
 81        lib_in_query = 0
 82        for peak in lib_mzs:
 83            if np.any(np.isclose(query_mz_arr, peak, atol=mz_tol_da)):
 84                lib_in_query += 1
 85        lib_in_query_fract = lib_in_query / len(lib_mzs)
 86
 87        if include_fragment_types:
 88            # Check that fragment types are present in the library entry
 89            if "fragment_types" not in lib_entry.keys():
 90                raise ValueError(
 91                    "Flash entropy library entry must have 'fragment_types' key to include fragment types in output."
 92                )
 93
 94            # Get types of fragments in the lib entry and convert it to a list on ", "
 95            lib_frags = [x.strip() for x in lib_entry["fragment_types"].split(",")]
 96            # make list of the fragment types that are present in the query spectrum
 97            lib_in_query_ids = list(
 98                set(
 99                    [
100                        ind
101                        for ind, x in enumerate(lib_mzs)
102                        if len(np.where(np.isclose(query_mz_arr, x, atol=mz_tol_da))[0])
103                        > 0
104                    ]
105                )
106            )
107            query_frags = list(set([lib_frags[x] for x in lib_in_query_ids]))
108            lib_frags = list(set(lib_frags))
109
110        else:
111            query_frags = None
112            lib_frags = None
113
114        return (
115            query_in_lib,
116            query_in_lib_fract,
117            lib_in_query,
118            lib_in_query_fract,
119            query_frags,
120            lib_frags,
121        )

Return additional match qualities between query and library entry.

Parameters
  • query_mz_arr (np.array): Array of query spectrum. Shape (N, 2), with m/z in the first column and abundance in the second.
  • lib_entry (dict): Library spectrum entry, with 'mz' key containing the spectrum in the format (mz, abundance),(mz, abundance), i.e. from MetabRef.
  • mz_tol_da (float, optional): Tolerance in Da for matching peaks (in MS2). Default is 0.1.
  • include_fragment_types (bool, optional): If True, include fragment type comparisons in output. Defaults to False.
Returns
  • tuple: Tuple of (query_in_lib, query_in_lib_fract, lib_in_query, lib_in_query_fract, query_frags, lib_frags, lib_precursor_mz).
Notes

query_in_lib : int Number of peaks in query that are present in the library entry (within mz_tol_da). query_in_lib_fract : float Fraction of peaks in query that are present in the library entry (within mz_tol_da). lib_in_query : int Number of peaks in the library entry that are present in the query (within mz_tol_da). lib_in_query_fract : float Fraction of peaks in the library entry that are present in the query (within mz_tol_da). query_frags : list List of unique fragment types present in the query, generally 'MLF' or 'LSF' or both. lib_frags : list List of unique fragment types present in the library entry, generally 'MLF' or 'LSF' or both.

Raises
  • ValueError: If library entry does not have 'fragment_types' key and include_fragment_types is True.