corems.molecular_id.search.molecularFormulaSearch

   1__author__ = "Yuri E. Corilo"
   2__date__ = "Jul 29, 2019"
   3
   4
   5from typing import Callable, List, NamedTuple, Optional, Tuple, Union
   6
   7import tqdm
   8
   9from corems import chunks, timeit
  10from corems.encapsulation.constant import Labels
  11from corems.encapsulation.factory.processingSetting import validate_used_atoms_keys
  12from corems.molecular_formula.factory.MolecularFormulaFactory import MolecularFormula
  13from corems.molecular_id.factory.MolecularLookupTable import MolecularCombinations
  14from corems.molecular_id.factory.molecularSQL import MolForm_SQL
  15from corems.ms_peak.factory.MSPeakClasses import _MSPeak
  16
  17last_error = 0
  18last_dif = 0
  19closest_error = 0
  20error_average = 0
  21nbValues = 0
  22
  23
  24class FormulaSearchJob(NamedTuple):
  25    """One ion-type candidate batch for molecular formula search.
  26
  27    Parameters
  28    ----------
  29    ion_charge : int
  30        Signed search charge for this batch.
  31    ion_type : str
  32        ``Labels`` ion type (protonated/deprotonated, radical, or adduct).
  33    adduct_atom : str or None
  34        Adduct element symbol when ``ion_type`` is adduct; otherwise None.
  35    candidate_formulas : dict
  36        Mapping of nominal m/z to formula objects, as consumed by ``run_search``.
  37    progress : str
  38        tqdm description for this class / ion type / charge.
  39    """
  40
  41    ion_charge: int
  42    ion_type: str
  43    adduct_atom: Optional[str]
  44    candidate_formulas: dict
  45    progress: str
  46
  47
  48class SearchMolecularFormulas:
  49    """Class for searching molecular formulas in a mass spectrum.
  50
  51    Parameters
  52    ----------
  53    mass_spectrum_obj : MassSpectrum
  54        The mass spectrum object.
  55    sql_db : MolForm_SQL, optional
  56        The SQL database object, by default None.
  57    first_hit : bool, optional
  58        Flag to indicate whether to skip peaks that already have a molecular formula assigned, by default False.
  59    find_isotopologues : bool, optional
  60        Flag to indicate whether to find isotopologues, by default True.
  61
  62    Attributes
  63    ----------
  64    mass_spectrum_obj : MassSpectrum
  65        The mass spectrum object.
  66    sql_db : MolForm_SQL
  67        The SQL database object.
  68    first_hit : bool
  69        Flag to indicate whether to skip peaks that already have a molecular formula assigned.
  70    find_isotopologues : bool
  71        Flag to indicate whether to find isotopologues.
  72
  73
  74    Methods
  75    -------
  76    * run_search().
  77        Run the molecular formula search.
  78    * run_worker_mass_spectrum().
  79        Run the molecular formula search on the mass spectrum object.
  80    * run_worker_ms_peaks().
  81        Run the molecular formula search on the given list of mass spectrum peaks.
  82    * database_to_dict().
  83        Convert the database results to a dictionary.
  84    * run_molecular_formula().
  85        Run the molecular formula search on the given list of mass spectrum peaks.
  86    * search_mol_formulas().
  87        Search for molecular formulas in the mass spectrum.
  88    * ion_charges_for_search().
  89        Build the signed charge tuple from polarity and absolute min/max settings.
  90
  91    """
  92
  93    def __init__(
  94        self,
  95        mass_spectrum_obj,
  96        sql_db=None,
  97        first_hit: bool = False,
  98        find_isotopologues: bool = True,
  99    ):
 100        self.first_hit = first_hit
 101
 102        self.find_isotopologues = find_isotopologues
 103
 104        self.mass_spectrum_obj = mass_spectrum_obj
 105
 106        # Catch in-place usedAtoms mutations after settings construction
 107        validate_used_atoms_keys(mass_spectrum_obj.molecular_search_settings.usedAtoms)
 108
 109        if not sql_db:
 110            self.sql_db = MolForm_SQL(
 111                url=mass_spectrum_obj.molecular_search_settings.url_database
 112            )
 113
 114        else:
 115            self.sql_db = sql_db
 116
 117    def __enter__(self):
 118        """Open the SQL database connection."""
 119        return self
 120
 121    def __exit__(self, exc_type, exc_val, exc_tb):
 122        """Close the SQL database connection."""
 123        self.sql_db.close()
 124
 125        return False
 126
 127    @staticmethod
 128    def ion_charges_for_search(
 129        polarity: Union[int, str],
 130        min_ion_charge: int = 1,
 131        max_ion_charge: int = 1,
 132    ) -> Tuple[int, ...]:
 133        """Signed ion charges to use for molecular formula search.
 134
 135        Absolute charges from ``min_ion_charge`` through ``max_ion_charge``
 136        (inclusive) are multiplied by the polarity sign. No charge-state
 137        determination is performed; every selected peak is searched at each
 138        returned charge and candidates are ranked by formula scores.
 139
 140        Parameters
 141        ----------
 142        polarity : int or str
 143            Mass-spectrum polarity as ``+1`` / ``-1``, or LCMS polarity string
 144            ``"positive"`` / ``"negative"``.
 145        min_ion_charge : int, optional
 146            Minimum absolute charge (default 1).
 147        max_ion_charge : int, optional
 148            Maximum absolute charge (default 1). Must be >= ``min_ion_charge``.
 149
 150        Returns
 151        -------
 152        tuple of int
 153            Signed charges sorted by ``(abs(z), z)`` so lower |z| is first
 154            (relevant when ``first_hit=True``).
 155
 156        Raises
 157        ------
 158        ValueError
 159            If polarity is invalid, min < 1, or max < min.
 160        """
 161        min_z = int(min_ion_charge)
 162        max_z = int(max_ion_charge)
 163        if min_z < 1:
 164            raise ValueError("min_ion_charge must be >= 1")
 165        if max_z < min_z:
 166            raise ValueError("max_ion_charge must be >= min_ion_charge")
 167
 168        if isinstance(polarity, str):
 169            p = polarity.strip().lower()
 170            if p == "positive":
 171                sign = 1
 172            elif p == "negative":
 173                sign = -1
 174            else:
 175                raise ValueError(
 176                    "polarity string must be 'positive' or 'negative', "
 177                    f"got {polarity!r}"
 178                )
 179        else:
 180            pol = int(polarity)
 181            if pol == 0:
 182                raise ValueError("polarity must be non-zero")
 183            sign = 1 if pol > 0 else -1
 184
 185        charges = [sign * z for z in range(min_z, max_z + 1)]
 186        return tuple(sorted(charges, key=lambda z: (abs(z), z)))
 187
 188    @staticmethod
 189    def _iter_ion_type_jobs(settings, dict_res, classe_str, ion_charge):
 190        """Yield ion-type candidate jobs for one heteroatom class at one charge.
 191
 192        Encodes shared search policy: protonated and radical follow the
 193        corresponding settings flags; adduct search runs only when adducts
 194        are enabled and ``abs(ion_charge) == 1``. Empty candidate lists and
 195        missing ``dict_res`` keys are skipped.
 196
 197        Parameters
 198        ----------
 199        settings : MolecularFormulaSearchSettings
 200            Molecular formula search settings (ion-type flags).
 201        dict_res : dict
 202            Result of ``database_to_dict`` for this charge: protonated/radical
 203            map class string to a nominal-m/z query dict; adduct maps adduct
 204            atom to that class mapping.
 205        classe_str : str
 206            Heteroatom class key as stored in ``dict_res``.
 207        ion_charge : int
 208            Signed ion charge for this batch.
 209
 210        Yields
 211        ------
 212        FormulaSearchJob
 213            One job per enabled ion type (and per adduct atom) with candidates.
 214        """
 215        if settings.isProtonated:
 216            by_class = (dict_res.get(Labels.protonated_de_ion) or {}).get(classe_str)
 217            if by_class:
 218                yield FormulaSearchJob(
 219                    ion_charge=ion_charge,
 220                    ion_type=Labels.protonated_de_ion,
 221                    adduct_atom=None,
 222                    candidate_formulas=by_class,
 223                    progress=(
 224                        "Started molecular formula search for class %s, "
 225                        "(de)protonated z=%s " % (classe_str, ion_charge)
 226                    ),
 227                )
 228
 229        if settings.isRadical:
 230            by_class = (dict_res.get(Labels.radical_ion) or {}).get(classe_str)
 231            if by_class:
 232                yield FormulaSearchJob(
 233                    ion_charge=ion_charge,
 234                    ion_type=Labels.radical_ion,
 235                    adduct_atom=None,
 236                    candidate_formulas=by_class,
 237                    progress=(
 238                        "Started molecular formula search for class %s, "
 239                        "radical z=%s " % (classe_str, ion_charge)
 240                    ),
 241                )
 242
 243        # Adduct search only at |z|==1; multi-charge adducts are not modeled.
 244        if settings.isAdduct and abs(ion_charge) == 1:
 245            dict_atoms_formulas = dict_res.get(Labels.adduct_ion) or {}
 246            for adduct_atom, dict_by_class in dict_atoms_formulas.items():
 247                by_class = dict_by_class.get(classe_str)
 248                if by_class:
 249                    yield FormulaSearchJob(
 250                        ion_charge=ion_charge,
 251                        ion_type=Labels.adduct_ion,
 252                        adduct_atom=adduct_atom,
 253                        candidate_formulas=by_class,
 254                        progress=(
 255                            "Started molecular formula search for class %s, "
 256                            "adduct z=%s " % (classe_str, ion_charge)
 257                        ),
 258                    )
 259
 260    @staticmethod
 261    def _run_formula_search_jobs(
 262        classes,
 263        nominal_mzs,
 264        mf_search_settings,
 265        search_charges,
 266        sql_db,
 267        apply_fn: Callable[[FormulaSearchJob], None],
 268    ):
 269        """Run charge × class × ion-type search, applying each candidate batch.
 270
 271        Loads the formula database once per charge per class chunk. Does not
 272        close ``sql_db``. tqdm iterates class chunks; the description is the
 273        job ``progress`` string.
 274
 275        Parameters
 276        ----------
 277        classes : list
 278            Heteroatom class tuples from ``MolecularCombinations.runworker``.
 279        nominal_mzs : list
 280            Nominal m/z values to load from the database.
 281        mf_search_settings : MolecularFormulaSearchSettings
 282            Search settings (chunk size, ion-type flags, verbosity).
 283        search_charges : sequence of int
 284            Signed charges from ``ion_charges_for_search``.
 285        sql_db : MolForm_SQL
 286            Open formula database connection.
 287        apply_fn : callable
 288            Called with each ``FormulaSearchJob``. Path-specific: DI uses
 289            ``run_search``; LC uses ``search_spectra_against_candidates``.
 290        """
 291        verbose = mf_search_settings.verbose_processing
 292        for classe_chunk in chunks(classes, mf_search_settings.db_chunk_size):
 293            classes_str_list = [class_tuple[0] for class_tuple in classe_chunk]
 294            for ion_charge in search_charges:
 295                dict_res = SearchMolecularFormulas.database_to_dict(
 296                    classes_str_list,
 297                    nominal_mzs,
 298                    mf_search_settings,
 299                    ion_charge,
 300                    sql_db=sql_db,
 301                )
 302                pbar = tqdm.tqdm(classe_chunk, disable=not verbose)
 303                for classe_tuple in pbar:
 304                    classe_str = classe_tuple[0]
 305                    for job in SearchMolecularFormulas._iter_ion_type_jobs(
 306                        mf_search_settings, dict_res, classe_str, ion_charge
 307                    ):
 308                        pbar.set_description_str(desc=job.progress, refresh=True)
 309                        apply_fn(job)
 310
 311    def run_search(
 312        self,
 313        mspeaks: list,
 314        query: dict,
 315        min_abundance: float,
 316        ion_type: str,
 317        ion_charge: int,
 318        adduct_atom=None,
 319    ):
 320        """Run the molecular formula search.
 321
 322        Parameters
 323        ----------
 324        mspeaks : list of MSPeak
 325            The list of mass spectrum peaks.
 326        query : dict
 327            The query dictionary containing the possible molecular formulas.
 328        min_abundance : float
 329            The minimum abundance threshold.
 330        ion_type : str
 331            The ion type.
 332        ion_charge : int
 333            The ion charge.
 334        adduct_atom : str, optional
 335            The adduct atom, by default None.
 336        """
 337
 338        def get_formulas(nominal_overlay: float = 0.1):
 339            """
 340            Get the list of formulas based on the nominal overlay.
 341
 342            Parameters
 343            ----------
 344            nominal_overlay : float, optional
 345                The nominal overlay, by default 0.1.
 346
 347            Returns
 348            -------
 349            list
 350                The list of formulas.
 351            """
 352            nominal_mz = ms_peak.nominal_mz_exp
 353
 354            defect_mass = ms_peak.mz_exp - nominal_mz
 355            nominal_masses = [nominal_mz]
 356
 357            if (defect_mass) >= 1 - nominal_overlay:
 358                nominal_masses.append(nominal_mz + 1)
 359            elif (defect_mass) <= nominal_overlay:
 360                nominal_masses.append(nominal_mz - 1)
 361
 362            list_formulas_candidates = []
 363
 364            for nominal_mass in nominal_masses:
 365                if nominal_mass in query.keys():
 366                    list_formulas_candidates.extend(query.get(nominal_mass))
 367
 368            return list_formulas_candidates
 369
 370        all_assigned_indexes = list()
 371
 372        # molecular_search_settings = self.mass_spectrum_obj.molecular_search_settings
 373
 374        search_molfrom = SearchMolecularFormulaWorker(
 375            find_isotopologues=self.find_isotopologues
 376        )
 377
 378        for ms_peak in mspeaks:
 379            # already assigned a molecular formula
 380            if self.first_hit:
 381                if ms_peak.is_assigned:
 382                    continue
 383
 384            ms_peak_indexes = search_molfrom.find_formulas(
 385                get_formulas(),
 386                min_abundance,
 387                self.mass_spectrum_obj,
 388                ms_peak,
 389                ion_type,
 390                ion_charge,
 391                adduct_atom,
 392            )
 393
 394            all_assigned_indexes.extend(ms_peak_indexes)
 395
 396        # all_assigned_indexes = MolecularFormulaSearchFilters().filter_isotopologue(all_assigned_indexes, self.mass_spectrum_obj)
 397
 398        # all_assigned_indexes = MolecularFormulaSearchFilters().filter_kendrick(all_assigned_indexes, self.mass_spectrum_obj)
 399
 400        # MolecularFormulaSearchFilters().check_min_peaks(all_assigned_indexes, self.mass_spectrum_obj)
 401        # filter per min peaks per mono isotopic class
 402
 403    def run_worker_mass_spectrum(self):
 404        """Run the molecular formula search on the mass spectrum object."""
 405        self.run_molecular_formula(
 406            self.mass_spectrum_obj.sort_by_abundance(),
 407            print_time=self.mass_spectrum_obj.molecular_search_settings.verbose_processing,
 408        )
 409
 410    def run_worker_ms_peaks(self, ms_peaks):
 411        """Run the molecular formula search on the given list of mass spectrum peaks.
 412
 413        Parameters
 414        ----------
 415        ms_peaks : list of MSPeak
 416            The list of mass spectrum peaks.
 417        """
 418        self.run_molecular_formula(
 419            ms_peaks,
 420            print_time=self.mass_spectrum_obj.molecular_search_settings.verbose_processing,
 421        )
 422
 423    @staticmethod
 424    def database_to_dict(
 425        classe_str_list,
 426        nominal_mzs,
 427        mf_search_settings,
 428        ion_charge,
 429        sql_db=None,
 430    ):
 431        """Convert the database results to a dictionary.
 432
 433        Parameters
 434        ----------
 435        classe_str_list : list
 436            The list of class strings.
 437        nominal_mzs : list
 438            The list of nominal m/z values.
 439        mf_search_settings : MolecularFormulaSearchSettings
 440            The molecular formula search settings.
 441        ion_charge : int
 442            The ion charge.
 443        sql_db : MolForm_SQL, optional
 444            The SQL database object, by default None. If None, a new MolForm_SQL object will be created.
 445
 446        Returns
 447        -------
 448        dict
 449            The dictionary containing the database results.
 450        """
 451        owns_db = sql_db is None
 452        if owns_db:
 453            sql_db = MolForm_SQL(url=mf_search_settings.url_database)
 454
 455        try:
 456            dict_res = {}
 457
 458            if mf_search_settings.isProtonated:
 459                dict_res[Labels.protonated_de_ion] = sql_db.get_dict_by_classes(
 460                    classe_str_list,
 461                    Labels.protonated_de_ion,
 462                    nominal_mzs,
 463                    ion_charge,
 464                    mf_search_settings,
 465                )
 466
 467            if mf_search_settings.isRadical:
 468                dict_res[Labels.radical_ion] = sql_db.get_dict_by_classes(
 469                    classe_str_list,
 470                    Labels.radical_ion,
 471                    nominal_mzs,
 472                    ion_charge,
 473                    mf_search_settings,
 474                )
 475
 476            if mf_search_settings.isAdduct:
 477                adduct_list = (
 478                    mf_search_settings.adduct_atoms_neg
 479                    if ion_charge < 0
 480                    else mf_search_settings.adduct_atoms_pos
 481                )
 482                dict_res[Labels.adduct_ion] = sql_db.get_dict_by_classes(
 483                    classe_str_list,
 484                    Labels.adduct_ion,
 485                    nominal_mzs,
 486                    ion_charge,
 487                    mf_search_settings,
 488                    adducts=adduct_list,
 489                )
 490
 491            return dict_res
 492        finally:
 493            if owns_db:
 494                sql_db.close()
 495
 496    @timeit(print_time=True)
 497    def run_molecular_formula(self, ms_peaks, **kwargs):
 498        """Run the molecular formula search on the given list of mass spectrum peaks.
 499
 500        Parameters
 501        ----------
 502        ms_peaks : list of MSPeak
 503            The list of mass spectrum peaks.
 504        **kwargs
 505            Additional keyword arguments.
 506            Most notably, print_time, which is a boolean flag to indicate whether to print the time
 507            and passed to the timeit decorator.
 508        """
 509        settings = self.mass_spectrum_obj.molecular_search_settings
 510        # Re-check after possible post-construction mutation of settings fields
 511        settings.validate_ion_charge_settings()
 512        # Polarity must come from the mass spectrum (DI often -1 from data).
 513        # Do not default to +1; LC already requires polarity on the LCMS object.
 514        polarity_raw = self.mass_spectrum_obj.polarity
 515        if polarity_raw is None:
 516            raise ValueError(
 517                "Mass spectrum polarity must be set from the data before "
 518                "formula search (e.g. +1 or -1). There is no default polarity."
 519            )
 520        try:
 521            polarity = int(polarity_raw)
 522        except (TypeError, ValueError) as exc:
 523            raise ValueError(
 524                "Mass spectrum polarity must be an integer +1 or -1 "
 525                f"(got {polarity_raw!r})."
 526            ) from exc
 527        if polarity == 0:
 528            raise ValueError("Mass spectrum polarity must be non-zero (+1 or -1).")
 529        # Multi-charge formula search: expand absolute min..max (default 1..1).
 530        # Does not rewrite MSPeak.polarity (remains ±1).
 531        search_charges = self.ion_charges_for_search(
 532            polarity,
 533            min_ion_charge=settings.min_ion_charge,
 534            max_ion_charge=settings.max_ion_charge,
 535        )
 536        min_abundance = self.mass_spectrum_obj.min_abundance
 537        nominal_mzs = self.mass_spectrum_obj.nominal_mz
 538
 539        # reset average error, only relevant is average mass error method is being used
 540        SearchMolecularFormulaWorker(
 541            find_isotopologues=self.find_isotopologues
 542        ).reset_error(self.mass_spectrum_obj)
 543
 544        # check database for all possible molecular formula combinations based on the setting passed to self.mass_spectrum_obj.molecular_search_settings
 545        classes = MolecularCombinations(self.sql_db).runworker(
 546            settings, print_time=settings.verbose_processing
 547        )
 548
 549        def apply_job(job):
 550            self.run_search(
 551                ms_peaks,
 552                job.candidate_formulas,
 553                min_abundance,
 554                job.ion_type,
 555                job.ion_charge,
 556                adduct_atom=job.adduct_atom,
 557            )
 558
 559        self._run_formula_search_jobs(
 560            classes,
 561            nominal_mzs,
 562            settings,
 563            search_charges,
 564            self.sql_db,
 565            apply_job,
 566        )
 567        self.sql_db.close()
 568
 569    def search_mol_formulas(
 570        self,
 571        possible_formulas_list: List[MolecularFormula],
 572        ion_type: str,
 573        neutral_molform=True,
 574        find_isotopologues=True,
 575        adduct_atom=None,
 576    ) -> List[_MSPeak]:
 577        """Search for molecular formulas in the mass spectrum.
 578
 579        Parameters
 580        ----------
 581        possible_formulas_list : list of MolecularFormula
 582            The list of possible molecular formulas.
 583        ion_type : str
 584            The ion type.
 585        neutral_molform : bool, optional
 586            Flag to indicate whether the molecular formulas are neutral, by default True.
 587        find_isotopologues : bool, optional
 588            Flag to indicate whether to find isotopologues, by default True.
 589        adduct_atom : str, optional
 590            The adduct atom, by default None.
 591
 592        Returns
 593        -------
 594        list of MSPeak
 595            The list of mass spectrum peaks with assigned molecular formulas.
 596        """
 597        # neutral_molform: some reference files already present the formula on ion mode, for instance, bruker reference files
 598        #    if that is the case than turn neutral_molform off
 599
 600        SearchMolecularFormulaWorker(find_isotopologues=find_isotopologues).reset_error(
 601            self.mass_spectrum_obj
 602        )
 603
 604        initial_min_peak_bool = (
 605            self.mass_spectrum_obj.molecular_search_settings.use_min_peaks_filter
 606        )
 607        initial_runtime_kendrick_filter = (
 608            self.mass_spectrum_obj.molecular_search_settings.use_runtime_kendrick_filter
 609        )
 610
 611        # Are the following 3 lines redundant?
 612        self.mass_spectrum_obj.molecular_search_settings.use_min_peaks_filter = False
 613        self.mass_spectrum_obj.molecular_search_settings.use_min_peaks_filter = (
 614            False  # TODO check this line
 615        )
 616        self.mass_spectrum_obj.molecular_search_settings.use_runtime_kendrick_filter = (
 617            False
 618        )
 619
 620        possible_formulas_dict_nm = {}
 621
 622        for mf in possible_formulas_list:
 623            if neutral_molform:
 624                # Bin by the nominal m/z of the requested ion hypothesis so
 625                # candidates land in the same bin run_search looks them up under.
 626                if ion_type == Labels.radical_ion:
 627                    nm = int(mf.radical_mz)
 628                elif ion_type == Labels.adduct_ion and adduct_atom:
 629                    nm = int(mf.adduct_mz(adduct_atom))
 630                else:
 631                    nm = int(mf.protonated_mz)
 632            else:
 633                nm = int(mf.mz_nominal_calc)
 634
 635            if nm in possible_formulas_dict_nm.keys():
 636                possible_formulas_dict_nm[nm].append(mf)
 637
 638            else:
 639                possible_formulas_dict_nm[nm] = [mf]
 640
 641        min_abundance = self.mass_spectrum_obj.min_abundance
 642
 643        ion_type = ion_type
 644
 645        self.run_search(
 646            self.mass_spectrum_obj,
 647            possible_formulas_dict_nm,
 648            min_abundance,
 649            ion_type,
 650            self.mass_spectrum_obj.polarity,
 651            adduct_atom=adduct_atom,
 652        )
 653
 654        self.mass_spectrum_obj.molecular_search_settings.use_min_peaks_filter = (
 655            initial_min_peak_bool
 656        )
 657        self.mass_spectrum_obj.molecular_search_settings.use_runtime_kendrick_filter = (
 658            initial_runtime_kendrick_filter
 659        )
 660
 661        mspeaks = [mspeak for mspeak in self.mass_spectrum_obj if mspeak.is_assigned]
 662
 663        self.sql_db.close()
 664
 665        return mspeaks
 666
 667
 668class SearchMolecularFormulaWorker:
 669    """Class for searching molecular formulas in a mass spectrum.
 670
 671    Parameters
 672    ----------
 673    find_isotopologues : bool, optional
 674        Flag to indicate whether to find isotopologues, by default True.
 675
 676    Attributes
 677    ----------
 678    find_isotopologues : bool
 679        Flag to indicate whether to find isotopologues.
 680
 681    Methods
 682    -------
 683    * reset_error().
 684        Reset the error variables.
 685    * set_last_error().
 686        Set the last error.
 687    * find_formulas().
 688        Find the formulas.
 689    * calc_error().
 690        Calculate the error.
 691    """
 692
 693    # TODO add reset error function
 694    # needs this wraper to pass the class to multiprocessing
 695
 696    def __init__(self, find_isotopologues=True):
 697        self.find_isotopologues = find_isotopologues
 698
 699    def __call__(self, args):
 700        """Call the find formulas function.
 701
 702        Parameters
 703        ----------
 704        args : tuple
 705            The arguments.
 706
 707        Returns
 708        -------
 709        list
 710            The list of mass spectrum peaks with assigned molecular formulas.
 711        """
 712        return self.find_formulas(*args)  # ,args[1]
 713
 714    def reset_error(self, mass_spectrum_obj):
 715        """Reset the error variables.
 716
 717        Parameters
 718        ----------
 719        mass_spectrum_obj : MassSpectrum
 720            The mass spectrum object.
 721
 722        Notes
 723        -----
 724        This function resets the error variables for the given mass spectrum object.
 725        """
 726        global last_error, last_dif, closest_error, error_average, nbValues
 727        last_error, last_dif, closest_error, nbValues = 0.0, 0.0, 0.0, 0.0
 728
 729    def set_last_error(self, error, mass_spectrum_obj):
 730        """Set the last error.
 731
 732        Parameters
 733        ----------
 734        error : float
 735            The error.
 736        mass_spectrum_obj : MassSpectrum
 737            The mass spectrum object.
 738        """
 739        # set the changes to the global variables, not internal ones
 740        global last_error, last_dif, closest_error, error_average, nbValues
 741
 742        if mass_spectrum_obj.molecular_search_settings.error_method == "distance":
 743            dif = error - last_error
 744            if dif < last_dif:
 745                last_dif = dif
 746                closest_error = error
 747                mass_spectrum_obj.molecular_search_settings.min_ppm_error = (
 748                    closest_error
 749                    - mass_spectrum_obj.molecular_search_settings.mz_error_range
 750                )
 751                mass_spectrum_obj.molecular_search_settings.max_ppm_error = (
 752                    closest_error
 753                    + mass_spectrum_obj.molecular_search_settings.mz_error_range
 754                )
 755
 756        elif mass_spectrum_obj.molecular_search_settings.error_method == "lowest":
 757            if error < last_error:
 758                mass_spectrum_obj.molecular_search_settings.min_ppm_error = (
 759                    error - mass_spectrum_obj.molecular_search_settings.mz_error_range
 760                )
 761                mass_spectrum_obj.molecular_search_settings.max_ppm_error = (
 762                    error + mass_spectrum_obj.molecular_search_settings.mz_error_range
 763                )
 764                last_error = error
 765
 766        elif mass_spectrum_obj.molecular_search_settings.error_method == "symmetrical":
 767            mass_spectrum_obj.molecular_search_settings.min_ppm_error = (
 768                mass_spectrum_obj.molecular_search_settings.mz_error_average
 769                - mass_spectrum_obj.molecular_search_settings.mz_error_range
 770            )
 771            mass_spectrum_obj.molecular_search_settings.max_ppm_error = (
 772                mass_spectrum_obj.molecular_search_settings.mz_error_average
 773                + mass_spectrum_obj.molecular_search_settings.mz_error_range
 774            )
 775
 776        elif mass_spectrum_obj.molecular_search_settings.error_method == "average":
 777            nbValues += 1
 778            error_average = error_average + ((error - error_average) / nbValues)
 779            mass_spectrum_obj.molecular_search_settings.min_ppm_error = (
 780                error_average
 781                - mass_spectrum_obj.molecular_search_settings.mz_error_range
 782            )
 783            mass_spectrum_obj.molecular_search_settings.max_ppm_error = (
 784                error_average
 785                + mass_spectrum_obj.molecular_search_settings.mz_error_range
 786            )
 787
 788        else:
 789            # using set mass_spectrum_obj.molecular_search_settings.min_ppm_error  and max_ppm_error range
 790            pass
 791
 792        # returns the error based on the selected method at mass_spectrum_obj.molecular_search_settings.method
 793
 794    @staticmethod
 795    def calc_error(mz_exp, mz_calc, method="ppm"):
 796        """Calculate the error.
 797
 798        Parameters
 799        ----------
 800        mz_exp : float
 801            The experimental m/z value.
 802        mz_calc : float
 803            The calculated m/z value.
 804        method : str, optional
 805            The method, by default 'ppm'.
 806
 807        Raises
 808        -------
 809        Exception
 810            If the method is not ppm or ppb.
 811
 812        Returns
 813        -------
 814        float
 815            The error.
 816        """
 817
 818        if method == "ppm":
 819            multi_factor = 1_000_000
 820
 821        elif method == "ppb":
 822            multi_factor = 1_000_000_000
 823
 824        elif method == "perc":
 825            multi_factor = 100
 826
 827        else:
 828            raise Exception(
 829                "method needs to be ppm or ppb, you have entered %s" % method
 830            )
 831
 832        if mz_exp:
 833            return ((mz_exp - mz_calc) / mz_calc) * multi_factor
 834
 835        else:
 836            raise Exception("Please set mz_calc first")
 837
 838    def find_formulas(
 839        self,
 840        formulas,
 841        min_abundance,
 842        mass_spectrum_obj,
 843        ms_peak,
 844        ion_type,
 845        ion_charge,
 846        adduct_atom=None,
 847    ):
 848        """Find the formulas.
 849
 850        Parameters
 851        ----------
 852        formulas : list of MolecularFormula
 853            The list of molecular formulas.
 854        min_abundance : float
 855            The minimum abundance threshold.
 856        mass_spectrum_obj : MassSpectrum
 857            The mass spectrum object.
 858        ms_peak : MSPeak
 859            The mass spectrum peak.
 860        ion_type : str
 861            The ion type.
 862        ion_charge : int
 863            The ion charge.
 864        adduct_atom : str, optional
 865            The adduct atom, by default None.
 866
 867        Returns
 868        -------
 869        list of MSPeak
 870            The list of mass spectrum peaks with assigned molecular formulas.
 871
 872        Notes
 873        -----
 874        Uses the closest error the next search (this is not ideal, it needs to use confidence
 875        metric to choose the right candidate then propagate the error using the error from the best candidate).
 876        It needs to add s/n to the equation.
 877        It need optimization to define the mz_error_range within a m/z unit since it is directly proportional
 878        with the mass, and inversely proportional to the rp. It's not linear, i.e., sigma mass.
 879        The idea it to correlate sigma to resolving power, signal to noise and sample complexity per mz unit.
 880        Method='distance'
 881        """
 882        mspeak_assigned_index = list()
 883
 884        min_ppm_error = mass_spectrum_obj.molecular_search_settings.min_ppm_error
 885        max_ppm_error = mass_spectrum_obj.molecular_search_settings.max_ppm_error
 886
 887        min_abun_error = mass_spectrum_obj.molecular_search_settings.min_abun_error
 888        max_abun_error = mass_spectrum_obj.molecular_search_settings.max_abun_error
 889
 890        # f = open("abundance_error.txt", "a+")
 891        ms_peak_mz_exp, ms_peak_abundance = ms_peak.mz_exp, ms_peak.abundance
 892        # min_error = min([pmf.mz_error for pmf in possible_formulas])
 893
 894        def mass_by_ion_type(possible_formula_obj):
 895            if ion_type == Labels.protonated_de_ion:
 896                return possible_formula_obj._protonated_mz(ion_charge)
 897
 898            elif ion_type == Labels.radical_ion:
 899                return possible_formula_obj._radical_mz(ion_charge)
 900
 901            elif ion_type == Labels.adduct_ion and adduct_atom:
 902                # Positional argument order differs between formula objects
 903                # and SQL links; keywords are unambiguous.
 904                return possible_formula_obj._adduct_mz(
 905                    ion_charge=ion_charge, adduct_atom=adduct_atom
 906                )
 907
 908            else:
 909                # will return externally calculated mz if is set, #use on Bruker Reference list import
 910                # if the ion type is known the ion mass based on molecular formula ion type
 911                # if ion type is unknow will return neutral mass
 912                return possible_formula_obj.mz_calc
 913
 914        for possible_formula in formulas:
 915            if possible_formula:
 916                error = self.calc_error(
 917                    ms_peak_mz_exp, mass_by_ion_type(possible_formula)
 918                )
 919
 920                # error = possible_formula.mz_error
 921
 922                if min_ppm_error <= error <= max_ppm_error:
 923                    # update the error
 924
 925                    self.set_last_error(error, mass_spectrum_obj)
 926
 927                    # add molecular formula match to ms_peak
 928
 929                    # get molecular formula dict from sql obj
 930                    # formula_dict = pickle.loads(possible_formula.mol_formula)
 931
 932                    formula_dict = possible_formula.to_dict()
 933                    # create the molecular formula obj to be stored
 934                    molecular_formula = MolecularFormula(
 935                        formula_dict,
 936                        ion_charge,
 937                        ion_type=ion_type,
 938                        adduct_atom=adduct_atom,
 939                    )
 940                    # add the molecular formula obj to the mspeak obj
 941                    # add the mspeak obj and it's index for tracking next assignment step
 942
 943                    if self.find_isotopologues:
 944                        # calculates isotopologues
 945                        isotopologues = molecular_formula.isotopologues(
 946                            min_abundance,
 947                            ms_peak_abundance,
 948                            mass_spectrum_obj.dynamic_range,
 949                        )
 950
 951                        # search for isotopologues
 952                        for isotopologue_formula in isotopologues:
 953                            molecular_formula.expected_isotopologues.append(
 954                                isotopologue_formula
 955                            )
 956                            # move this outside to improve preformace
 957                            # we need to increase the search space to -+1 m_z
 958                            first_index, last_index = (
 959                                mass_spectrum_obj.get_nominal_mz_first_last_indexes(
 960                                    isotopologue_formula.mz_nominal_calc
 961                                )
 962                            )
 963
 964                            for ms_peak_iso in mass_spectrum_obj[
 965                                first_index:last_index
 966                            ]:
 967                                error = self.calc_error(
 968                                    ms_peak_iso.mz_exp, isotopologue_formula.mz_calc
 969                                )
 970
 971                                if min_ppm_error <= error <= max_ppm_error:
 972                                    # need to define error distribution for abundance measurements
 973
 974                                    # if mass_spectrum_obj.is_centroid:
 975
 976                                    abundance_error = self.calc_error(
 977                                        isotopologue_formula.abundance_calc,
 978                                        ms_peak_iso.abundance,
 979                                        method="perc",
 980                                    )
 981
 982                                    # area_error = self.calc_error(ms_peak.area, ms_peak_iso.area, method='perc')
 983
 984                                    # margin of error was set empirically/ needs statistical calculation
 985                                    #  of margin of error for the measurement of the abundances
 986                                    if (
 987                                        min_abun_error
 988                                        <= abundance_error
 989                                        <= max_abun_error
 990                                    ):
 991                                        # update the error
 992
 993                                        self.set_last_error(error, mass_spectrum_obj)
 994
 995                                        # isotopologue_formula.mz_error = error
 996
 997                                        # isotopologue_formula.area_error = area_error
 998
 999                                        # isotopologue_formula.abundance_error = abundance_error
