corems.molecular_formula.calc.MolecularFormulaCalc

   1__author__ = "Yuri E. Corilo"
   2__date__ = "Jun 24, 2019"
   3
   4
   5import re
   6import warnings
   7
   8import IsoSpecPy
   9from numpy import array, exp, isnan, nextafter, power
  10
  11# this is to handle both versions of IsoSpecPy, 2.0.2 and 2.2.2
  12# TODO in a future release remove support for legacy isospecpy
  13from packaging import version
  14
  15from corems.encapsulation.constant import (
  16    ADDUCT_ALIASES as _ADDUCT_ALIASES,
  17    Atoms,
  18    ION_TYPE_DICT,
  19    Labels,
  20)
  21from corems.encapsulation.factory.parameters import MSParameters
  22from corems.molecular_id.calc.SpectralSimilarity import SpectralSimilarity
  23
  24isospec_version = IsoSpecPy.__version__
  25if version.parse(isospec_version) > version.parse("2.0.2"):
  26    legacy_isospec = False
  27else:
  28    legacy_isospec = True
  29if legacy_isospec:
  30    from IsoSpecPy import IsoSpecPy
  31
  32    warnings.warn(
  33        f"IsoSpecPy version {isospec_version} is installed, and support is deprecated. Please update to 2.2.2",
  34        DeprecationWarning,
  35    )
  36
  37
  38class MolecularFormulaCalc:
  39    """Class of calculations related to molecular formula
  40
  41    This class is not intended to be used directly, but rather to be inherited by other classes in the molecular_formula/factory module like MolecularFormula and MolecularFormulaIsotopologue
  42
  43    Attributes
  44    ----------
  45    mz_calc : float
  46        The m/z value of the molecular formula.
  47    neutral_mass : float
  48        The neutral mass of the molecular formula.
  49    ion_charge : int
  50        The ion charge of the molecular formula.
  51    _external_mz : float
  52        The externally provided m/z value of the molecular formula.
  53    _d_molecular_formula : dict
  54        The dictionary representation of the molecular formula.
  55    _mspeak_parent : object
  56        The parent MS peak object associated with the molecular formula.
  57    _assignment_mass_error : float
  58        The mass error of the molecular formula.
  59
  60    Methods
  61    -------
  62    * _calc_resolving_power_low_pressure(B, T)
  63        Calculate the resolving power at low pressure.
  64    * _calc_resolving_power_high_pressure(B, T)
  65        Calculate the resolving power at high pressure.
  66    * _adduct_mz(adduct_atom, ion_charge)
  67        Get the m/z value of an adducted ion version of the molecular formula.
  68    * _protonated_mz(ion_charge)
  69        Get the m/z value of a protonated or deprotonated ion version of the molecular formula.
  70    * _radical_mz(ion_charge)
  71        Get the m/z value of a radical ion version of the molecular formula.
  72    * normalize_ion_type(ion_type)
  73        Map alternate LCMS adduct labels onto ``ion_type_dict`` keys.
  74    * charge_from_adduct(adduct)
  75        Parse integer charge from an LCMS adduct string.
  76    * get_ion_formula(neutral_formula, ion_type)
  77        Build ion formula string from neutral formula and LCMS adduct / ion type.
  78    * precursor_mz_from_formula(neutral_formula, ion_type, charge=None)
  79        Calculate precursor m/z from neutral formula and LCMS adduct / ion type.
  80    * _neutral_mass()
  81        Get the neutral mass of the molecular formula.
  82    * _calc_mz()
  83        Get the m/z value of the molecular formula.
  84    * _calc_assignment_mass_error(method='ppm')
  85        Calculate the mass error of the molecular formula.
  86    * _calc_mz_confidence(mean=0)
  87        Calculate the m/z confidence of the molecular formula.
  88    * _calc_isotopologue_confidence()
  89        Calculate the isotopologue confidence of the molecular formula.
  90    * normalize_distance(dist, dist_range)
  91        Normalize the distance value.
  92    * subtract_formula(formula_obj, formated=True)
  93        Subtract a formula from the current formula object.
  94    * _calc_average_mz_score()
  95        Calculate the average m/z error score of the molecular formula identification, including the isotopologues.
  96
  97    Notes
  98    -----
  99    LCMS adduct helpers use :data:`~corems.encapsulation.constant.ION_TYPE_DICT`
 100    (``[M+H]+`` / ``[M-H]-``, etc.). Distinct from FTMS assignment labels in
 101    :class:`~corems.encapsulation.constant.Labels` (``DE_OR_PROTONATED``,
 102    ``RADICAL``, ``ADDUCT``) used by ``_calc_mz``.
 103    """
 104
 105    # Back-compat aliases → corems.encapsulation.constant source of truth
 106    ion_type_dict = ION_TYPE_DICT
 107    ADDUCT_ALIASES = _ADDUCT_ALIASES
 108
 109    def _calc_resolving_power_low_pressure(self, B, T):
 110        """
 111        Calculate the resolving power at low pressure.
 112
 113        Parameters
 114        ----------
 115        B : float
 116            Magnetic Strength (Testa).
 117        T : float
 118            Transient time (seconds).
 119
 120        """
 121        return (1.274 * 10000000 * B * T) * (1 / self.mz_calc)
 122
 123    def _calc_resolving_power_high_pressure(self, B, T):
 124        """
 125        Calculate the resolving power at high pressure.
 126
 127        Parameters
 128        ----------
 129        B : float
 130            Magnetic Strength (Testa).
 131        T : float
 132            Transient time (seconds).
 133
 134        """
 135        return (2.758 * 10000000 * B * T) * (1 / self.mz_calc)
 136
 137    def _adduct_mz(self, adduct_atom, ion_charge):
 138        """Get the m/z value of an adducted ion version of the molecular formula.
 139
 140        Parameters
 141        ----------
 142        adduct_atom : str
 143            The adduct atom.
 144        ion_charge : int
 145            The ion charge.
 146
 147        """
 148        return (
 149            self.neutral_mass
 150            + (Atoms.atomic_masses.get(adduct_atom))
 151            + (ion_charge * -1 * Atoms.electron_mass)
 152        ) / abs(ion_charge)
 153
 154    def _protonated_mz(self, ion_charge):
 155        """Get the m/z value of a protonated or deprotonated ion version of the molecular formula.
 156
 157        Parameters
 158        ----------
 159        ion_charge : int
 160            The ion charge.
 161        """
 162        return (
 163            self.neutral_mass
 164            + (ion_charge * Atoms.atomic_masses.get("H"))
 165            + (ion_charge * -1 * Atoms.electron_mass)
 166        ) / abs(ion_charge)
 167
 168    def _radical_mz(self, ion_charge):
 169        """Get the m/z value of a radical ion version of the molecular formula.
 170
 171        Parameters
 172        ----------
 173        ion_charge : int
 174            The ion charge.
 175        """
 176        return (self.neutral_mass + (ion_charge * -1 * Atoms.electron_mass)) / abs(
 177            ion_charge
 178        )
 179
 180    @staticmethod
 181    def normalize_ion_type(ion_type: str) -> str:
 182        """Map alternate LCMS adduct labels onto ``ion_type_dict`` keys."""
 183        return MolecularFormulaCalc.ADDUCT_ALIASES.get(ion_type, ion_type)
 184
 185    @staticmethod
 186    def charge_from_adduct(adduct: str) -> int:
 187        """Parse integer charge from an adduct string (e.g. ``[M+H]+``, ``[M+2H]2+``)."""
 188        a = (adduct or "").strip()
 189        m = re.search(r"(\d+)([+-])$", a)
 190        if m:
 191            return int(m.group(1)) * (1 if m.group(2) == "+" else -1)
 192        if a.endswith("+"):
 193            return 1
 194        if a.endswith("-"):
 195            return -1
 196        raise ValueError(f"cannot parse charge from adduct {adduct!r}")
 197
 198    @staticmethod
 199    def get_ion_formula(neutral_formula, ion_type):
 200        """From a neutral formula and an ion type, return the formula of the ion.
 201
 202        Parameters
 203        ----------
 204        neutral_formula : str
 205            Neutral formula as MolecularFormula-style (``'C2 H4 O2'``) or compact
 206            (``'C2H4O2'``; isotopes only handled in the spaced form).
 207        ion_type : str
 208            Ion type / adduct key (see ``ion_type_dict``), or a known alias.
 209            LCMS adduct notation (e.g. ``[M+H]+``), not FTMS ``Labels`` ion types.
 210
 211        Returns
 212        -------
 213        str or None
 214            Ion formula string (e.g. ``'C2 H5 O2'``), or None if
 215            ``neutral_formula`` is not a string.
 216        """
 217        # Lazy import avoids circular dependency: factory inherits this class.
 218        from corems.molecular_formula.factory.MolecularFormulaFactory import (
 219            MolecularFormula,
 220        )
 221
 222        if not isinstance(neutral_formula, str):
 223            return None
 224
 225        ion_type = MolecularFormulaCalc.normalize_ion_type(ion_type)
 226        if ion_type not in MolecularFormulaCalc.ion_type_dict:
 227            raise KeyError(f"unsupported ion type {ion_type!r}")
 228
 229        if re.search(r"\s", neutral_formula):
 230            formula_obj = MolecularFormula(neutral_formula, ion_charge=0)
 231        else:
 232            form_pre = re.sub(r"([A-Z])", r" \1", neutral_formula)[1:]
 233            elements = [re.findall(r"[A-Z][a-z]*", x) for x in form_pre.split()]
 234            counts = [re.findall(r"\d+", x) for x in form_pre.split()]
 235            formula_obj = MolecularFormula(
 236                dict(
 237                    zip(
 238                        [x[0] for x in elements],
 239                        [int(x[0]) if x else 1 for x in counts],
 240                    )
 241                ),
 242                ion_charge=0,
 243            )
 244        neutral_formula_dict = formula_obj.to_dict().copy()
 245
 246        entry = MolecularFormulaCalc.ion_type_dict[ion_type]
 247        adduct_add_dict = entry["add"]
 248        for key in adduct_add_dict:
 249            if key in neutral_formula_dict:
 250                neutral_formula_dict[key] += adduct_add_dict[key]
 251            else:
 252                neutral_formula_dict[key] = adduct_add_dict[key]
 253
 254        adduct_subtract = entry["sub"]
 255        for key in adduct_subtract:
 256            neutral_formula_dict[key] -= adduct_subtract[key]
 257
 258        return MolecularFormula(neutral_formula_dict, ion_charge=0).string
 259
 260    @staticmethod
 261    def polarity_from_ion_type(ion_type: str) -> str:
 262        """Return ``'positive'`` or ``'negative'`` from ``ION_TYPE_DICT``."""
 263        ion_type = MolecularFormulaCalc.normalize_ion_type(ion_type)
 264        if ion_type not in ION_TYPE_DICT:
 265            raise KeyError(f"unsupported ion type {ion_type!r}")
 266        return ION_TYPE_DICT[ion_type]["polarity"]
 267
 268    @staticmethod
 269    def precursor_mz_from_formula(neutral_formula, ion_type, charge=None):
 270        """Calculated precursor m/z from neutral formula and adduct / ion type.
 271
 272        Parameters
 273        ----------
 274        neutral_formula : str
 275            Neutral molecular formula.
 276        ion_type : str
 277            Adduct or ion type (``ion_type_dict`` key or alias).
 278            LCMS adduct notation (e.g. ``[M-H]-``), not FTMS ``Labels`` ion types.
 279        charge : int, optional
 280            Ion charge. If omitted, inferred from ``ion_type``.
 281
 282        Returns
 283        -------
 284        float
 285            Precursor m/z of the ion.
 286
 287        Raises
 288        ------
 289        KeyError
 290            If ``ion_type`` is not supported.
 291        ValueError
 292            If the ion formula cannot be built or charge cannot be inferred.
 293        """
 294        # Lazy import avoids circular dependency: factory inherits this class.
 295        from corems.molecular_formula.factory.MolecularFormulaFactory import (
 296            MolecularFormula,
 297        )
 298
 299        ion_form = MolecularFormulaCalc.get_ion_formula(neutral_formula, ion_type)
 300        if ion_form is None:
 301            raise ValueError(
 302                f"could not build ion formula for {neutral_formula!r} {ion_type!r}"
 303            )
 304        if charge is None:
 305            charge = MolecularFormulaCalc.charge_from_adduct(
 306                MolecularFormulaCalc.normalize_ion_type(ion_type)
 307            )
 308        if charge == 0:
 309            raise ValueError("charge must be non-zero for precursor m/z")
 310        mf = MolecularFormula(ion_form, ion_charge=0)
 311        return float(
 312            (mf.neutral_mass + (charge * -1 * Atoms.electron_mass)) / abs(charge)
 313        )
 314
 315    def _neutral_mass(self):
 316        """Get the neutral mass of the molecular formula."""
 317
 318        mass = 0
 319
 320        for each_atom in self._d_molecular_formula.keys():
 321            if each_atom != Labels.ion_type and each_atom != "HC":
 322                try:
 323                    mass = mass + Atoms.atomic_masses[
 324                        each_atom
 325                    ] * self._d_molecular_formula.get(each_atom)
 326
 327                except:
 328                    print(Labels.ion_type, each_atom)
 329
 330        return mass
 331
 332    def _calc_mz(self):
 333        """Get the m/z value of the molecular formula, based on the ion charge and ion type."""
 334
 335        if self.ion_charge is not None:
 336            if self._external_mz:
 337                return self._external_mz
 338
 339            else:
 340                ion_type = self._d_molecular_formula.get(Labels.ion_type)
 341
 342                if ion_type == Labels.protonated_de_ion:
 343                    return self.protonated_mz
 344
 345                elif ion_type == Labels.radical_ion or ion_type == Labels.adduct_ion:
 346                    return self.radical_mz
 347
 348                elif ion_type == Labels.neutral:
 349                    return self.neutral_mass
 350
 351                elif self.ion_charge == 0:
 352                    return self.neutral_mass
 353
 354                else:
 355                    # formula is probably ion form used for bruker ref list
 356                    return self.neutral_mass
 357
 358        else:
 359            raise Exception("Please set ion charge first")
 360
 361    def _calc_assignment_mass_error(self, method="ppm"):
 362        """Calculate the mass error of the molecular formula, based on the experimental m/z and the calculated m/z.
 363
 364        Parameters
 365        ----------
 366        method : str, optional
 367            The method to calculate the mass error, by default 'ppm', but can be 'ppb'
 368
 369        Raises
 370        ------
 371        Exception
 372            If the method is not 'ppm' or 'ppb'.
 373        Exception
 374            If there is no ms peak associated with the molecular formula instance.
 375        """
 376
 377        if method == "ppm":
 378            multi_factor = 1000000
 379
 380        elif method == "ppb":
 381            multi_factor = 1000000
 382
 383        else:
 384            raise Exception(
 385                "method needs to be ppm or ppb, you have entered %s" % method
 386            )
 387
 388        if self._mspeak_parent.mz_exp:
 389            self._assignment_mass_error = (
 390                (self._mspeak_parent.mz_exp - self.mz_calc) / self.mz_calc
 391            ) * multi_factor
 392
 393            return (
 394                (self._mspeak_parent.mz_exp - self.mz_calc) / self.mz_calc
 395            ) * multi_factor
 396
 397        else:
 398            raise Exception(
 399                "No ms peak associated with the molecular formula instance %s", self
 400            )
 401
 402    def _calc_mz_confidence(self, mean=0):
 403        """Calculate the monoisotopic peak Gaussian mass-accuracy score.
 404
 405        Implementation for the public :attr:`mz_error_score` property.
 406        See that property and the ``corems.molecular_formula`` package
 407        documentation for the equation, ``predicted_std`` / 1.66 ppm
 408        fallback, and interpretation.
 409
 410        Parameters
 411        ----------
 412        mean : float, optional
 413            Mean mass error (ppm) for the Gaussian; default 0 (calibrated
 414            spectrum with near-zero mean error).
 415        """
 416
 417        # predicted_std not set: fallback width (ppm) for the Gaussian score
 418        if not self._mspeak_parent.predicted_std:
 419            self._mspeak_parent.predicted_std = 1.66
 420
 421        return exp(
 422            -1
 423            * (
 424                power((self.mz_error - mean), 2)
 425                / (2 * power(self._mspeak_parent.predicted_std, 2))
 426            )
 427        )
 428
 429    def _calc_isotopologue_confidence(self):
 430        """Calculate isotopologue abundance pattern similarity.
 431
 432        Implementation for the public :attr:`isotopologue_similarity`
 433        property. See that property and the ``corems.molecular_formula``
 434        package documentation for interpretation.
 435
 436        Returns
 437        -------
 438        float
 439            Similarity in about [0, 1], or 0.0 when no isotopologues are
 440            expected.
 441        """
 442
 443        if self.is_isotopologue:
 444            # confidence of isotopologue is pure mz error
 445            # TODO add more features here
 446
 447            mformula_index = self.mono_isotopic_formula_index
 448            mspeak_index = self.mspeak_index_mono_isotopic
 449            parent_ms = self._mspeak_parent._ms_parent
 450            if (
 451                mformula_index is None
 452                or mspeak_index is None
 453                or mspeak_index < 0
 454                or mspeak_index >= len(parent_ms)
 455            ):
 456                return 0.0
 457            mspeak = parent_ms[mspeak_index]
 458            if mformula_index < 0 or mformula_index >= len(mspeak):
 459                # Parent formula list was filtered (e.g. ion-type constrain);
 460                # stale isotopologue indexes should not fail export.
 461                return 0.0
 462
 463            expected_isotopologues = mspeak[mformula_index].expected_isotopologues
 464
 465            mono_mz = mspeak[mformula_index].mz_calc
 466            mono_abundance = mspeak.abundance
 467
 468        else:
 469            mono_mz = self.mz_calc
 470            mono_abundance = self._mspeak_parent.abundance
 471
 472            expected_isotopologues = self.expected_isotopologues
 473            # has isotopologues based on current dinamic range
 474
 475        if expected_isotopologues:
 476            dict_mz_abund_ref = {"mz": [mono_mz], "abundance": [mono_abundance]}
 477
 478            # get reference data
 479            for mf in expected_isotopologues:
 480                dict_mz_abund_ref["abundance"].append(mf.abundance_calc)
 481                dict_mz_abund_ref["mz"].append(mf.mz_calc)
 482
 483            dict_mz_abund_exp = {mono_mz: mono_abundance}
 484
 485            # get experimental data
 486            for mf in expected_isotopologues:
 487                # molecular formula has been assigned to a peak
 488                if mf._mspeak_parent:
 489                    # stores mspeak abundance
 490                    dict_mz_abund_exp[mf.mz_calc] = mf._mspeak_parent.abundance
 491
 492                else:
 493                    # fill missing mz with abundance 0 and mz error score of 0
 494                    dict_mz_abund_exp[mf.mz_calc] = nextafter(0, 1)
 495
 496            distance = SpectralSimilarity(
 497                dict_mz_abund_exp, dict_mz_abund_ref
 498            ).manhattan_distance()
 499            correlation = 1 - self.normalize_distance(distance, [0, 2])
 500            # correlation = dwt_correlation(dict_mz_abund_exp, dict_mz_abund_ref)
 501            # correlation = cosine_correlation(dict_mz_abund_exp, dict_mz_abund_ref)
 502
 503            if isnan(correlation):
 504                correlation = 0.00001
 505
 506        else:
 507            # no isotopologue expected giving a correlation score of 0.0 but it needs optimization
 508            correlation = 0.0
 509
 510        return correlation
 511
 512    def normalize_distance(self, dist, dist_range):
 513        """
 514        Normalize the distance value.
 515
 516        Parameters
 517        ----------
 518        dist : float
 519            The distance value to be normalized.
 520        dist_range : list
 521            The range of the distance value.
 522
 523        """
 524        result = (dist - dist_range[0]) / (dist_range[1] - dist_range[0])
 525
 526        if result < 0:
 527            result = 0.0
 528        elif result > 1:
 529            result = 1.0
 530
 531        return result
 532
 533    def subtract_formula(self, formula_obj, formated=True):
 534        """Subtract a formula from the current formula object
 535
 536        Parameters
 537        ----------
 538        formula_obj : MolecularFormula
 539            MolecularFormula object to be subtracted from the current formula object
 540        formated : bool, optional
 541            If True, returns the formula in string format, by default True
 542
 543        """
 544        subtraction = {}
 545        for atom, value in self.to_dict().items():
 546            if atom != Labels.ion_type:
 547                if formula_obj.get(atom):
 548                    # value_subtraction = value - formula_obj.get(atom)
 549                    if value - formula_obj.get(atom) > 0:
 550                        subtraction[atom] = value - formula_obj.get(atom)
 551                else:
 552                    subtraction[atom] = value
 553        if formated:
 554            SUB = str.maketrans("0123456789", "₀₁₂₃₄₅₆₇₈₉")
 555            SUP = str.maketrans("0123456789", "⁰¹²³⁴⁵⁶⁷⁸⁹")
 556        else:
 557            SUB = str.maketrans("0123456789", "0123456789")
 558            SUP = str.maketrans("0123456789", "0123456789")
 559        formula_srt = ""
 560        for atom in Atoms.atoms_order:
 561            if atom in subtraction.keys():
 562                formula_srt += atom.translate(SUP) + str(
 563                    int(subtraction.get(atom))
 564                ).translate(SUB)
 565
 566        return formula_srt
 567
 568    def _calc_average_mz_score(self):
 569        """Average mass-accuracy score over mono and expected isotopologues.
 570
 571        Implementation for the public :attr:`average_mz_error_score`
 572        property. See that property and the ``corems.molecular_formula``
 573        package documentation for interpretation.
 574        """
 575        if self.is_isotopologue:
 576            # confidence of isotopologue is pure mz error
 577            # TODO add more features here
 578
 579            mformula_index = self.mono_isotopic_formula_index
 580            mspeak_index = self.mspeak_index_mono_isotopic
 581
 582            mspeak = self._mspeak_parent._ms_parent[mspeak_index]
 583
 584            expected_isotopologues = mspeak[mformula_index].expected_isotopologues
 585
 586        else:
 587            expected_isotopologues = self.expected_isotopologues
 588            # has isotopologues based on current dinamic range
 589
 590        accumulated_mz_score = [self.mz_error_score]
 591
 592        if expected_isotopologues:
 593            for mf in expected_isotopologues:
 594                # molecular formula has been assigned to a peak
 595                if mf._mspeak_parent:
 596                    # stores mspeak abundance
 597                    accumulated_mz_score.append(mf.mz_error_score)
 598                else:
 599                    # fill missing mz with abundance 0 and mz error score of 0
 600                    accumulated_mz_score.append(0.0)
 601
 602        average_mz_score = sum(accumulated_mz_score) / len(accumulated_mz_score)
 603
 604        if isnan(average_mz_score):
 605            average_mz_score = 0.0
 606
 607        return average_mz_score
 608
 609    def _calc_confidence_score(self):
 610        """Composite confidence score from mass error and isotopologues.
 611
 612        Implementation for the public :attr:`confidence_score` property:
 613
 614        ```
 615        CS = (w_err * average_mz_error_score) + (w_iso * isotopologue_similarity)
 616        ```
 617
 618        with weights from the parent spectrum's molecular search settings.
 619        Assumes a calibrated spectrum with near-zero mean mass error.
 620        See the public property and the ``corems.molecular_formula`` package
 621        documentation for the full formulation and literature reference.
 622        """
 623
 624        isotopologue_correlation = self.isotopologue_similarity
 625        average_mz_score = self.average_mz_error_score
 626        a = self._mspeak_parent._ms_parent.molecular_search_settings.mz_error_score_weight
 627        b = self._mspeak_parent._ms_parent.molecular_search_settings.isotopologue_score_weight
 628
 629        score = (isotopologue_correlation * b) + (average_mz_score * a)
 630
 631        return score
 632
 633    def _calc_abundance_error(self, method="percentile"):
 634        """Calculate the abundance error of the molecular formula, based on the experimental abundance and the calculated abundance.
 635
 636        Parameters
 637        ----------
 638        method : str, optional
 639            The method to calculate the abundance error, by default 'percentile', but can be 'ppm' or 'ppb'
 640
 641        Returns
 642        -------
 643        float
 644            The abundance error of the molecular formula.
 645
 646        Raises
 647        ------
 648        Exception
 649            If isotopologues were not calculated.
 650        """
 651
 652        mult_factor = 100
 653
 654        iso_abundance = self._mspeak_parent.abundance
 655        mono_abundance = self._mspeak_parent._ms_parent[
 656            self.mspeak_index_mono_isotopic
 657        ].abundance
 658
 659        if self.prob_ratio:
 660            theor_abundance = mono_abundance * self.prob_ratio
 661            # self.parent need to have a MassSpecPeak associated with the MolecularFormula class
 662            return ((theor_abundance - iso_abundance) / theor_abundance) * mult_factor
 663
 664        else:
 665            raise Exception("Please calc_isotopologues")
 666
 667    def _calc_area_error(self, method="percentile"):
 668        """Calculate the area error of the molecular formula, based on the experimental area and the calculated area.
 669
 670        Parameters
 671        ----------
 672        method : str, optional
 673            The method to calculate the area error, by default 'percentile', but can be 'ppm' or 'ppb'
 674
 675        Returns
 676        -------
 677        float
 678            The area error of the molecular formula.
 679
 680        Raises
 681        ------
 682        Exception
 683            If isotopologues were not calculated.
 684        """
 685
 686        mult_factor = 100
 687
 688        iso_area = self._mspeak_parent.area
 689        mono_area = self._mspeak_parent._ms_parent[self.mspeak_index_mono_isotopic].area
 690
 691        if self.prob_ratio:
 692            if mono_area and iso_area:
 693                # exp_ratio = iso_area/mono_area
 694
 695                area_calc = mono_area * self.prob_ratio
 696
 697                # self.parent need to have a MassSpecPeak associated with the MolecularFormula class
 698                return ((area_calc - iso_area) / area_calc) * mult_factor
 699                # return ((self.prob_ratio - exp_ratio )/self.prob_ratio)*mult_factor
 700
 701            else:
 702                # centroid mass spectrum
 703                return 0
 704        else:
 705            raise Exception("Please calc_isotopologues")
 706
 707    def _calc_aromaticity_index_mod(self):
 708        """Calculate the modified aromaticity index of the molecular formula.
 709
 710        Returns
 711        -------
 712        float
 713            The aromaticity index of the molecular formula.
 714
 715        Notes
 716        -----
 717        Source Koch and Dittmar, 2006 https://doi.org/10.1002/rcm.2386
 718        corrected in https://doi.org/10.1002/rcm.7433
 719        """
 720        # Prepare empty dictionary to store the number of atoms of each element
 721        ai_es = {"C": 0, "H": 0, "O": 0, "N": 0, "S": 0}
 722
 723        # Count the number of atoms of each element in the molecular formula, inclusive of isotopes
 724        for element in ai_es:
 725            elements_w_iso = [element] + Atoms.isotopes.get(element)[1]
 726            for element_w_iso in elements_w_iso:
 727                if element_w_iso in self._d_molecular_formula:
 728                    ai_es[element] += self._d_molecular_formula[element_w_iso]
 729
 730        ai_n = (
 731            1
 732            + ai_es["C"]
 733            - (0.5 * ai_es["O"])
 734            - ai_es["S"]
 735            - (0.5 * (ai_es["N"] + ai_es["H"]))
 736        )
 737        ai_d = ai_es["C"] - (0.5 * ai_es["O"]) - ai_es["N"] - ai_es["S"]
 738
 739        ai = ai_n / ai_d
 740
 741        if ai < 0:
 742            ai = 0
 743        if ai > 1:
 744            ai = 1
 745
 746        return ai
 747
 748    def _calc_aromaticity_index(self):
 749        """Calculate the aromaticity index of the molecular formula.
 750
 751        Returns
 752        -------
 753        float
 754            The aromaticity index of the molecular formula.
 755
 756        Notes
 757        -----
 758        Source Koch and Dittmar, 2006 https://doi.org/10.1002/rcm.2386
 759        corrected in https://doi.org/10.1002/rcm.7433
 760        """
 761        # Prepare empty dictionary to store the number of atoms of each element
 762        ai_es = {"C": 0, "H": 0, "O": 0, "N": 0, "S": 0}
 763
 764        # Count the number of atoms of each element in the molecular formula, inclusive of isotopes
 765        for element in ai_es:
 766            elements_w_iso = [element] + Atoms.isotopes.get(element)[1]
 767            for element_w_iso in elements_w_iso:
 768                if element_w_iso in self._d_molecular_formula:
 769                    ai_es[element] += self._d_molecular_formula[element_w_iso]
 770
 771                ai_n = (
 772                    1
 773                    + ai_es["C"]
 774                    - (ai_es["O"])
 775                    - ai_es["S"]
 776                    - (0.5 * (ai_es["N"] + ai_es["H"]))
 777                )
 778        ai_d = ai_es["C"] - (ai_es["O"]) - ai_es["N"] - ai_es["S"]
 779
 780        ai = ai_n / ai_d
 781
 782        if ai < 0:
 783            ai = 0
 784        if ai > 1:
 785            ai = 1
 786
 787        return ai
 788
 789    def _calc_nosc(self):
 790        """Calculate the average nominal oxidation state of carbon
 791
 792        Returns
 793        -------
 794        float
 795            The average nominal oxidation state of carbon
 796
 797        Notes
 798        -----
 799        Source LaRowe and Van Cappellen, 2011 https://doi.org/10.1016/j.gca.2011.01.020
 800        """
 801        # Prepare empty dictionary to store the number of atoms of each element
 802        nosc_es = {"C": 0, "H": 0, "O": 0, "N": 0, "S": 0, "P": 0}
 803
 804        # Count the number of atoms of each element in the molecular formula, inclusive of isotopes
 805        for element in nosc_es:
 806            elements_w_iso = [element] + Atoms.isotopes.get(element)[1]
 807            for element_w_iso in elements_w_iso:
 808                if element_w_iso in self._d_molecular_formula:
 809                    nosc_es[element] += self._d_molecular_formula[element_w_iso]
 810
 811        nosc = (
 812            -(
 813                (
 814                    4 * nosc_es["C"]
 815                    + nosc_es["H"]
 816                    - 3 * nosc_es["N"]
 817                    - 2 * nosc_es["O"]
 818                    + 5 * nosc_es["P"]
 819                    - 2 * nosc_es["S"]
 820                )
 821                / nosc_es["C"]
 822            )
 823            + 4
 824        )
 825
 826        # If nosc is infinite or negative infinity, set it to nan
 827        if nosc == float("inf") or nosc == float("-inf"):
 828            nosc = float("nan")
 829
 830        return nosc
 831
 832    @property
 833    def dbe_ai(self):
 834        """Calculate the double bond equivalent (DBE) of the molecular formula, based on the number of carbons, hydrogens, and oxygens."""
 835
 836        carbons = self._d_molecular_formula.get("C")
 837        hydrogens = self._d_molecular_formula.get("H")
 838        oxygens = self._d_molecular_formula.get("O")
 839        return 1 + (((2 * carbons) - hydrogens - (2 * oxygens)) * 0.5)
 840
 841    def _calc_dbe(self):
 842        """Calculate the double bond equivalent (DBE) of the molecular formula"""
 843
 844        individual_dbe = 0
 845
 846        for atom in self._d_molecular_formula.keys():
 847            if atom != Labels.ion_type:
 848                n_atom = int(self._d_molecular_formula.get(atom))
 849
 850                clean_atom = "".join([i for i in atom if not i.isdigit()])
 851
 852                if self._mspeak_parent:
 853                    valencia = self._mspeak_parent._ms_parent.molecular_search_settings.used_atom_valences.get(
 854                        clean_atom
 855                    )
 856                else:
 857                    valencia = MSParameters.molecular_search.used_atom_valences.get(
 858                        clean_atom
 859                    )
 860                # valencia = Atoms.atoms_covalence.get(atom)
 861
 862                if type(valencia) is tuple:
 863                    valencia = valencia[0]
 864                if valencia > 0:
 865                    # print atom, valencia, n_atom, individual_dbe
 866                    individual_dbe = individual_dbe + (n_atom * (valencia - 2))
 867                else:
 868                    continue
 869
 870        dbe = 1 + (0.5 * individual_dbe)
 871
 872        if self.ion_type == Labels.adduct_ion:
 873            dbe = dbe + 0.5
 874
 875        return dbe
 876
 877    def _calc_kmd(self, dict_base):
 878        """Calculate the Kendrick mass defect (KMD) of the molecular formula, based on the monoisotopic mass and the Kendrick mass.
 879
 880        Parameters
 881        ----------
 882        dict_base : dict
 883            The dictionary of the base formula, e.g. {'C':1, 'H':2}
 884
 885        Returns
 886        -------
 887        tuple
 888            The tuple of the KMD, Kendrick mass, and nominal Kendrick mass.
 889        """
 890        mass = 0
 891        for atom in dict_base.keys():
 892            mass = mass + Atoms.atomic_masses.get(atom) * dict_base.get(atom)
 893
 894        kendrick_mass = (int(mass) / mass) * self.mz_calc
 895
 896        nominal_km = int(kendrick_mass)
 897
 898        kmd = (nominal_km - kendrick_mass) * 100
 899
 900        # kmd = (nominal_km - km) * 1
 901        kmd = round(kmd, 0)
 902
 903        return kmd, kendrick_mass, nominal_km
 904
 905    def _cal_isotopologues(
 906        self, formula_dict, min_abundance, current_abundance, ms_dynamic_range
 907    ):
 908        """Calculate the isotopologues for a given molecular formula.
 909
 910        Parameters
 911        ----------
 912        formula_dict : dict
 913            The dictionary of the molecular formula. Example: {'C':10, 'H', 20, 'O', 2}
 914        min_abundance : float
 915            The minimum abundance.
 916        current_abundance : float
 917            The current monoisotopic abundance.
 918        ms_dynamic_range : float
 919            The dynamic range.
 920
 921
 922        Notes
 923        -----
 924        This is the primary function to look for isotopologues based on a monoisotopic molecular formula.
 925        It needs to be expanded to include the calculation of resolving power and plot the results.
 926        Use this function at runtime during the molecular identification algorithm only when a positive ID is observed to the monoisotopic ion.
 927        Use this function to simulate mass spectrum (needs resolving power calculation to be fully operational).
 928        It might break when adding non-conventional atoms (not yet tested).
 929        This function employs the IsoSpecPy library https://github.com/MatteoLacki/IsoSpec.
 930
 931
 932        """
 933
 934        # last update on 05-26-2020, Yuri E. Corilo
 935
 936        # updated it to reflect min possible mass peak abundance
 937        cut_off_to_IsoSpeccPy = 1 - (1 / ms_dynamic_range)
 938
 939        # print("cut_off_to_IsoSpeccPy", cut_off_to_IsoSpeccPy, current_abundance, min_abundance, ms_dynamic_range)
 940        # print(cut_off_to_IsoSpeccPy)
 941        atoms_labels = (
 942            atom
 943            for atom in formula_dict.keys()
 944            if atom != Labels.ion_type and atom != "H"
 945        )
 946
 947        atoms_count = []
 948        masses_list_tuples = []
 949        props_list_tuples = []
 950        all_atoms_list = []
 951
 952        for atom_label in atoms_labels:
 953            if Atoms.isotopes.get(atom_label)[1][0] is None:
 954                "This atom_label has no heavy isotope"
 955                atoms_count.append(formula_dict.get(atom_label))
 956                mass = Atoms.atomic_masses.get(atom_label)
 957                prop = Atoms.isotopic_abundance.get(atom_label)
 958                masses_list_tuples.append([mass])
 959                props_list_tuples.append([prop])
 960                all_atoms_list.append(atom_label)
 961
 962            else:
 963                isotopes_label_list = Atoms.isotopes.get(atom_label)[1]
 964
 965                if len(isotopes_label_list) > 1:
 966                    "This atom_label has two or more heavy isotope"
 967                    isotopos_labels = [i for i in isotopes_label_list]
 968                else:
 969                    "This atom_label only has one heavy isotope"
 970                    isotopos_labels = [isotopes_label_list[0]]
 971
 972                # all_atoms_list.extend(isotopos_labels)
 973                isotopos_labels = [atom_label] + isotopos_labels
 974
 975                all_atoms_list.extend(isotopos_labels)
 976
 977                masses = [
 978                    Atoms.atomic_masses.get(atom_label)
 979                    for atom_label in isotopos_labels
 980                ]
 981                props = [
 982                    Atoms.isotopic_abundance.get(atom_label)
 983                    for atom_label in isotopos_labels
 984                ]
 985
 986                atoms_count.append(formula_dict.get(atom_label))
 987                masses_list_tuples.append(masses)
 988                props_list_tuples.append(props)
 989        if legacy_isospec:
 990            iso = IsoSpecPy.IsoSpec(
 991                atoms_count,
 992                masses_list_tuples,
 993                props_list_tuples,
 994                cut_off_to_IsoSpeccPy,
 995            )
 996            conf = iso.getConfs()
 997            masses = conf[0]
 998            probs = exp(conf[1])
 999            molecular_formulas = conf[2]
