corems.molecular_id.factory.classification

  1__author__ = "Yuri E. Corilo"
  2__date__ = "Jan 31, 2020"
  3
  4from collections.abc import Mapping
  5
  6from matplotlib import pyplot as plt
  7from numpy import linspace
  8
  9from corems.encapsulation.constant import Atoms, Labels
 10
 11flatten_list = lambda l: [item for sublist in l for item in sublist]
 12
 13
 14class HeteroatomsClassification(Mapping):
 15    """Class for grouping mass spectrum data by heteroatom classes (Nn, Oo, Ss, NnOo, NnSs, etc..)
 16
 17    Parameters
 18    ----------
 19    mass_spectrum : MassSpectrum
 20        The mass spectrum object.
 21    choose_molecular_formula : bool, optional
 22        If True, the molecular formula with the highest abundance is chosen. If False, all molecular formulas are considered. Default is True.
 23
 24    Raises
 25    ------
 26    Exception
 27        If no molecular formula is associated with any mspeak objects.
 28
 29    Attributes
 30    ----------
 31    _ms_grouped_class : dict
 32        A dictionary of classes and a list of ms_peak objects.
 33    choose_mf : bool
 34        If True, the molecular formula with the highest abundance is chosen. If False, all molecular formulas are considered.
 35    total_peaks : int
 36        The total number of peaks.
 37    sum_abundance : float
 38        The sum of the abundance of all peaks.
 39    min_max_mz : tuple
 40        The minimum and maximum mz values.
 41    min_max_abundance : tuple
 42        The minimum and maximum abundance values.
 43    min_ppm_error : float
 44        The minimum ppm error.
 45    max_ppm_error : float
 46        The maximum ppm error.
 47    all_identified_atoms : list
 48        A list of all identified atoms.
 49
 50    Methods
 51    -------
 52    * __len__().
 53        Returns the number of classes.
 54    * __getitem__(classe)
 55        Returns the molecular formula list for specified class.
 56    * __iter__().
 57        Returns an iterator over the keys of the dictionary.
 58    * get_classes(threshold_perc=1, isotopologue=True).
 59        Returns a list of classes with abundance percentile above threshold.
 60    * molecular_formula_string(classe).
 61        Returns a list of molecular formula string for specified class.
 62    * molecular_formula(classe).
 63        Returns a list of molecular formula for specified class.
 64    * carbon_number(classe).
 65        Returns a list of carbon number for specified class.
 66    * atom_count(atom, classe).
 67        Returns a list of atom count for specified class.
 68    * dbe(classe).
 69        Returns a list of DBE for specified class.
 70    * atoms_ratio(classe, numerator, denominator).
 71        Returns a list of atoms ratio for specified class.
 72    * mz_exp(classe).
 73        Returns a list of experimental mz for specified class.
 74    * abundance(classe).
 75        Returns a list of abundance for specified class.
 76    * mz_error(classe).
 77        Returns a list of mz error for specified class.
 78    * mz_calc(classe).
 79        Returns a list of calculated mz for specified class.
 80    * peaks_count_percentile(classe).
 81        Returns the peaks count percentile of a specific class.
 82    * abundance_count_percentile(classe).
 83        Returns the abundance percentile of a specific class.
 84    * mz_exp_assigned().
 85        Returns a list of experimental mz for all assigned classes.
 86    * abundance_assigned().
 87        Returns a list of abundance for all classes.
 88    * mz_exp_all().
 89        Returns a list of mz for all classes.
 90
 91    """
 92
 93    # Group mass spectrum data by heteroatom classes (Nn, Oo, Ss, NnOo, NnSs, etc..)
 94
 95    #   class obj behaves as a dictionary of classes and return a list of ms_peak obj
 96
 97    def __init__(self, mass_spectrum, choose_molecular_formula=True):
 98        def sort_atoms_method(atom):
 99            """Sort atoms by order of appearance in the Atoms class"""