1000
1001                                        isotopologue_formula.mspeak_index_mono_isotopic = ms_peak.index
1002
1003                                        mono_isotopic_formula_index = len(ms_peak)
1004
1005                                        isotopologue_formula.mspeak_index_mono_isotopic = ms_peak.index
1006
1007                                        isotopologue_formula.mono_isotopic_formula_index = mono_isotopic_formula_index
1008
1009                                        # add mspeaks isotopologue index to the mono isotopic MolecularFormula obj and the respective formula position
1010
1011                                        # add molecular formula match to ms_peak
1012                                        x = ms_peak_iso.add_molecular_formula(
1013                                            isotopologue_formula
1014                                        )
1015
1016                                        molecular_formula.mspeak_mf_isotopologues_indexes.append(
1017                                            (ms_peak_iso.index, x)
1018                                        )
1019                                        # add mspeaks mono isotopic index to the isotopologue MolecularFormula obj
1020
1021                    y = ms_peak.add_molecular_formula(molecular_formula)
1022
1023                    mspeak_assigned_index.append((ms_peak.index, y))
1024
1025        return mspeak_assigned_index
1026
1027
1028class SearchMolecularFormulasLC:
1029    """Class for searching molecular formulas in a LC object.
1030
1031    Parameters
1032    ----------
1033    lcms_obj : LCMSBase
1034        The LCMSBase object.
1035    sql_db : MolForm_SQL, optional
1036        The SQL database object, by default None.
1037    first_hit : bool, optional
1038        Flag to indicate whether to skip peaks that already have a molecular formula assigned, by default False.
1039    find_isotopologues : bool, optional
1040        Flag to indicate whether to find isotopologues, by default True.
1041
1042    Methods
1043    -------
1044
1045    * search_spectra_against_candidates().
1046        Search a list of mass spectra against a list of candidate formulas with a given ion type, charge, and optional adduct_atom.
1047    * bulk_run_molecular_formula_search().
1048        Run the molecular formula search on the given list of mass spectra.
1049        Pulls the settings from the LCMSBase object to set ion type and charge to search for.
1050    * run_mass_feature_search().
1051        Run the molecular formula search on mass features.
1052        Calls bulk_run_molecular_formula_search() with specified mass spectra and mass peaks.
1053    * run_untargeted_worker_ms1().
1054        Run untargeted molecular formula search on the ms1 mass spectrum.
1055        DEPRECATED: use run_mass_feature_search() or bulk_run_molecular_formula_search() instead.
1056    * run_target_worker_ms1().
1057        Run targeted molecular formula search on the ms1 mass spectrum.
1058        DEPRECATED: use run_mass_feature_search() or bulk_run_molecular_formula_search() instead.
1059    """
1060
1061    def __init__(self, lcms_obj, sql_db=None, first_hit=False, find_isotopologues=True):
1062        self.first_hit = first_hit
1063
1064        self.find_isotopologues = find_isotopologues
1065
1066        self.lcms_obj = lcms_obj
1067
1068        # Catch in-place usedAtoms mutations after settings construction
1069        validate_used_atoms_keys(
1070            self.lcms_obj.parameters.mass_spectrum["ms1"].molecular_search.usedAtoms
1071        )
1072
1073        if not sql_db:
1074            self.sql_db = MolForm_SQL(
1075                url=self.lcms_obj.parameters.mass_spectrum[
1076                    "ms1"
1077                ].molecular_search.url_database
1078            )
1079
1080        else:
1081            self.sql_db = sql_db
1082
1083    def search_spectra_against_candidates(
1084        self,
1085        mass_spectrum_list,
1086        ms_peaks_list,
1087        candidate_formulas,
1088        ion_type,
1089        ion_charge,
1090        adduct_atom=None,
1091    ):
1092        """Search a list of mass spectra against a list of candidate formulas with a given ion type and charge.
1093
1094        Parameters
1095        ----------
1096        mass_spectrum_list : list of MassSpectrum
1097            The list of mass spectra to perform the search on.
1098        ms_peaks_list : list of lists of MSPeak objects
1099            The list of mass spectrum peaks to search within each mass spectrum.
1100        candidate_formulas : dict
1101            The candidate formulas.
1102        ion_type : str
1103            The ion type.
1104        ion_charge : int
1105            signed ion charge used for the candidate formulas.
1106        adduct_atom : str, optional
1107            The adduct atom, by default None.
1108
1109        Notes
1110        -----
1111        This function is designed to be used with the bulk_run_molecular_formula_search function.
1112        """
1113        for mass_spectrum, ms_peaks in zip(mass_spectrum_list, ms_peaks_list):
1114            single_ms_search = SearchMolecularFormulas(
1115                mass_spectrum,
1116                sql_db=self.sql_db,
1117                first_hit=self.first_hit,
1118                find_isotopologues=self.find_isotopologues,
1119            )
1120            single_ms_search.run_search(
1121                ms_peaks,
1122                candidate_formulas,
1123                mass_spectrum.min_abundance,
1124                ion_type,
1125                ion_charge,
1126                adduct_atom=adduct_atom,
1127            )
1128
1129    def bulk_run_molecular_formula_search(
1130        self, mass_spectrum_list, ms_peaks_list, mass_spectrum_setting_key="ms1"
1131    ):
1132        """Run the molecular formula search on the given list of mass spectra.
1133
1134        Uses ``min_ion_charge`` / ``max_ion_charge`` from molecular search settings
1135        to search each peak at every signed charge in that absolute range
1136        (defaults 1..1 = single-charge). Does not determine or rewrite peak
1137        ``polarity``; assignment charge is on each ``MolecularFormula``.
1138
1139        Parameters
1140        ----------
1141        mass_spectrum_list : list of MassSpectrum
1142            The list of mass spectra to search.
1143        ms_peaks_list : list of lists of MSPeak objects
1144            The mass peaks to perform molecular formula search within each mass spectrum
1145        mass_spectrum_setting_key : str, optional
1146            The mass spectrum setting key, by default 'ms1'.
1147            This is used to get the appropriate molecular search settings from the LCMSBase object
1148        """
1149        if self.lcms_obj.polarity not in ("positive", "negative"):
1150            raise ValueError("Polarity must be either 'positive' or 'negative'")
1151
1152        # Check that the length of the mass spectrum list and the ms_peaks list are the same
1153        if len(mass_spectrum_list) != len(ms_peaks_list):
1154            raise ValueError(
1155                "The length of the mass spectrum list and the ms_peaks list must be the same"
1156            )
1157
1158        mol_settings = self.lcms_obj.parameters.mass_spectrum[
1159            mass_spectrum_setting_key
1160        ].molecular_search
1161        # Re-check after possible post-construction mutation of settings fields
1162        mol_settings.validate_ion_charge_settings()
1163        search_charges = SearchMolecularFormulas.ion_charges_for_search(
1164            self.lcms_obj.polarity,
1165            min_ion_charge=mol_settings.min_ion_charge,
1166            max_ion_charge=mol_settings.max_ion_charge,
1167        )
1168
1169        nominal_mzs = [x.nominal_mz for x in mass_spectrum_list]
1170        nominal_mzs = list(set([item for sublist in nominal_mzs for item in sublist]))
1171
1172        # reset average error, only relevant if average mass error method is being used
1173        SearchMolecularFormulaWorker(
1174            find_isotopologues=self.find_isotopologues
1175        ).reset_error(mass_spectrum_list[0])
1176
1177        # check database for all possible molecular formula combinations based on the setting passed to self.mass_spectrum_obj.molecular_search_settings
1178        classes = MolecularCombinations(self.sql_db).runworker(
1179            mol_settings, print_time=mol_settings.verbose_processing
1180        )
1181
1182        try:
1183
1184            def apply_job(job):
1185                self.search_spectra_against_candidates(
1186                    mass_spectrum_list=mass_spectrum_list,
1187                    ms_peaks_list=ms_peaks_list,
1188                    candidate_formulas=job.candidate_formulas,
1189                    ion_type=job.ion_type,
1190                    ion_charge=job.ion_charge,
1191                    adduct_atom=job.adduct_atom,
1192                )
1193
1194            SearchMolecularFormulas._run_formula_search_jobs(
1195                classes,
1196                nominal_mzs,
1197                mol_settings,
1198                search_charges,
1199                self.sql_db,
1200                apply_job,
1201            )
1202        finally:
1203            self.sql_db.close()
1204
1205    def run_mass_feature_search(self):
1206        """Run the molecular formula search on the mass features.
1207
1208        Calls bulk_run_molecular_formula_search() with specified mass spectra and mass peaks.
1209        Multi-charge range comes from molecular search ``min_ion_charge`` /
1210        ``max_ion_charge``. Peak ``polarity`` is not rewritten.
1211        """
1212        mass_features_df = self.lcms_obj.mass_features_to_df()
1213
1214        # Get the list of mass spectrum (and peaks to search with each mass spectrum) for all mass features
1215        scan_list = mass_features_df.apex_scan.unique()
1216        mass_spectrum_list = [self.lcms_obj._ms[x] for x in scan_list]
1217        ms_peaks = []
1218        for scan in scan_list:
1219            mf_df_scan = mass_features_df[mass_features_df.apex_scan == scan]
1220            peaks_to_search = [
1221                self.lcms_obj.mass_features[x].ms1_peak
1222                for x in mf_df_scan.index.tolist()
1223            ]
1224            ms_peaks.append(peaks_to_search)
1225
1226        # Run the molecular formula search
1227        self.bulk_run_molecular_formula_search(mass_spectrum_list, ms_peaks)
1228
1229    def run_untargeted_worker_ms1(self):
1230        """Run untargeted molecular formula search on the ms1 mass spectrum."""
1231        raise NotImplementedError(
1232            "run_untargeted_worker_ms1 search is not implemented in CoreMS 3.0 and greater"
1233        )
1234
1235    def run_target_worker_ms1(self):
1236        """Run targeted molecular formula search on the ms1 mass spectrum."""
1237        raise NotImplementedError(
1238            "run_target_worker_ms1 formula search is not yet implemented in CoreMS 3.0 and greater"
1239        )