1000            # print('conf', conf)
1001            # print('probs', conf[1])
1002        else:
1003            # This syntax in IsoSpecPy 2.2.2 yields the same information as the legacy approach
1004            iso = IsoSpecPy.IsoTotalProb(
1005                atomCounts=atoms_count,
1006                isotopeMasses=masses_list_tuples,
1007                isotopeProbabilities=props_list_tuples,
1008                prob_to_cover=cut_off_to_IsoSpeccPy,
1009                get_confs=True,
1010            )
1011            masses = list(iso.masses)
1012            probs = array(list(iso.probs))
1013            confs = list(iso.confs)
1014
1015            molecular_formulas = []
1016            for x in confs:
1017                tmplist = []
1018                for y in x:
1019                    tmplist.extend(list(y))
1020                molecular_formulas.append(tmplist)
1021
1022        new_formulas = []
1023
1024        for isotopologue_index in range(len(iso)):
1025            # skip_mono_isotopic
1026
1027            formula_list = molecular_formulas[isotopologue_index]
1028            new_formula_dict = dict(zip(all_atoms_list, formula_list))
1029            new_formula_dict[Labels.ion_type] = formula_dict.get(Labels.ion_type)
1030            if formula_dict.get("H"):
1031                new_formula_dict["H"] = formula_dict.get("H")
1032
1033            new_formulas.append({x: y for x, y in new_formula_dict.items() if y != 0})
1034
1035        # formula_dict in new_formulas check if monoisotopic is being returned
1036        if new_formulas:  # and formula_dict in new_formulas:
1037            # print(conf)
1038            # print(new_formulas)
1039            # print(atoms_count)
1040            # print(all_atoms_list)
1041            # print(masses_list_tuples)
1042            # print(props_list_tuples)
1043            # find where monoisotopic is
1044            index_mono = new_formulas.index(formula_dict)
1045            # calculate ratio iso/mono
1046            probs = list(probs / probs[index_mono])
1047
1048            # delete the monoisotopic
1049            del probs[index_mono]
1050            del new_formulas[index_mono]
1051
1052            # print('probs_exp', probs)
1053            for formulas, prob in zip(new_formulas, probs):
1054                theor_abundance = current_abundance * prob
1055                if theor_abundance > min_abundance:
1056                    # print(prob, theor_abundance, current_abundance)
1057                    yield (formulas, prob)
1058            # return zip(new_formulas, probs )
1059
1060        # else:
1061        #    return []