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, max_hc=2.5, max_oc=2, ticks_number=5, color="viridis"
703    ):
704        """Plot Van Krevelen Diagram
705
706        Parameters
707        ----------
708        classe : str
709            Class name
710        max_hc : float, optional
711            Max H/C ratio, by default 2.5
712        max_oc : float, optional
713            Max O/C ratio, by default 2
714        ticks_number : int, optional
715            Number of ticks, by default 5
716        color : str, optional
717            Matplotlib color, by default "viridis"
718
719        Returns
720        -------
721        ax : matplotlib.axes
722            Matplotlib axes object
723        abun_perc : float
724            Class percentile of the relative abundance
725        """
726        if classe != Labels.unassigned:
727            # get data
728            abun_perc = self.abundance_count_percentile(classe)
729            hc = self.atoms_ratio(classe, "H", "C")
730            oc = self.atoms_ratio(classe, "O", "C")
731            abundance = self.abundance(classe)
732
733            # plot data
734            ax = plt.gca()
735
736            ax.scatter(oc, hc, c=abundance, alpha=0.5, cmap=color)
737
738            # ax.scatter(carbon_number, dbe, c=color, alpha=0.5)
739
740            title = "%s, %.2f %%" % (classe, abun_perc)
741            ax.set_title(title)
742            ax.set_xlabel("O/C", fontsize=16)
743            ax.set_ylabel("H/C", fontsize=16)
744            ax.tick_params(axis="both", which="major", labelsize=18)
745            ax.set_xticks(linspace(0, max_oc, ticks_number, endpoint=True))
746            ax.set_yticks(linspace(0, max_hc, ticks_number, endpoint=True))
747
748            # returns matplot axes obj and the class percentile of the relative abundance
749
750            return ax, abun_perc
751
752    def plot_dbe_vs_carbon_number(
753        self, classe, max_c=50, max_dbe=40, dbe_incr=5, c_incr=10, color="viridis"
754    ):
755        """Plot DBE vs Carbon Number
756
757        Parameters
758        ----------
759        classe : str
760            Class name
761        max_c : int, optional
762            Max Carbon Number, by default 50
763        max_dbe : int, optional
764            Max DBE, by default 40
765        dbe_incr : int, optional
766            DBE increment, by default 5
767        c_incr : int, optional
768            Carbon Number increment, by default 10
769        color : str, optional
770            Matplotlib color, by default "viridis"
771
772        Returns
773        -------
774        ax : matplotlib.axes
775            Matplotlib axes object
776        abun_perc : float
777            Class percentile of the relative abundance
778        """
779        if classe != Labels.unassigned:
780            # get data
781            abun_perc = self.abundance_count_percentile(classe)
782            carbon_number = self.carbon_number(classe)
783            dbe = self.dbe(classe)
784            abundance = self.abundance(classe)
785
786            # plot data
787            ax = plt.gca()
788
789            ax.scatter(carbon_number, dbe, c=abundance, alpha=0.5, cmap=color)
790
791            # ax.scatter(carbon_number, dbe, c=color, alpha=0.5)
792
793            title = "%s, %.2f %%" % (classe, abun_perc)
794            ax.set_title(title)
795            ax.set_xlabel("Carbon number", fontsize=16)
796            ax.set_ylabel("DBE", fontsize=16)
797            ax.tick_params(axis="both", which="major", labelsize=18)
798            ax.set_xticks(range(0, max_c, c_incr))
799            ax.set_yticks(range(0, max_dbe, dbe_incr))
800
801            # returns matplot axes obj and the class percentile of the relative abundance
802
803            return ax, abun_perc