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                            "ion_charge": ms_peak.ion_charge,
487                            "mass_error": m_formula.mz_error,
488                            "DBE": m_formula.dbe,
489                            "class": classe,
490                            "HC": m_formula.H_C,
491                            "OC": m_formula.O_C,
492                            "ion_type": str(m_formula.ion_type.lower().encode("utf-8")),
493                            "is_isotopologue": int(m_formula.is_isotopologue),
494                            "class_abundance": percent_abundance,
495                            "class_count": peaks_count_percentile,
496                        }
497
498                        for atom in formula_dict.keys():
499                            dict_result[atom] = formula_dict.get(atom)
500
501                    dict_data_list.append(dict_result)
502
503                else:
504                    if not include_unassigned:
505                        continue
506
507                    dict_result = {
508                        "mz": ms_peak._mz_exp,
509                        "calibrated_mz": ms_peak.mz_exp,
510                        "abundance": ms_peak.abundance,
511                        "resolving_power": ms_peak.resolving_power,
512                        "sn": ms_peak.signal_to_noise,
513                        "ion_charge": ms_peak.ion_charge,
514                        "class": classe,
515                        "class_abundance": percent_abundance,
516                        "class_count": percent_abundance,
517                    }
518
519                    dict_data_list.append(dict_result)
520
521        columns = columns_labels + self.all_identified_atoms
522
523        return DataFrame(dict_data_list, columns=columns)
524
525    def plot_ms_assigned_unassigned(self, assigned_color="b", unassigned_color="r"):
526        """Plot stick mass spectrum of all classes
527
528        Parameters
529        ----------
530        assigned_color : str, optional
531            Matplotlib color for the assigned peaks, by default "b"
532        unassigned_color : str, optional
533            Matplotlib color for the unassigned peaks, by default "r"
534
535        Returns
536        -------
537        ax : matplotlib.axes
538            Matplotlib axes object
539        """
540        mz_assigned = self.mz_exp_assigned()
541        abundance_assigned = self.abundance_assigned()
542
543        mz_not_assigned = self.mz_exp(Labels.unassigned)
544        abundance_not_assigned = self.abundance(Labels.unassigned)
545
546        ax = plt.gca()
547
548        for plot_obj in ax.stem(
549            mz_assigned,
550            abundance_assigned,
551            linefmt="-",
552            markerfmt=" ",
553            label="Assigned",
554        ):
555            plt.setp(plot_obj, "color", assigned_color, "linewidth", 2)
556
557        for plot_obj in ax.stem(
558            mz_not_assigned,
559            abundance_not_assigned,
560            linefmt="-",
561            markerfmt=" ",
562            label="Unassigned",
563        ):
564            plt.setp(plot_obj, "color", unassigned_color, "linewidth", 2)
565
566        ax.set_xlabel("$\t{m/z}$", fontsize=12)
567        ax.set_ylabel("Abundance", fontsize=12)
568        ax.tick_params(axis="both", which="major", labelsize=12)
569
570        ax.axes.spines["top"].set_visible(False)
571        ax.axes.spines["right"].set_visible(False)
572
573        ax.get_yaxis().set_visible(False)
574        ax.spines["left"].set_visible(False)
575        plt.legend()
576
577        return ax
578
579    def plot_mz_error(self, color="g"):
580        """Plot mz error scatter plot of all classes
581
582        Parameters
583        ----------
584        color : str, optional
585            Matplotlib color, by default "g"
586
587        Returns
588        -------
589        ax : matplotlib.axes
590            Matplotlib axes object
591        """
592        ax = plt.gca()
593
594        mz_assigned = self.mz_exp_all()
595        mz_error = self.mz_error_all()
596
597        ax.scatter(mz_assigned, mz_error, c=color)
598
599        ax.set_xlabel("$\t{m/z}$", fontsize=12)
600        ax.set_ylabel("Error (ppm)", fontsize=12)
601        ax.tick_params(axis="both", which="major", labelsize=12)
602
603        ax.axes.spines["top"].set_visible(True)
604        ax.axes.spines["right"].set_visible(True)
605
606        ax.get_yaxis().set_visible(True)
607        ax.spines["left"].set_visible(True)
608
609        ax.set_xlim(self.min_max_mz)
610        ax.set_ylim(self.min_ppm_error, self.max_ppm_error)
611
612        return ax
613
614    def plot_mz_error_class(self, classe, color="g"):
615        """Plot mz error scatter plot of a specific class
616
617        Parameters
618        ----------
619        classe : str
620            Class name
621        color : str, optional
622            Matplotlib color, by default "g"
623
624        Returns
625        -------
626        ax : matplotlib.axes
627            Matplotlib axes object
628
629        """
630        if classe != Labels.unassigned:
631            ax = plt.gca()
632
633            abun_perc = self.abundance_count_percentile(classe)
634            mz_assigned = self.mz_exp(classe)
635            mz_error = self.mz_error(classe)
636
637            ax.scatter(mz_assigned, mz_error, c=color)
638
639            title = "%s, %.2f %%" % (classe, abun_perc)
640            ax.set_title(title)
641            ax.set_xlabel("$\t{m/z}$", fontsize=12)
642            ax.set_ylabel("Error (ppm)", fontsize=12)
643            ax.tick_params(axis="both", which="major", labelsize=12)
644
645            ax.axes.spines["top"].set_visible(True)
646            ax.axes.spines["right"].set_visible(True)
647
648            ax.get_yaxis().set_visible(True)
649            ax.spines["left"].set_visible(True)
650
651            ax.set_xlim(self.min_max_mz)
652            ax.set_ylim(self.min_ppm_error, self.max_ppm_error)
653
654            return ax
655
656    def plot_ms_class(self, classe, color="g"):
657        """Plot stick mass spectrum of a specific class
658
659        Parameters
660        ----------
661        classe : str
662            Class name
663        color : str, optional
664            Matplotlib color, by default "g"
665
666        Returns
667        -------
668        ax : matplotlib.axes
669            Matplotlib axes object
670
671        """
672        if classe != Labels.unassigned:
673            ax = plt.gca()
674
675            abun_perc = self.abundance_count_percentile(classe)
676            mz_assigned = self.mz_exp(classe)
677            abundance_assigned = self.abundance(classe)
678
679            for plot_obj in ax.stem(
680                mz_assigned, abundance_assigned, linefmt="-", markerfmt=" "
681            ):
682                plt.setp(plot_obj, "color", color, "linewidth", 2)
683
684            title = "%s, %.2f %%" % (classe, abun_perc)
685            ax.set_title(title)
686            ax.set_xlabel("$\t{m/z}$", fontsize=12)
687            ax.set_ylabel("Abundance", fontsize=12)
688            ax.tick_params(axis="both", which="major", labelsize=12)
689
690            ax.axes.spines["top"].set_visible(False)
691            ax.axes.spines["right"].set_visible(False)
692
693            ax.get_yaxis().set_visible(False)
694            ax.spines["left"].set_visible(False)
695
696            ax.set_xlim(self.min_max_mz)
697            ax.set_ylim(self.min_max_abundance)
698
699            return ax
700
701    def plot_van_krevelen(
702        self, classe=None, max_hc=2.5, max_oc=2, ticks_number=5, color="viridis",
703        alpha=0.5, log_abundance=False
704    ):
705        """Plot Van Krevelen Diagram for a single class or all assigned classes
706
707        Parameters
708        ----------
709        classe : str, optional
710            Class name or None to plot all assigned classes, by default None
711        max_hc : float, optional
712            Max H/C ratio, by default 2.5
713        max_oc : float, optional
714            Max O/C ratio, by default 2
715        ticks_number : int, optional
716            Number of ticks, by default 5
717        color : str, optional
718            Matplotlib color/colormap, by default "viridis"
719        alpha : float, optional
720            Transparency of points, by default 0.5
721        log_abundance : bool, optional
722            If True, use log10 scale for abundance values, by default False
723
724        Returns
725        -------
726        ax : matplotlib.axes
727            Matplotlib axes object
728        abun_perc : float or None
729            Class percentile of the relative abundance (if classe specified)
730        """
731        import numpy as np
732        ax = plt.gca()
733
734        if classe is not None and classe != Labels.unassigned:
735            # Single class plot
736            abun_perc = self.abundance_count_percentile(classe)
737            hc = self.atoms_ratio(classe, "H", "C")
738            oc = self.atoms_ratio(classe, "O", "C")
739            abundance = self.abundance(classe)
740
741            if log_abundance:
742                abundance = [np.log10(a + 1e-10) for a in abundance]
743                colorbar_label = 'log\u2081\u2080(Abundance)'
744            else:
745                colorbar_label = 'Abundance'
746
747            # Sort by abundance so higher values are plotted on top
748            sorted_indices = sorted(range(len(abundance)), key=lambda i: abundance[i])
749            hc = [hc[i] for i in sorted_indices]
750            oc = [oc[i] for i in sorted_indices]
751            abundance = [abundance[i] for i in sorted_indices]
752
753            scatter = ax.scatter(oc, hc, c=abundance, alpha=alpha, cmap=color)
754            plt.colorbar(scatter, label=colorbar_label)
755
756            title = "%s, %.2f %%" % (classe, abun_perc)
757            ax.set_title(title)
758
759            return_val = ax, abun_perc
760        else:
761            # All assigned classes plot
762            hc = self.atoms_ratio_all("H", "C")
763            oc = self.atoms_ratio_all("O", "C")
764            abundance = self.abundance_assigned()
765
766            if log_abundance:
767                abundance = [np.log10(a + 1e-10) for a in abundance]
768                colorbar_label = 'log\u2081\u2080(Abundance)'
769            else:
770                colorbar_label = 'Abundance'
771
772            sorted_indices = sorted(range(len(abundance)), key=lambda i: abundance[i])
773            hc = [hc[i] for i in sorted_indices]
774            oc = [oc[i] for i in sorted_indices]
775            abundance = [abundance[i] for i in sorted_indices]
776
777            scatter = ax.scatter(oc, hc, c=abundance, alpha=alpha, cmap=color)
778            plt.colorbar(scatter, label=colorbar_label)
779
780            ax.set_title("Van Krevelen Diagram - All Assigned Classes")
781
782            return_val = ax
783
784        ax.set_xlabel("O/C", fontsize=16)
785        ax.set_ylabel("H/C", fontsize=16)
786        ax.tick_params(axis="both", which="major", labelsize=18)
787        ax.set_xticks(linspace(0, max_oc, ticks_number, endpoint=True))
788        ax.set_yticks(linspace(0, max_hc, ticks_number, endpoint=True))
789        ax.grid(alpha=0.3, linestyle='--')
790
791        return return_val
792
793    def plot_dbe_vs_carbon_number(
794        self, classe=None, max_c=50, max_dbe=40, dbe_incr=5, c_incr=10, color="viridis",
795        alpha=0.5, log_abundance=False
796    ):
797        """Plot DBE vs Carbon Number for a single class or all assigned classes
798
799        Parameters
800        ----------
801        classe : str, optional
802            Class name or None to plot all assigned classes, by default None
803        max_c : int, optional
804            Max Carbon Number, by default 50
805        max_dbe : int, optional
806            Max DBE, by default 40
807        dbe_incr : int, optional
808            DBE increment, by default 5
809        c_incr : int, optional
810            Carbon Number increment, by default 10
811        color : str, optional
812            Matplotlib color/colormap, by default "viridis"
813        alpha : float, optional
814            Transparency of points, by default 0.5
815        log_abundance : bool, optional
816            If True, use log10 scale for abundance values, by default False
817
818        Returns
819        -------
820        ax : matplotlib.axes
821            Matplotlib axes object
822        abun_perc : float or None
823            Class percentile of the relative abundance (if classe specified)
824        """
825        import numpy as np
826        ax = plt.gca()
827
828        if classe is not None and classe != Labels.unassigned:
829            # Single class plot
830            abun_perc = self.abundance_count_percentile(classe)
831            carbon_number = self.carbon_number(classe)
832            dbe = self.dbe(classe)
833            abundance = self.abundance(classe)
834
835            if log_abundance:
836                abundance = [np.log10(a + 1e-10) for a in abundance]
837                colorbar_label = 'log\u2081\u2080(Abundance)'
838            else:
839                colorbar_label = 'Abundance'
840
841            sorted_indices = sorted(range(len(abundance)), key=lambda i: abundance[i])
842            carbon_number = [carbon_number[i] for i in sorted_indices]
843            dbe = [dbe[i] for i in sorted_indices]
844            abundance = [abundance[i] for i in sorted_indices]
845
846            scatter = ax.scatter(carbon_number, dbe, c=abundance, alpha=alpha, cmap=color)
847            plt.colorbar(scatter, label=colorbar_label)
848
849            title = "%s, %.2f %%" % (classe, abun_perc)
850            ax.set_title(title)
851
852            return_val = ax, abun_perc
853        else:
854            # All assigned classes plot
855            carbon_number = self.carbon_number_all()
856            dbe = self.dbe_all()
857            abundance = self.abundance_assigned()
858
859            if log_abundance:
860                abundance = [np.log10(a + 1e-10) for a in abundance]
861                colorbar_label = 'log\u2081\u2080(Abundance)'
862            else:
863                colorbar_label = 'Abundance'
864
865            sorted_indices = sorted(range(len(abundance)), key=lambda i: abundance[i])
866            carbon_number = [carbon_number[i] for i in sorted_indices]
867            dbe = [dbe[i] for i in sorted_indices]
868            abundance = [abundance[i] for i in sorted_indices]
869
870            scatter = ax.scatter(carbon_number, dbe, c=abundance, alpha=alpha, cmap=color)
871            plt.colorbar(scatter, label=colorbar_label)
872
873            ax.set_title("DBE vs Carbon Number - All Assigned Classes")
874
875            return_val = ax
876
877        ax.set_xlabel("Carbon number", fontsize=16)
878        ax.set_ylabel("DBE", fontsize=16)
879        ax.tick_params(axis="both", which="major", labelsize=18)
880        ax.set_xticks(range(0, max_c, c_incr))
881        ax.set_yticks(range(0, max_dbe, dbe_incr))
882        ax.grid(alpha=0.3, linestyle='--')
883
884        return return_val