100            return [Atoms.atoms_order.index(atom)]
101
102        self._ms_grouped_class = dict()
103
104        self.choose_mf = choose_molecular_formula
105
106        # mapping for ms peaks without any molecular formula associated
107        self._ms_grouped_class[Labels.unassigned] = list()
108
109        self.total_peaks = 0
110
111        self.sum_abundance = 0
112
113        self.min_max_mz = (mass_spectrum.min_mz_exp, mass_spectrum.max_mz_exp)
114
115        self.min_max_abundance = (
116            mass_spectrum.min_abundance,
117            mass_spectrum.max_abundance,
118        )
119
120        self.min_ppm_error = mass_spectrum.molecular_search_settings.min_ppm_error
121
122        self.max_ppm_error = mass_spectrum.molecular_search_settings.max_ppm_error
123
124        check_assign = False
125
126        all_used_atoms = set()
127
128        for ms_peak in mass_spectrum:
129            self.total_peaks += 1
130
131            self.sum_abundance += ms_peak.abundance
132
133            if not ms_peak.is_assigned:
134                self._ms_grouped_class.get(Labels.unassigned).append(ms_peak)
135
136            else:
137                check_assign = True
138
139                if choose_molecular_formula:
140                    mf = ms_peak.best_molecular_formula_candidate
141
142                    classes = [mf.class_label]
143
144                    for atom in mf.atoms:
145                        all_used_atoms.add(atom)
146
147                else:
148                    classes = []
149
150                    for mf in ms_peak:
151                        classes.append(mf.class_label)
152
153                        for atom in mf.atoms:
154                            all_used_atoms.add(atom)
155
156                for classe in classes:
157                    if classe in self._ms_grouped_class.keys():
158                        self._ms_grouped_class.get(classe).append(ms_peak)
159
160                    else:
161                        self._ms_grouped_class[classe] = [ms_peak]
162
163        self.all_identified_atoms = sorted(all_used_atoms, key=sort_atoms_method)
164
165        if not check_assign:
166            raise Exception("No molecular formula associated with any mspeak objects")
167
168    def __len__(self):
169        """Return the number of classes"""
170        return len(self._ms_grouped_class)
171
172    def __getitem__(self, classe):
173        """Return the molecular formula list for specified class"""
174        return self._ms_grouped_class.get(classe)
175
176    def __iter__(self):
177        """Return an iterator over the keys of the dictionary."""
178        return iter(self._ms_grouped_class)
179
180    def get_classes(self, threshold_perc=1, isotopologue=True):
181        """Return a list of classes with abundance percentile above threshold"""
182        classes = list()
183        for classe in self.keys():
184            if classe != Labels.unassigned:
185                if self.abundance_count_percentile(classe) > threshold_perc:
186                    if classe != Labels.unassigned:
187                        # access first molecular formula inside the first ms peak and check isotopologue
188                        if not isotopologue and self.get(classe)[0][0].is_isotopologue:
189                            continue
190
191                    classes.append(classe)
192        # TODO sort classes chemically here too
193        return classes
194
195    def molecular_formula_string(
196        self,
197        classe,
198    ):
199        """Return a list of molecular formula string for specified class"""
200        if self.choose_mf:
201            return [mspeak.best_molecular_formula_candidate for mspeak in self[classe]]
202        else:
203            return [
204                mf
205                for mspeak in self[classe]
206                for mf in mspeak
207                if mf.class_label == classe
208            ]
209
210    def molecular_formula(
211        self,
212        classe,
213    ):
214        """Return a list of molecular formula for specified class"""
215        if self.choose_mf:
216            return [mspeak.best_molecular_formula_candidate for mspeak in self[classe]]
217        else:
218            return [
219                mf
220                for mspeak in self[classe]
221                for mf in mspeak
222                if mf.class_label == classe
223            ]
224
225    def carbon_number(self, classe):
226        """Return a list of carbon number for specified class"""
227        if self.choose_mf:
228            return [
229                mspeak.best_molecular_formula_candidate.get("C")
230                for mspeak in self[classe]
231            ]
232        else:
233            return [
234                mf.get("C")
235                for mspeak in self[classe]
236                for mf in mspeak
237                if mf.class_label == classe
238            ]
239
240    def atom_count(self, atom, classe):
241        """Return a list of atom count for specified class"""
242
243        if self.choose_mf:
244            return [
245                mspeak.best_molecular_formula_candidate.get(atom)
246                for mspeak in self[classe]
247            ]
248        else:
249            return [
250                mf.get(atom)
251                for mspeak in self[classe]
252                for mf in mspeak
253                if mf.class_label == classe
254            ]
255
256    def dbe(self, classe):
257        """Return a list of DBE for specified class"""
258        if self.choose_mf:
259            return [
260                mspeak.best_molecular_formula_candidate.dbe for mspeak in self[classe]
261            ]
262        else:
263            return [
264                mf.dbe
265                for mspeak in self[classe]
266                for mf in mspeak
267                if mf.class_label == classe
268            ]
269
270    def atoms_ratio(self, classe, numerator, denominator):
271        """Return a list of atoms ratio for specified class"""
272        return [
273            mf.get(numerator) / mf.get(denominator)
274            for mf in self.molecular_formula(classe)
275        ]
276
277    def mz_exp(self, classe):
278        """Return a list of experimental mz for specified class"""
279        if self.choose_mf or classe == Labels.unassigned:
280            return [mspeak.mz_exp for mspeak in self[classe]]
281
282        else:
283            return [
284                mspeak.mz_exp
285                for mspeak in self[classe]
286                for mf in mspeak
287                if mf.class_label == classe
288            ]
289
290    def abundance(self, classe):
291        """Return a list of abundance for specified class"""
292        if self.choose_mf or classe == Labels.unassigned:
293            return [mspeak.abundance for mspeak in self[classe]]
294
295        else:
296            return [
297                mspeak.abundance
298                for mspeak in self[classe]
299                for mf in mspeak
300                if mf.class_label == classe
301            ]
302
303    def mz_error(self, classe):
304        """Return a list of mz error for specified class"""
305        if classe != Labels.unassigned:
306            if self.choose_mf:
307                return [
308                    mspeak.best_molecular_formula_candidate.mz_error
309                    for mspeak in self[classe]
310                ]
311
312            else:
313                return [
314                    mf.mz_error
315                    for mspeak in self[classe]
316                    for mf in mspeak
317                    if mf.class_label == classe
318                ]
319
320    def mz_calc(self, classe):
321        """Return a list of calculated mz for specified class"""
322        if self.choose_mf:
323            return [
324                mspeak.best_molecular_formula_candidate.mz_calc
325                for mspeak in self[classe]
326            ]
327
328        else:
329            return [
330                mf.mz_calc
331                for mspeak in self[classe]
332                for mf in mspeak
333                if mf.class_label == classe
334            ]
335
336    def peaks_count_percentile(self, classe):
337        """Return the peaks count percentile of a specific class"""
338        return (len(self[classe]) / self.total_peaks) * 100
339
340    def abundance_count_percentile(self, classe):
341        """Return the abundance percentile of a specific class"""
342        return (
343            sum([mspeak.abundance for mspeak in self[classe]]) / self.sum_abundance
344        ) * 100
345
346    def mz_exp_assigned(self):
347        """Return a list of experimental mz for all assigned classes"""
348        classes = self.keys()
349
350        return [
351            mspeak.mz_exp
352            for classe in classes
353            for mspeak in self[classe]
354            if classe != Labels.unassigned
355        ]
356
357    def abundance_assigned(self):
358        """Return a list of abundance for all classes"""
359        classes = self.keys()
360
361        return [
362            mspeak.abundance
363            for classe in classes
364            for mspeak in self[classe]
365            if classe != Labels.unassigned
366        ]
367
368    def mz_exp_all(self):
369        """Return a list of mz for all classes"""
370        classes = self.keys()
371
372        return flatten_list(
373            [self.mz_exp(classe) for classe in classes if classe != Labels.unassigned]
374        )
375
376    def mz_error_all(self):
377        """Return a list of mz error for all classes"""
378        classes = self.keys()
379
380        return flatten_list(
381            [self.mz_error(classe) for classe in classes if classe != Labels.unassigned]
382        )
383
384    def carbon_number_all(self):
385        """Return a list of carbon number for all classes"""
386        classes = self.keys()
387
388        return flatten_list(
389            [
390                self.carbon_number(classe)
391                for classe in classes
392                if classe != Labels.unassigned
393            ]
394        )
395
396    def dbe_all(self):
397        """Return a list of DBE for all classes"""
398        classes = self.keys()
399
400        return flatten_list(
401            [self.dbe(classe) for classe in classes if classe != Labels.unassigned]
402        )
403
404    def atoms_ratio_all(self, numerator, denominator):
405        """Return a list of atoms ratio for all classes"""
406        classes = self.keys()
407
408        return flatten_list(
409            [
410                self.atoms_ratio(classe, numerator, denominator)
411                for classe in classes
412                if classe != Labels.unassigned
413            ]
414        )
415
416    def to_dataframe(
417        self,
418        include_isotopologue=False,
419        abundance_perc_threshold=5,
420        include_unassigned=False,
421    ):
422        """Return a pandas dataframe with all the data from the class
423
424        Parameters
425        ----------
426        include_isotopologue : bool, optional
427            Include isotopologues, by default False
428        abundance_perc_threshold : int, optional
429            Abundance percentile threshold, by default 5
430        include_unassigned : bool, optional
431            Include unassigned peaks, by default False
432
433        Returns
434        -------
435        DataFrame
436            Pandas dataframe with all the data from the class
437        """
438        from pandas import DataFrame
439
440        columns_labels = [
441            "mz",
442            "calibrated_mz",
443            "calculated_m_z",
444            "abundance",
445            "resolving_power",
446            "sn",
447            "ion_charge",
448            "mass_error",
449            "DBE",
450            "class",
451            "HC",
452            "OC",
453            "ion_type",
454            "is_isotopologue",
455            "class_abundance",
456            "class_count",
457        ]
458
459        dict_data_list = []
460
461        for classe, list_mspeaks in self.items():
462            percent_abundance = self.abundance_count_percentile(classe)
463
464            # ignores low abundant classes
465            if abundance_perc_threshold < abundance_perc_threshold:
466                continue
467
468            peaks_count_percentile = self.peaks_count_percentile(classe)
469
470            for ms_peak in list_mspeaks:
471                if ms_peak.is_assigned:
472                    for m_formula in ms_peak:
473                        # ignores isotopologues
474                        if not include_isotopologue and m_formula.is_isotopologue:
475                            continue
476
477                        formula_dict = m_formula.to_dict()
478
479                        dict_result = {
480                            "mz": ms_peak._mz_exp,
481                            "calibrated_mz": ms_peak.mz_exp,
482                            "calculated_mz": m_formula.mz_calc,
483                            "abundance": ms_peak.abundance,
484                            "resolving_power": ms_peak.resolving_power,
485                            "sn": ms_peak.signal_to_noise,
486                            # Assigned: formula charge (may be multi-charge), not peak polarity
487                            "ion_charge": m_formula.ion_charge,
488                            "mass_error": m_formula.mz_error,
489                            "DBE": m_formula.dbe,
490                            "class": classe,
491                            "HC": m_formula.H_C,
492                            "OC": m_formula.O_C,
493                            "ion_type": str(m_formula.ion_type.lower().encode("utf-8")),
494                            "is_isotopologue": int(m_formula.is_isotopologue),
495                            "class_abundance": percent_abundance,
496                            "class_count": peaks_count_percentile,
497                        }
498
499                        for atom in formula_dict.keys():
500                            dict_result[atom] = formula_dict.get(atom)
501
502                    dict_data_list.append(dict_result)
503
504                else:
505                    if not include_unassigned:
506                        continue
507
508                    dict_result = {
509                        "mz": ms_peak._mz_exp,
510                        "calibrated_mz": ms_peak.mz_exp,
511                        "abundance": ms_peak.abundance,
512                        "resolving_power": ms_peak.resolving_power,
513                        "sn": ms_peak.signal_to_noise,
514                        # Peak polarity (±1); no formula charge for unassigned peaks
515                        "ion_charge": ms_peak.polarity,
516                        "class": classe,
517                        "class_abundance": percent_abundance,
518                        "class_count": percent_abundance,
519                    }
520
521                    dict_data_list.append(dict_result)
522
523        columns = columns_labels + self.all_identified_atoms
524
525        return DataFrame(dict_data_list, columns=columns)
526
527    def plot_ms_assigned_unassigned(self, assigned_color="b", unassigned_color="r"):
528        """Plot stick mass spectrum of all classes
529
530        Parameters
531        ----------
532        assigned_color : str, optional
533            Matplotlib color for the assigned peaks, by default "b"
534        unassigned_color : str, optional
535            Matplotlib color for the unassigned peaks, by default "r"
536
537        Returns
538        -------
539        ax : matplotlib.axes
540            Matplotlib axes object
541        """
542        mz_assigned = self.mz_exp_assigned()
543        abundance_assigned = self.abundance_assigned()
544
545        mz_not_assigned = self.mz_exp(Labels.unassigned)
546        abundance_not_assigned = self.abundance(Labels.unassigned)
547
548        ax = plt.gca()
549
550        for plot_obj in ax.stem(
551            mz_assigned,
552            abundance_assigned,
553            linefmt="-",
554            markerfmt=" ",
555            label="Assigned",
556        ):
557            plt.setp(plot_obj, "color", assigned_color, "linewidth", 2)
558
559        for plot_obj in ax.stem(
560            mz_not_assigned,
561            abundance_not_assigned,
562            linefmt="-",
563            markerfmt=" ",
564            label="Unassigned",
565        ):
566            plt.setp(plot_obj, "color", unassigned_color, "linewidth", 2)
567
568        ax.set_xlabel("$\t{m/z}$", fontsize=12)
569        ax.set_ylabel("Abundance", fontsize=12)
570        ax.tick_params(axis="both", which="major", labelsize=12)
571
572        ax.axes.spines["top"].set_visible(False)
573        ax.axes.spines["right"].set_visible(False)
574
575        ax.get_yaxis().set_visible(False)
576        ax.spines["left"].set_visible(False)
577        plt.legend()
578
579        return ax
580
581    def plot_mz_error(self, color="g"):
582        """Plot mz error scatter plot of all classes
583
584        Parameters
585        ----------
586        color : str, optional
587            Matplotlib color, by default "g"
588
589        Returns
590        -------
591        ax : matplotlib.axes
592            Matplotlib axes object
593        """
594        ax = plt.gca()
595
596        mz_assigned = self.mz_exp_all()
597        mz_error = self.mz_error_all()
598
599        ax.scatter(mz_assigned, mz_error, c=color)
600
601        ax.set_xlabel("$\t{m/z}$", fontsize=12)
602        ax.set_ylabel("Error (ppm)", fontsize=12)
603        ax.tick_params(axis="both", which="major", labelsize=12)
604
605        ax.axes.spines["top"].set_visible(True)
606        ax.axes.spines["right"].set_visible(True)
607
608        ax.get_yaxis().set_visible(True)
609        ax.spines["left"].set_visible(True)
610
611        ax.set_xlim(self.min_max_mz)
612        ax.set_ylim(self.min_ppm_error, self.max_ppm_error)
613
614        return ax
615
616    def plot_mz_error_class(self, classe, color="g"):
617        """Plot mz error scatter plot of a specific class
618
619        Parameters
620        ----------
621        classe : str
622            Class name
623        color : str, optional
624            Matplotlib color, by default "g"
625
626        Returns
627        -------
628        ax : matplotlib.axes
629            Matplotlib axes object
630
631        """
632        if classe != Labels.unassigned:
633            ax = plt.gca()
634
635            abun_perc = self.abundance_count_percentile(classe)
636            mz_assigned = self.mz_exp(classe)
637            mz_error = self.mz_error(classe)
638
639            ax.scatter(mz_assigned, mz_error, c=color)
640
641            title = "%s, %.2f %%" % (classe, abun_perc)
642            ax.set_title(title)
643            ax.set_xlabel("$\t{m/z}$", fontsize=12)
644            ax.set_ylabel("Error (ppm)", fontsize=12)
645            ax.tick_params(axis="both", which="major", labelsize=12)
646
647            ax.axes.spines["top"].set_visible(True)
648            ax.axes.spines["right"].set_visible(True)
649
650            ax.get_yaxis().set_visible(True)
651            ax.spines["left"].set_visible(True)
652
653            ax.set_xlim(self.min_max_mz)
654            ax.set_ylim(self.min_ppm_error, self.max_ppm_error)
655
656            return ax
657
658    def plot_ms_class(self, classe, color="g"):
659        """Plot stick mass spectrum of a specific class
660
661        Parameters
662        ----------
663        classe : str
664            Class name
665        color : str, optional
666            Matplotlib color, by default "g"
667
668        Returns
669        -------
670        ax : matplotlib.axes
671            Matplotlib axes object
672
673        """
674        if classe != Labels.unassigned:
675            ax = plt.gca()
676
677            abun_perc = self.abundance_count_percentile(classe)
678            mz_assigned = self.mz_exp(classe)
679            abundance_assigned = self.abundance(classe)
680
681            for plot_obj in ax.stem(
682                mz_assigned, abundance_assigned, linefmt="-", markerfmt=" "
683            ):
684                plt.setp(plot_obj, "color", color, "linewidth", 2)
685
686            title = "%s, %.2f %%" % (classe, abun_perc)
687            ax.set_title(title)
688            ax.set_xlabel("$\t{m/z}$", fontsize=12)
689            ax.set_ylabel("Abundance", fontsize=12)
690            ax.tick_params(axis="both", which="major", labelsize=12)
691
692            ax.axes.spines["top"].set_visible(False)
693            ax.axes.spines["right"].set_visible(False)
694
695            ax.get_yaxis().set_visible(False)
696            ax.spines["left"].set_visible(False)
697
698            ax.set_xlim(self.min_max_mz)
699            ax.set_ylim(self.min_max_abundance)
700
701            return ax
702
703    def plot_van_krevelen(
704        self,
705        classe=None,
706        max_hc=2.5,
707        max_oc=2,
708        ticks_number=5,
709        color="viridis",
710        alpha=0.5,
711        log_abundance=False,
712    ):
713        """Plot Van Krevelen Diagram for a single class or all assigned classes
714
715        Parameters
716        ----------
717        classe : str, optional
718            Class name or None to plot all assigned classes, by default None
719        max_hc : float, optional
720            Max H/C ratio, by default 2.5
721        max_oc : float, optional
722            Max O/C ratio, by default 2
723        ticks_number : int, optional
724            Number of ticks, by default 5
725        color : str, optional
726            Matplotlib color/colormap, by default "viridis"
727        alpha : float, optional
728            Transparency of points, by default 0.5
729        log_abundance : bool, optional
730            If True, use log10 scale for abundance values, by default False
731
732        Returns
733        -------
734        ax : matplotlib.axes
735            Matplotlib axes object
736        abun_perc : float or None
737            Class percentile of the relative abundance (if classe specified)
738        """
739        import numpy as np
740
741        ax = plt.gca()
742
743        if classe is not None and classe != Labels.unassigned:
744            # Single class plot
745            abun_perc = self.abundance_count_percentile(classe)
746            hc = self.atoms_ratio(classe, "H", "C")
747            oc = self.atoms_ratio(classe, "O", "C")
748            abundance = self.abundance(classe)
749
750            if log_abundance:
751                abundance = [np.log10(a + 1e-10) for a in abundance]
752                colorbar_label = "log\u2081\u2080(Abundance)"
753            else:
754                colorbar_label = "Abundance"
755
756            # Sort by abundance so higher values are plotted on top
757            sorted_indices = sorted(range(len(abundance)), key=lambda i: abundance[i])
758            hc = [hc[i] for i in sorted_indices]
759            oc = [oc[i] for i in sorted_indices]
760            abundance = [abundance[i] for i in sorted_indices]
761
762            scatter = ax.scatter(oc, hc, c=abundance, alpha=alpha, cmap=color)
763            plt.colorbar(scatter, label=colorbar_label)
764
765            title = "%s, %.2f %%" % (classe, abun_perc)
766            ax.set_title(title)
767
768            return_val = ax, abun_perc
769        else:
770            # All assigned classes plot
771            hc = self.atoms_ratio_all("H", "C")
772            oc = self.atoms_ratio_all("O", "C")
773            abundance = self.abundance_assigned()
774
775            if log_abundance:
776                abundance = [np.log10(a + 1e-10) for a in abundance]
777                colorbar_label = "log\u2081\u2080(Abundance)"
778            else:
779                colorbar_label = "Abundance"
780
781            sorted_indices = sorted(range(len(abundance)), key=lambda i: abundance[i])
782            hc = [hc[i] for i in sorted_indices]
783            oc = [oc[i] for i in sorted_indices]
784            abundance = [abundance[i] for i in sorted_indices]
785
786            scatter = ax.scatter(oc, hc, c=abundance, alpha=alpha, cmap=color)
787            plt.colorbar(scatter, label=colorbar_label)
788
789            ax.set_title("Van Krevelen Diagram - All Assigned Classes")
790
791            return_val = ax
792
793        ax.set_xlabel("O/C", fontsize=16)
794        ax.set_ylabel("H/C", fontsize=16)
795        ax.tick_params(axis="both", which="major", labelsize=18)
796        ax.set_xticks(linspace(0, max_oc, ticks_number, endpoint=True))
797        ax.set_yticks(linspace(0, max_hc, ticks_number, endpoint=True))
798        ax.grid(alpha=0.3, linestyle="--")
799
800        return return_val
801
802    def plot_dbe_vs_carbon_number(
803        self,
804        classe=None,
805        max_c=50,
806        max_dbe=40,
807        dbe_incr=5,
808        c_incr=10,
809        color="viridis",
810        alpha=0.5,
811        log_abundance=False,
812    ):
813        """Plot DBE vs Carbon Number for a single class or all assigned classes
814
815        Parameters
816        ----------
817        classe : str, optional
818            Class name or None to plot all assigned classes, by default None
819        max_c : int, optional
820            Max Carbon Number, by default 50
821        max_dbe : int, optional
822            Max DBE, by default 40
823        dbe_incr : int, optional
824            DBE increment, by default 5
825        c_incr : int, optional
826            Carbon Number increment, by default 10
827        color : str, optional
828            Matplotlib color/colormap, by default "viridis"
829        alpha : float, optional
830            Transparency of points, by default 0.5
831        log_abundance : bool, optional
832            If True, use log10 scale for abundance values, by default False
833
834        Returns
835        -------
836        ax : matplotlib.axes
837            Matplotlib axes object
838        abun_perc : float or None
839            Class percentile of the relative abundance (if classe specified)
840        """
841        import numpy as np
842
843        ax = plt.gca()
844
845        if classe is not None and classe != Labels.unassigned:
846            # Single class plot
847            abun_perc = self.abundance_count_percentile(classe)
848            carbon_number = self.carbon_number(classe)
849            dbe = self.dbe(classe)
850            abundance = self.abundance(classe)
851
852            if log_abundance:
853                abundance = [np.log10(a + 1e-10) for a in abundance]
854                colorbar_label = "log\u2081\u2080(Abundance)"
855            else:
856                colorbar_label = "Abundance"
857
858            sorted_indices = sorted(range(len(abundance)), key=lambda i: abundance[i])
859            carbon_number = [carbon_number[i] for i in sorted_indices]
860            dbe = [dbe[i] for i in sorted_indices]
861            abundance = [abundance[i] for i in sorted_indices]
862
863            scatter = ax.scatter(
864                carbon_number, dbe, c=abundance, alpha=alpha, cmap=color
865            )
866            plt.colorbar(scatter, label=colorbar_label)
867
868            title = "%s, %.2f %%" % (classe, abun_perc)
869            ax.set_title(title)
870
871            return_val = ax, abun_perc
872        else:
873            # All assigned classes plot
874            carbon_number = self.carbon_number_all()
875            dbe = self.dbe_all()
876            abundance = self.abundance_assigned()
877
878            if log_abundance:
879                abundance = [np.log10(a + 1e-10) for a in abundance]
880                colorbar_label = "log\u2081\u2080(Abundance)"
881            else:
882                colorbar_label = "Abundance"
883
884            sorted_indices = sorted(range(len(abundance)), key=lambda i: abundance[i])
885            carbon_number = [carbon_number[i] for i in sorted_indices]
886            dbe = [dbe[i] for i in sorted_indices]
887            abundance = [abundance[i] for i in sorted_indices]
888
889            scatter = ax.scatter(
890                carbon_number, dbe, c=abundance, alpha=alpha, cmap=color
891            )
892            plt.colorbar(scatter, label=colorbar_label)
893
894            ax.set_title("DBE vs Carbon Number - All Assigned Classes")
895
896            return_val = ax
897
898        ax.set_xlabel("Carbon number", fontsize=16)
899        ax.set_ylabel("DBE", fontsize=16)
900        ax.tick_params(axis="both", which="major", labelsize=18)
901        ax.set_xticks(range(0, max_c, c_incr))
902        ax.set_yticks(range(0, max_dbe, dbe_incr))
903        ax.grid(alpha=0.3, linestyle="--")
904
905        return return_val