corems.chroma_peak.factory.chroma_peak_classes

   1__author__ = "Yuri E. Corilo"
   2__date__ = "Jun 12, 2019"
   3
   4import matplotlib.pyplot as plt
   5import numpy as np
   6import pandas as pd
   7import copy
   8
   9from corems.chroma_peak.calc.ChromaPeakCalc import (
  10    GCPeakCalculation,
  11    LCMSMassFeatureCalculation,
  12)
  13from corems.encapsulation.plot_utils import _finalize_plot
  14from corems.mass_spectra.factory.chromat_data import EIC_Data
  15from corems.molecular_id.factory.EI_SQL import LowResCompoundRef
  16
  17
  18class ChromaPeakBase:
  19    """Base class for chromatographic peak (ChromaPeak) objects.
  20
  21    Parameters
  22    -------
  23    chromatogram_parent : Chromatogram
  24        The parent chromatogram object.
  25    mass_spectrum_obj : MassSpectrum
  26        The mass spectrum object.
  27    start_index : int
  28        The start index of the peak.
  29    index : int
  30        The index of the peak.
  31    final_index : int
  32        The final index of the peak.
  33
  34    Attributes
  35    --------
  36    start_scan : int
  37        The start scan of the peak.
  38    final_scan : int
  39        The final scan of the peak.
  40    apex_scan : int
  41        The apex scan of the peak.
  42    chromatogram_parent : Chromatogram
  43        The parent chromatogram object.
  44    mass_spectrum : MassSpectrum
  45        The mass spectrum object.
  46    _area : float
  47        The area of the peak.
  48
  49    Properties
  50    --------
  51    * retention_time : float.
  52        The retention time of the peak.
  53    * tic : float.
  54        The total ion current of the peak.
  55    * area : float.
  56        The area of the peak.
  57    * rt_list : list.
  58        The list of retention times within the peak.
  59    * tic_list : list.
  60        The list of total ion currents within the peak.
  61
  62    Methods
  63    --------
  64    * None
  65    """
  66
  67    def __init__(
  68        self, chromatogram_parent, mass_spectrum_obj, start_index, index, final_index
  69    ):
  70        self.start_scan = start_index
  71        self.final_scan = final_index
  72        self.apex_scan = int(index)
  73        self.chromatogram_parent = chromatogram_parent
  74        self.mass_spectrum = mass_spectrum_obj
  75        self._area = None
  76
  77    @property
  78    def retention_time(self):
  79        """Retention Time"""
  80        return self.mass_spectrum.retention_time
  81
  82    @property
  83    def tic(self):
  84        """Total Ion Current"""
  85        return self.mass_spectrum.tic
  86
  87    @property
  88    def area(self):
  89        """Peak Area"""
  90        return self._area
  91
  92    @property
  93    def rt_list(self):
  94        """Retention Time List"""
  95        return [
  96            self.chromatogram_parent.retention_time[i]
  97            for i in range(self.start_scan, self.final_scan + 1)
  98        ]
  99
 100    @property
 101    def tic_list(self):
 102        """Total Ion Current List"""
 103        return [
 104            self.chromatogram_parent.tic[i]
 105            for i in range(self.start_scan, self.final_scan + 1)
 106        ]
 107
 108
 109class LCMSMassFeature(ChromaPeakBase, LCMSMassFeatureCalculation):
 110    """Class representing a mass feature in a liquid chromatography (LC) chromatogram.
 111
 112    Parameters
 113    -------
 114    lcms_parent : LCMS
 115        The parent LCMSBase object.
 116    mz : float
 117        The observed mass to charge ratio of the feature.
 118    retention_time : float
 119        The retention time of the feature (in minutes), at the apex.
 120    intensity : float
 121        The intensity of the feature.
 122    apex_scan : int
 123        The scan number of the apex of the feature.
 124    persistence : float, optional
 125        The persistence of the feature. Default is None.
 126
 127    Attributes
 128    --------
 129    _mz_exp : float
 130        The observed mass to charge ratio of the feature.
 131    _mz_cal : float
 132        The calibrated mass to charge ratio of the feature.
 133    _retention_time : float
 134        The retention time of the feature (in minutes), at the apex.
 135    _apex_scan : int
 136        The scan number of the apex of the feature.
 137    _intensity : float
 138        The intensity of the feature.
 139    _persistence : float
 140        The persistence of the feature.
 141    _eic_data : EIC_Data
 142        The EIC data object associated with the feature.
 143    _eic_mz : float
 144        The m/z value used to extract the EIC data,
 145        sometimes different from the observed m/z due to calibration, centroiding, or other processing.
 146    _dispersity_index : float
 147        The dispersity index of the feature, in minutes.
 148    _normalized_dispersity_index : float
 149        The normalized dispersity index of the feature (unitless, fraction of total window used to calculate dispersity index).
 150    _half_height_width : numpy.ndarray
 151        The half height width of the feature (in minutes, as an array of min and max values).
 152    _tailing_factor : float
 153        The tailing factor of the feature.
 154        > 1 indicates tailing, < 1 indicates fronting, = 1 indicates symmetrical peak.
 155    _noise_score : tuple
 156        The noise score of the feature, as a tuple of (left, right) scores.
 157        Each score is a float, with higher values indicating better signal to noise.
 158    _gaussian_similarity : float
 159        The Gaussian similarity of the feature, as a float between 0 and 1.
 160        1 indicates a perfect Gaussian shape, 0 indicates a non-Gaussian shape.
 161    _ms_deconvoluted_idx : [int]
 162        The indexes of the mass_spectrum attribute in the deconvoluted mass spectrum.
 163    _type : str
 164        The type of mass feature. Default is "untargeted".
 165        Can be "untargeted", "targeted", or another customized type.
 166    is_calibrated : bool
 167        If True, the feature has been calibrated. Default is False.
 168    monoisotopic_mf_id : int
 169        Mass feature id that is the monoisotopic version of self.
 170        If self.id, then self is the monoisotopic feature). Default is None.
 171    isotopologue_type : str
 172        The isotopic class of the feature, i.e. "13C1", "13C2", "13C1 37Cl1" etc.
 173        Default is None.
 174    ms2_scan_numbers : list
 175        List of scan numbers of the MS2 spectra associated with the feature.
 176        Default is an empty list.
 177    ms2_mass_spectra : dict
 178        Dictionary of MS2 spectra associated with the feature (key = scan number for DDA).
 179        Default is an empty dictionary.
 180    ms2_similarity_results : list
 181        List of MS2 similarity results associated with the mass feature.
 182        Default is an empty list.
 183    id : int
 184        The ID of the feature, also the key in the parent LCMS object's
 185        `mass_features` dictionary.
 186    mass_spectrum_deconvoluted_parent : bool
 187        If True, the mass feature corresponds to the most intense peak in the deconvoluted mass spectrum. Default is None.
 188    associated_mass_features_deconvoluted : list
 189        List of mass features associated with the deconvoluted mass spectrum. Default is an empty list.
 190
 191    """
 192
 193    def __init__(
 194        self,
 195        lcms_parent,
 196        mz: float,
 197        retention_time: float,
 198        intensity: float,
 199        apex_scan: int,
 200        persistence: float = None,
 201        id: int = None,
 202    ):
 203        super().__init__(
 204            chromatogram_parent=lcms_parent,
 205            mass_spectrum_obj=None,
 206            start_index=None,
 207            index=apex_scan,
 208            final_index=None,
 209        )
 210        # Core attributes, marked as private
 211        self._mz_exp: float = mz
 212        self._mz_cal: float = None
 213        self._retention_time: float = retention_time
 214        self._apex_scan: int = apex_scan
 215        self._intensity: float = intensity
 216        self._persistence: float = persistence
 217        self._eic_data: EIC_Data = None
 218        # m/z used for EIC extraction; set on integrate (may differ from mz after cal)
 219        self._eic_mz: float = None
 220        self._dispersity_index: float = None
 221        self._normalized_dispersity_index: float = None
 222        self._half_height_width: np.ndarray = None
 223        self._ms_deconvoluted_idx = None
 224        self._tailing_factor: float = None
 225        self._noise_score: tuple = None
 226        self._gaussian_similarity: float = None
 227        self._type: str = "untargeted"
 228
 229        # Additional attributes
 230        self.monoisotopic_mf_id = None
 231        self.isotopologue_type = None
 232        self.ms2_scan_numbers = []
 233        self.ms2_mass_spectra = {}
 234        self.ms2_similarity_results = []
 235        self.mass_spectrum_deconvoluted_parent: bool = None
 236        self.associated_mass_features_deconvoluted = []
 237
 238        if id:
 239            self.id = id
 240        else:
 241            # get the parent's mass feature keys and add 1 to the max value to get the new key
 242            self.id = (
 243                max(lcms_parent.mass_features.keys()) + 1
 244                if lcms_parent.mass_features.keys()
 245                else 0
 246            )
 247
 248    def update_mz(self):
 249        """Update the mass to charge ratio from the mass spectrum object."""
 250        if self.mass_spectrum is None:
 251            raise ValueError(
 252                "The mass spectrum object is not set, cannot update the m/z from the MassSpectrum object"
 253            )
 254        if len(self.mass_spectrum.mz_exp) == 0:
 255            raise ValueError(
 256                "The mass spectrum object has no m/z values, cannot update the m/z from the MassSpectrum object until it is processed"
 257            )
 258        new_mz = self.ms1_peak.mz_exp
 259
 260        # calculate the difference between the new and old m/z, only update if it is close
 261        mz_diff = new_mz - self.mz
 262        if abs(mz_diff) < 0.01:
 263            self._mz_exp = new_mz
 264
 265    def _plot_ms1_spectrum(self, ax, deconvoluted=False, sample_name=None):
 266        """Internal method to plot MS1 spectrum on a given axis.
 267
 268        Parameters
 269        ----------
 270        ax : matplotlib.axes.Axes
 271            The axis to plot on.
 272        deconvoluted : bool, optional
 273            If True and deconvoluted spectrum exists, plot both raw and deconvoluted. Default is False.
 274        sample_name : str, optional
 275            Sample name to include in title. Default is None.
 276        """
 277        if self.mass_spectrum is None:
 278            raise ValueError("MS1 spectrum is not available")
 279
 280        title_prefix = "MS1 (deconvoluted)" if deconvoluted else "MS1 (raw)"
 281        if sample_name:
 282            ax.set_title(f"{title_prefix} - {sample_name}", loc="left")
 283        else:
 284            ax.set_title(title_prefix, loc="left")
 285
 286        if deconvoluted and self._ms_deconvoluted_idx is not None:
 287            # Plot both raw and deconvoluted
 288            ax.vlines(
 289                self.mass_spectrum.mz_exp,
 290                0,
 291                self.mass_spectrum.abundance,
 292                color="k",
 293                alpha=0.2,
 294                label="Raw MS1",
 295            )
 296            ax.vlines(
 297                self.mass_spectrum_deconvoluted.mz_exp,
 298                0,
 299                self.mass_spectrum_deconvoluted.abundance,
 300                color="k",
 301                label="Deconvoluted MS1",
 302            )
 303            ax.set_xlim(
 304                self.mass_spectrum_deconvoluted.mz_exp.min() * 0.8,
 305                self.mass_spectrum_deconvoluted.mz_exp.max() * 1.1,
 306            )
 307            ax.set_ylim(0, self.mass_spectrum_deconvoluted.abundance.max() * 1.1)
 308        else:
 309            # Plot raw only
 310            ax.vlines(
 311                self.mass_spectrum.mz_exp,
 312                0,
 313                self.mass_spectrum.abundance,
 314                color="k",
 315                label="Raw MS1",
 316            )
 317            ax.set_xlim(
 318                self.mass_spectrum.mz_exp.min() * 0.8,
 319                self.mass_spectrum.mz_exp.max() * 1.1,
 320            )
 321            ax.set_ylim(bottom=0)
 322
 323        # Highlight the feature m/z if close enough
 324        if abs(self.ms1_peak.mz_exp - self.mz) < 0.01:
 325            ax.vlines(
 326                self.ms1_peak.mz_exp,
 327                0,
 328                self.ms1_peak.abundance,
 329                color="m",
 330                label="Feature m/z",
 331            )
 332        else:
 333            if self.chromatogram_parent.parameters.lc_ms.verbose_processing:
 334                print(
 335                    f"The m/z of the mass feature {self.id} is different from the m/z of MS1 peak, "
 336                    "the MS1 peak will not be plotted"
 337                )
 338
 339        ax.legend(loc="upper left")
 340        ax.set_ylabel("Intensity")
 341        ax.set_xlabel("m/z")
 342        # Combining tick_params(labelleft=False) with set_title(loc="left") makes
 343        # tight_layout() produce NaN axis positions on matplotlib 3.11/numpy 2.5.
 344        ax.set_yticklabels([])
 345
 346    def _plot_ms2_spectrum(self, ax, sample_name=None):
 347        """Internal method to plot MS2 spectrum on a given axis.
 348
 349        Parameters
 350        ----------
 351        ax : matplotlib.axes.Axes
 352            The axis to plot on.
 353        sample_name : str, optional
 354            Sample name to include in title. Default is None.
 355        """
 356        if len(self.ms2_mass_spectra) == 0:
 357            raise ValueError("MS2 spectrum is not available")
 358
 359        if sample_name:
 360            ax.set_title(f"MS2 - {sample_name}", loc="left")
 361        else:
 362            ax.set_title("MS2", loc="left")
 363
 364        ax.vlines(self.best_ms2.mz_exp, 0, self.best_ms2.abundance, color="k")
 365        ax.set_ylabel("Intensity")
 366        ax.set_xlabel("m/z")
 367        ax.set_ylim(bottom=0)
 368        ax.yaxis.get_major_formatter().set_scientific(False)
 369        ax.yaxis.get_major_formatter().set_useOffset(False)
 370
 371    def _plot_ms2_mirror(self, ax, molecular_metadata=None, spectral_library=None):
 372        """Internal method to plot MS2 mirror spectrum on a given axis.
 373
 374        Plots experimental MS2 on top (positive) and library MS2 on bottom (negative/mirrored)
 375        if MS2 similarity results are available. If no MS2 similarity results exist,
 376        falls back to regular MS2 plot.
 377
 378        Parameters
 379        ----------
 380        ax : matplotlib.axes.Axes
 381            The axis to plot on.
 382        molecular_metadata : dict, optional
 383            Dictionary mapping molecular IDs to MetaboliteMetadata objects.
 384            If provided, uses metadata for compound names.
 385            Default is None.
 386        spectral_library : FlashEntropySearch or list of FlashEntropySearch, optional
 387            FlashEntropy spectral library (or list of libraries) containing MS2 spectra.
 388            If provided, uses library to retrieve MS2 spectra by ref_ms_id.
 389            Default is None.
 390
 391        Raises
 392        ------
 393        ValueError
 394            If MS2 similarity results exist but molecular_metadata or spectral_library is None.
 395        """
 396        if len(self.ms2_mass_spectra) == 0:
 397            ax.text(
 398                0.5,
 399                0.5,
 400                "No MS2 data available",
 401                ha="center",
 402                va="center",
 403                transform=ax.transAxes,
 404                fontsize=12,
 405            )
 406            ax.set_xlabel("m/z", fontsize=10)
 407            ax.set_ylabel("Relative Intensity (%)", fontsize=10)
 408            return
 409
 410        # Check if we have MS2 similarity results - if not, fall back to regular MS2 plot
 411        if len(self.ms2_similarity_results) == 0:
 412            self._plot_ms2_spectrum(ax)
 413            return
 414
 415        # If we have MS2 similarity results, we need both molecular_metadata and spectral_library
 416        if molecular_metadata is None or spectral_library is None:
 417            raise ValueError(
 418                "MS2 mirror plot requires both 'molecular_metadata' and 'spectral_library' "
 419                "parameters when MS2 similarity results are present. "
 420                "Please provide both parameters to plot_cluster() or plot()."
 421            )
 422
 423        # Get experimental MS2
 424        sample_ms2 = self.best_ms2
 425        sample_mz = sample_ms2.mz_exp
 426        sample_int = sample_ms2.abundance
 427
 428        # Normalize sample MS2
 429        if len(sample_int) > 0 and max(sample_int) > 0:
 430            sample_int = sample_int / max(sample_int) * 100
 431
 432        # Plot sample MS2 on top (positive)
 433        ax.vlines(
 434            sample_mz,
 435            0,
 436            sample_int,
 437            colors="blue",
 438            alpha=0.7,
 439            linewidths=1.5,
 440            label="Sample MS2",
 441        )
 442
 443        # Check if we have MS2 similarity results
 444        library_ms2_peaks = None
 445        entropy_similarity = None
 446        molecule_name = None
 447        mol_id = None
 448
 449        if len(self.ms2_similarity_results) > 0:
 450            # Get all results as dataframes and find the best match
 451            results_df = [x.to_dataframe() for x in self.ms2_similarity_results]
 452            results_df = pd.concat(results_df)
 453            results_df = results_df.sort_values(
 454                by="entropy_similarity", ascending=False
 455            )
 456
 457            # Get the best match
 458            best_result = results_df.iloc[0]
 459            entropy_similarity = best_result["entropy_similarity"]
 460            mol_id = best_result.get("ref_mol_id", None)
 461            ref_ms_id = best_result.get("ref_ms_id", None)
 462
 463            # Get library spectrum from spectral_library using ref_ms_id
 464            if spectral_library is not None and ref_ms_id is not None:
 465                # Handle both single library and list of libraries
 466                libraries = (
 467                    spectral_library
 468                    if isinstance(spectral_library, list)
 469                    else [spectral_library]
 470                )
 471
 472                # Search through all libraries to find the ref_ms_id
 473                for library in libraries:
 474                    try:
 475                        # Get the IDs in the spectral library
 476                        fe_spec_index = [x["id"] for x in library].index(ref_ms_id)
 477                        library_ms2_peaks = library[fe_spec_index]["peaks"]
 478                        break  # Found the spectrum, exit the loop
 479                    except ValueError:
 480                        # ref_ms_id not found in this library, continue to next
 481                        continue
 482
 483                # If ref_ms_id was not found in any library, raise an error
 484                if library_ms2_peaks is None:
 485                    raise ValueError(
 486                        f"Reference MS ID '{ref_ms_id}' not found in any of the provided spectral libraries. "
 487                        f"Please ensure the spectral library contains the matching reference spectrum."
 488                    )
 489
 490            # Get compound name from molecular_metadata using mol_id
 491            if molecular_metadata is not None and mol_id is not None:
 492                if mol_id in molecular_metadata:
 493                    metadata = molecular_metadata[mol_id]
 494                    # Get compound name from metadata
 495                    molecule_name = getattr(
 496                        metadata, "common_name", getattr(metadata, "name", "Unknown")
 497                    )
 498
 499        # Plot library MS2 on bottom (negative/mirrored)
 500        if library_ms2_peaks is not None and len(library_ms2_peaks) > 0:
 501            lib_mz = library_ms2_peaks[:, 0]
 502            lib_int = library_ms2_peaks[:, 1]
 503            # Normalize
 504            if len(lib_int) > 0 and max(lib_int) > 0:
 505                lib_int = lib_int / max(lib_int) * 100
 506            # Mirror to negative
 507            lib_int_mirror = -lib_int
 508
 509            # Create label with molecule name and molecular ID
 510            lib_label = "Library MS2"
 511            if molecule_name:
 512                lib_label += f" ({molecule_name})"
 513            if mol_id:
 514                lib_label += f" [ID: {mol_id}]"
 515
 516            ax.vlines(
 517                lib_mz,
 518                0,
 519                lib_int_mirror,
 520                colors="red",
 521                alpha=0.7,
 522                linewidths=1.5,
 523                label=lib_label,
 524            )
 525
 526        ax.axhline(0, color="black", linewidth=0.5)
 527        ax.set_xlabel("m/z", fontsize=10)
 528        ax.set_ylabel("Relative Intensity (%)", fontsize=10)
 529        ax.legend(fontsize=8, loc="upper right")
 530        ax.grid(True, alpha=0.3)
 531
 532        # Set y-axis to symmetric range
 533        ax.set_ylim(-105, 105)
 534
 535        # Add entropy similarity to the title if available
 536        if entropy_similarity is not None:
 537            ax.set_title(
 538                f"MS2 Mirror Plot (Entropy Similarity: {entropy_similarity:.3f})",
 539                loc="left",
 540            )
 541        else:
 542            ax.set_title("MS2 Mirror Plot", loc="left")
 543
 544    def _plot_single_eic(
 545        self,
 546        ax,
 547        plot_smoothed=False,
 548        plot_datapoints=False,
 549        eic_buffer_time=None,
 550        show_ms2_scan=True,
 551    ):
 552        """Internal method to plot a single EIC on a given axis.
 553
 554        Parameters
 555        ----------
 556        ax : matplotlib.axes.Axes
 557            The axis to plot on.
 558        plot_smoothed : bool, optional
 559            If True, plot smoothed EIC. Default is False.
 560        plot_datapoints : bool, optional
 561            If True, plot EIC datapoints. Default is False.
 562        eic_buffer_time : float, optional
 563            Time buffer around the peak (minutes). If None, uses parameter setting. Default is None.
 564        show_ms2_scan : bool, optional
 565            If True and MS2 scans exist, show vertical line at MS2 scan time. Default is True.
 566        """
 567        if self._eic_data is None:
 568            raise ValueError("EIC data is not available")
 569
 570        if eic_buffer_time is None:
 571            eic_buffer_time = self.chromatogram_parent.parameters.lc_ms.eic_buffer_time
 572
 573        ax.set_title("EIC", loc="left")
 574        ax.plot(self._eic_data.time, self._eic_data.eic, c="tab:blue", label="EIC")
 575
 576        if plot_datapoints:
 577            ax.scatter(
 578                self._eic_data.time,
 579                self._eic_data.eic,
 580                c="tab:blue",
 581                label="EIC Data Points",
 582            )
 583
 584        if plot_smoothed and hasattr(self._eic_data, "eic_smoothed"):
 585            ax.plot(
 586                self._eic_data.time,
 587                self._eic_data.eic_smoothed,
 588                c="tab:red",
 589                label="Smoothed EIC",
 590            )
 591
 592        # Fill integrated area if available
 593        if self.start_scan is not None:
 594            ax.fill_between(self.eic_rt_list, self.eic_list, color="b", alpha=0.2)
 595        else:
 596            if self.chromatogram_parent.parameters.lc_ms.verbose_processing:
 597                print(
 598                    f"No start and final scan numbers were provided for mass feature {self.id}"
 599                )
 600
 601        ax.set_ylabel("Intensity")
 602        ax.set_xlabel("Time (minutes)")
 603        ax.set_ylim(0, self.eic_list.max() * 1.1)
 604        ax.set_xlim(
 605            self.retention_time - eic_buffer_time,
 606            self.retention_time + eic_buffer_time,
 607        )
 608        ax.axvline(x=self.retention_time, color="k", label="MS1 scan time (apex)")
 609
 610        # Show MS2 scan time if available and requested
 611        if show_ms2_scan and len(self.ms2_scan_numbers) > 0:
 612            ax.axvline(
 613                x=self.chromatogram_parent.get_time_of_scan_id(
 614                    self.best_ms2.scan_number
 615                ),
 616                color="grey",
 617                linestyle="--",
 618                label="MS2 scan time",
 619            )
 620
 621        ax.legend(loc="upper left")
 622        ax.yaxis.get_major_formatter().set_useOffset(False)
 623
 624    def plot(
 625        self,
 626        to_plot=["EIC", "MS1", "MS2"],
 627        return_fig=True,
 628        path=None,
 629        plot_smoothed_eic=False,
 630        plot_eic_datapoints=False,
 631        molecular_metadata=None,
 632        spectral_library=None,
 633    ):
 634        """Plot the mass feature.
 635
 636        Parameters
 637        ----------
 638        to_plot : list, optional
 639            List of strings specifying what to plot, any iteration of
 640            "EIC", "MS2", "MS2_mirror", and "MS1".
 641            Default is ["EIC", "MS1", "MS2"].
 642        return_fig : bool, optional
 643            If True, return the open figure (caller owns lifecycle).
 644            Default is True.
 645        path : str or path-like, optional
 646            If set, save the figure to this path. When ``return_fig`` is False,
 647            the figure is closed after saving and ``plt.show()`` is not called.
 648        plot_smoothed_eic : bool, optional
 649            If True, the smoothed EIC is plotted. Default is False.
 650        plot_eic_datapoints : bool, optional
 651            If True, the EIC data points are plotted. Default is False.
 652        molecular_metadata : dict, optional
 653            Dictionary mapping molecular IDs to MetaboliteMetadata objects.
 654            Required if "MS2_mirror" is in to_plot. Default is None.
 655        spectral_library : FlashEntropySearch, optional
 656            FlashEntropy spectral library containing MS2 spectra.
 657            Required if "MS2_mirror" is in to_plot. Default is None.
 658
 659        Returns
 660        -------
 661        matplotlib.figure.Figure or None
 662            The figure object if `return_fig` is True (left open for the caller).
 663            Otherwise None; the figure is displayed with ``plt.show()`` when
 664            ``path`` is not set.
 665        """
 666        # Adjust to_plot list if there are not spectra added to the mass features
 667        if self.mass_spectrum is None:
 668            to_plot = [x for x in to_plot if x != "MS1"]
 669        if len(self.ms2_mass_spectra) == 0:
 670            to_plot = [x for x in to_plot if x not in ["MS2", "MS2_mirror"]]
 671        if self._eic_data is None:
 672            to_plot = [x for x in to_plot if x != "EIC"]
 673
 674        # Check if MS2_mirror is requested without molecular_metadata
 675        if "MS2_mirror" in to_plot and molecular_metadata is None:
 676            raise ValueError(
 677                "molecular_metadata is required when 'MS2_mirror' is in to_plot"
 678            )
 679
 680        # Check if both MS2 and MS2_mirror are requested (not allowed)
 681        if "MS2" in to_plot and "MS2_mirror" in to_plot:
 682            # Remove regular MS2 if mirror is requested
 683            to_plot = [x for x in to_plot if x != "MS2"]
 684
 685        deconvoluted = self._ms_deconvoluted_idx is not None
 686
 687        fig, axs = plt.subplots(
 688            len(to_plot), 1, figsize=(9, len(to_plot) * 4), squeeze=False
 689        )
 690        fig.suptitle(
 691            f"Mass Feature {self.id}: m/z = {round(self.mz, ndigits=4)}; "
 692            f"time = {round(self.retention_time, ndigits=1)} minutes"
 693        )
 694
 695        i = 0
 696        # EIC plot
 697        if "EIC" in to_plot:
 698            self._plot_single_eic(
 699                axs[i][0],
 700                plot_smoothed=plot_smoothed_eic,
 701                plot_datapoints=plot_eic_datapoints,
 702            )
 703            i += 1
 704
 705        # MS1 plot
 706        if "MS1" in to_plot:
 707            self._plot_ms1_spectrum(axs[i][0], deconvoluted=deconvoluted)
 708            i += 1
 709
 710        # MS2 plot
 711        if "MS2" in to_plot:
 712            self._plot_ms2_spectrum(axs[i][0])
 713            i += 1
 714
 715        # MS2 mirror plot
 716        if "MS2_mirror" in to_plot:
 717            self._plot_ms2_mirror(
 718                axs[i][0],
 719                molecular_metadata=molecular_metadata,
 720                spectral_library=spectral_library,
 721            )
 722            i += 1
 723
 724        # Add space between subplots
 725        plt.tight_layout()
 726
 727        return _finalize_plot(fig, return_fig=return_fig, path=path)
 728
 729    @property
 730    def mz(self):
 731        """Mass to charge ratio of the mass feature"""
 732        # If the mass feature has been calibrated, return the calibrated m/z, otherwise return the measured m/z
 733        if self._mz_cal is not None:
 734            return self._mz_cal
 735        else:
 736            return self._mz_exp
 737
 738    @property
 739    def mass_spectrum_deconvoluted(self):
 740        """Returns the deconvoluted mass spectrum object associated with the mass feature, if deconvolution has been performed."""
 741        if self._ms_deconvoluted_idx is not None:
 742            ms_deconvoluted = copy.deepcopy(self.mass_spectrum)
 743            ms_deconvoluted.set_indexes(self._ms_deconvoluted_idx)
 744            return ms_deconvoluted
 745        else:
 746            raise ValueError(
 747                "Deconvolution has not been performed for mass feature " + str(self.id)
 748            )
 749
 750    @property
 751    def retention_time(self):
 752        """Retention time of the mass feature"""
 753        return self._retention_time
 754
 755    @retention_time.setter
 756    def retention_time(self, value):
 757        """Set the retention time of the mass feature"""
 758        if not isinstance(value, float):
 759            raise ValueError("The retention time of the mass feature must be a float")
 760        self._retention_time = value
 761
 762    @property
 763    def apex_scan(self):
 764        """Apex scan of the mass feature"""
 765        return self._apex_scan
 766
 767    @apex_scan.setter
 768    def apex_scan(self, value):
 769        """Set the apex scan of the mass feature"""
 770        if not isinstance(value, int):
 771            raise ValueError("The apex scan of the mass feature must be an integer")
 772        self._apex_scan = value
 773
 774    @property
 775    def intensity(self):
 776        """Intensity of the mass feature"""
 777        return self._intensity
 778
 779    @intensity.setter
 780    def intensity(self, value):
 781        """Set the intensity of the mass feature"""
 782        if not isinstance(value, float):
 783            raise ValueError("The intensity of the mass feature must be a float")
 784        self._intensity = value
 785
 786    @property
 787    def persistence(self):
 788        """Persistence of the mass feature"""
 789        return self._persistence
 790
 791    @persistence.setter
 792    def persistence(self, value):
 793        """Set the persistence of the mass feature"""
 794        if not isinstance(value, float):
 795            raise ValueError("The persistence of the mass feature must be a float")
 796        self._persistence = value
 797
 798    @property
 799    def eic_rt_list(self):
 800        """Retention time list between the beginning and end of the mass feature"""
 801        # Find index of the start and final scans in the EIC data
 802        start_index = self._eic_data.scans.tolist().index(self.start_scan)
 803        final_index = self._eic_data.scans.tolist().index(self.final_scan)
 804
 805        # Get the retention time list
 806        rt_list = self._eic_data.time[start_index : final_index + 1]
 807        return rt_list
 808
 809    @property
 810    def eic_list(self):
 811        """EIC List between the beginning and end of the mass feature"""
 812        # Find index of the start and final scans in the EIC data
 813        start_index = self._eic_data.scans.tolist().index(self.start_scan)
 814        final_index = self._eic_data.scans.tolist().index(self.final_scan)
 815
 816        # Get the retention time list
 817        eic = self._eic_data.eic[start_index : final_index + 1]
 818        return eic
 819
 820    @property
 821    def ms1_peak(self):
 822        """MS1 peak from associated mass spectrum that is closest to the mass feature's m/z"""
 823        # Find index array self.mass_spectrum.mz_exp that is closest to self.mz
 824        closest_mz = min(self.mass_spectrum.mz_exp, key=lambda x: abs(x - self.mz))
 825        closest_mz_index = self.mass_spectrum.mz_exp.tolist().index(closest_mz)
 826
 827        return self.mass_spectrum._mspeaks[closest_mz_index]
 828
 829    @property
 830    def tailing_factor(self):
 831        """Tailing factor of the mass feature"""
 832        return self._tailing_factor
 833
 834    @tailing_factor.setter
 835    def tailing_factor(self, value):
 836        """Set the tailing factor of the mass feature"""
 837        if not isinstance(value, float):
 838            raise ValueError("The tailing factor of the mass feature must be a float")
 839        self._tailing_factor = value
 840
 841    @property
 842    def dispersity_index(self):
 843        """Dispersity index of the mass feature"""
 844        return self._dispersity_index
 845
 846    @dispersity_index.setter
 847    def dispersity_index(self, value):
 848        """Set the dispersity index of the mass feature"""
 849        if not isinstance(value, float):
 850            raise ValueError("The dispersity index of the mass feature must be a float")
 851        self._dispersity_index = value
 852
 853    @property
 854    def normalized_dispersity_index(self):
 855        """Normalized dispersity index of the mass feature, unitless (fraction of total window used)"""
 856        return self._normalized_dispersity_index
 857
 858    @property
 859    def half_height_width(self):
 860        """Half height width of the mass feature, average of min and max values, in minutes"""
 861        return np.mean(self._half_height_width)
 862
 863    @property
 864    def noise_score(self):
 865        """Mean of left and right noise scores.
 866
 867        Returns
 868        -------
 869        float or np.nan
 870            Mean noise score, or np.nan if both sides are np.nan.
 871        """
 872        if self._noise_score is None:
 873            return np.nan
 874
 875        left, right = self._noise_score
 876        # Handle NaN values
 877        if np.isnan(left) and np.isnan(right):
 878            return np.nan
 879        elif np.isnan(left):
 880            return right
 881        elif np.isnan(right):
 882            return left
 883        else:
 884            return (left + right) / 2.0
 885
 886    @property
 887    def noise_score_min(self):
 888        """Minimum of left and right noise scores.
 889
 890        Returns
 891        -------
 892        float or np.nan
 893            Minimum noise score, or np.nan if both sides are np.nan.
 894        """
 895        if self._noise_score is None:
 896            return np.nan
 897
 898        left, right = self._noise_score
 899        # Handle NaN values - nanmin ignores NaN
 900        return np.nanmin([left, right])
 901
 902    @property
 903    def noise_score_max(self):
 904        """Maximum of left and right noise scores.
 905
 906        Returns
 907        -------
 908        float or np.nan
 909            Maximum noise score, or np.nan if both sides are np.nan.
 910        """
 911        if self._noise_score is None:
 912            return np.nan
 913
 914        left, right = self._noise_score
 915        # Handle NaN values - nanmax ignores NaN
 916        return np.nanmax([left, right])
 917
 918    @property
 919    def type(self):
 920        """Type of the mass feature.
 921
 922        Returns
 923        -------
 924        str
 925            The type of mass feature ("untargeted", "targeted", or "internal standard").
 926        """
 927        return self._type
 928
 929    @type.setter
 930    def type(self, value):
 931        """Set the type of the mass feature.
 932
 933        Parameters
 934        ----------
 935        value : str
 936            The type of mass feature. Should be one of: "untargeted", "targeted", "internal standard".
 937        """
 938        if not isinstance(value, str):
 939            raise ValueError("The type of the mass feature must be a string")
 940        self._type = value
 941
 942    @property
 943    def best_ms2(self):
 944        """Points to the best representative MS2 mass spectrum
 945
 946        Notes
 947        -----
 948        If there is only one MS2 mass spectrum, it will be returned
 949        If there are MS2 similarity results, this will return the MS2 mass spectrum with the highest entropy similarity score.
 950        If there are no MS2 similarity results, the best MS2 mass spectrum is determined by the closest scan time to the apex of the mass feature, with higher resolving power.  Checks for and disqualifies possible chimeric spectra.
 951
 952        Returns
 953        -------
 954        MassSpectrum or None
 955            The best MS2 mass spectrum.
 956        """
 957        if len(self.ms2_similarity_results) > 0:
 958            # the scan number with the highest similarity score
 959            results_df = [x.to_dataframe() for x in self.ms2_similarity_results]
 960            results_df = pd.concat(results_df)
 961            results_df = results_df.sort_values(
 962                by="entropy_similarity", ascending=False
 963            )
 964            best_scan_number = results_df.iloc[0]["query_spectrum_id"]
 965            return self.ms2_mass_spectra[best_scan_number]
 966
 967        ms2_scans = list(self.ms2_mass_spectra.keys())
 968        if len(ms2_scans) > 1:
 969            mz_diff_list = []  # List of mz difference between mz of mass feature and mass of nearest mz in each scan
 970            res_list = []  # List of maximum resolving power of peaks in each scan
 971            time_diff_list = []  # List of time difference between scan and apex scan in each scan
 972            for scan in ms2_scans:
 973                if len(self.ms2_mass_spectra[scan].mspeaks) > 0:
 974                    # Find mz closest to mass feature mz, return both the difference in mass and its resolution
 975                    closest_mz = min(
 976                        self.ms2_mass_spectra[scan].mz_exp,
 977                        key=lambda x: abs(x - self.mz),
 978                    )
 979                    if all(
 980                        np.isnan(self.ms2_mass_spectra[scan].resolving_power)
 981                    ):  # All NA for resolving power in peaks, not uncommon in CID spectra
 982                        res_list.append(2)  # Assumes very low resolving power
 983                    else:
 984                        res_list.append(
 985                            np.nanmax(self.ms2_mass_spectra[scan].resolving_power)
 986                        )
 987                    mz_diff_list.append(np.abs(closest_mz - self.mz))
 988                    time_diff_list.append(
 989                        np.abs(
 990                            self.chromatogram_parent.get_time_of_scan_id(scan)
 991                            - self.retention_time
 992                        )
 993                    )
 994                else:
 995                    res_list.append(np.nan)
 996                    mz_diff_list.append(np.nan)
 997                    time_diff_list.append(np.nan)
 998            # Convert diff_lists into logical scores (higher is better for each score)
 999            time_score = 1 - np.array(time_diff_list) / np.nanmax(
1000                np.array(time_diff_list)
1001            )
1002            res_score = np.array(res_list) / np.nanmax(np.array(res_list))
1003            # mz_score is 0 for possible chimerics, 1 for all others (already within mass tolerance before assigning)
1004            mz_score = np.zeros(len(ms2_scans))
1005            for i in np.arange(0, len(ms2_scans)):
1006                if mz_diff_list[i] < 0.8 and mz_diff_list[i] > 0.1:  # Possible chimeric
1007                    mz_score[i] = 0
1008                else:
1009                    mz_score[i] = 1
1010            # get the index of the best score and return the mass spectrum
1011            if len([np.nanargmax(time_score * res_score * mz_score)]) == 1:
1012                return self.ms2_mass_spectra[
1013                    ms2_scans[np.nanargmax(time_score * res_score * mz_score)]
1014                ]
1015            # remove the mz_score condition and try again
1016            elif len(np.argmax(time_score * res_score)) == 1:
1017                return self.ms2_mass_spectra[
1018                    ms2_scans[np.nanargmax(time_score * res_score)]
1019                ]
1020            else:
1021                raise ValueError(
1022                    "No best MS2 mass spectrum could be found for mass feature "
1023                    + str(self.id)
1024                )
1025        elif len(ms2_scans) == 1:  # if only one ms2 spectra, return it
1026            return self.ms2_mass_spectra[ms2_scans[0]]
1027        else:  # if no ms2 spectra, return None
1028            return None
1029
1030
1031class GCPeak(ChromaPeakBase, GCPeakCalculation):
1032    """Class representing a peak in a gas chromatography (GC) chromatogram.
1033
1034    Parameters
1035    ----------
1036    chromatogram_parent : Chromatogram
1037        The parent chromatogram object.
1038    mass_spectrum_obj : MassSpectrum
1039        The mass spectrum object associated with the peak.
1040    indexes : tuple
1041        The indexes of the peak in the chromatogram.
1042
1043    Attributes
1044    ----------
1045    _compounds : list
1046        List of compounds associated with the peak.
1047    _ri : float or None
1048        Retention index of the peak.
1049
1050    Methods
1051    -------
1052    * __len__(). Returns the number of compounds associated with the peak.
1053    * __getitem__(position).  Returns the compound at the specified position.
1054    * remove_compound(compounds_obj). Removes the specified compound from the peak.
1055    * clear_compounds(). Removes all compounds from the peak.
1056    * add_compound(compounds_dict, spectral_similarity_scores, ri_score=None, similarity_score=None). Adds a compound to the peak with the specified attributes.
1057    * ri().  Returns the retention index of the peak.
1058    * highest_ss_compound(). Returns the compound with the highest spectral similarity score.
1059    * highest_score_compound(). Returns the compound with the highest similarity score.
1060    * compound_names(). Returns a list of names of compounds associated with the peak.
1061    """
1062
1063    def __init__(self, chromatogram_parent, mass_spectrum_obj, indexes):
1064        self._compounds = []
1065        self._ri = None
1066        super().__init__(chromatogram_parent, mass_spectrum_obj, *indexes)
1067
1068    def __len__(self):
1069        return len(self._compounds)
1070
1071    def __getitem__(self, position):
1072        return self._compounds[position]
1073
1074    def remove_compound(self, compounds_obj):
1075        self._compounds.remove(compounds_obj)
1076
1077    def clear_compounds(self):
1078        self._compounds = []
1079
1080    def add_compound(
1081        self,
1082        compounds_dict,
1083        spectral_similarity_scores,
1084        ri_score=None,
1085        similarity_score=None,
1086    ):
1087        """Adds a compound to the peak with the specified attributes.
1088
1089        Parameters
1090        ----------
1091        compounds_dict : dict
1092            Dictionary containing the compound information.
1093        spectral_similarity_scores : dict
1094            Dictionary containing the spectral similarity scores.
1095        ri_score : float or None, optional
1096            The retention index score of the compound. Default is None.
1097        similarity_score : float or None, optional
1098            The similarity score of the compound. Default is None.
1099        """
1100        compound_obj = LowResCompoundRef(compounds_dict)
1101        compound_obj.spectral_similarity_scores = spectral_similarity_scores
1102        compound_obj.spectral_similarity_score = spectral_similarity_scores.get(
1103            "cosine_correlation"
1104        )
1105        # TODO check is the above line correct?
1106        compound_obj.ri_score = ri_score
1107        compound_obj.similarity_score = similarity_score
1108        self._compounds.append(compound_obj)
1109        if similarity_score:
1110            self._compounds.sort(key=lambda c: c.similarity_score, reverse=True)
1111        else:
1112            self._compounds.sort(
1113                key=lambda c: c.spectral_similarity_score, reverse=True
1114            )
1115
1116    @property
1117    def ri(self):
1118        """Returns the retention index of the peak.
1119
1120        Returns
1121        -------
1122        float or None
1123            The retention index of the peak.
1124        """
1125        return self._ri
1126
1127    @property
1128    def highest_ss_compound(self):
1129        """Returns the compound with the highest spectral similarity score.
1130
1131        Returns
1132        -------
1133        LowResCompoundRef or None
1134            The compound with the highest spectral similarity score.
1135        """
1136        if self:
1137            return max(self, key=lambda c: c.spectral_similarity_score)
1138        else:
1139            return None
1140
1141    @property
1142    def highest_score_compound(self):
1143        """Returns the compound with the highest similarity score.
1144
1145        Returns
1146        -------
1147        LowResCompoundRef or None
1148            The compound with the highest similarity score.
1149        """
1150        if self:
1151            return max(self, key=lambda c: c.similarity_score)
1152        else:
1153            return None
1154
1155    @property
1156    def compound_names(self):
1157        """Returns a list of names of compounds associated with the peak.
1158
1159        Returns
1160        -------
1161        list
1162            List of names of compounds associated with the peak.
1163        """
1164        if self:
1165            return [c.name for c in self]
1166        else:
1167            return []
1168
1169
1170class GCPeakDeconvolved(GCPeak):
1171    """Represents a deconvolved peak in a chromatogram.
1172
1173    Parameters
1174    ----------
1175    chromatogram_parent : Chromatogram
1176        The parent chromatogram object.
1177    mass_spectra : list
1178        List of mass spectra associated with the peak.
1179    apex_index : int
1180        Index of the apex mass spectrum in the `mass_spectra` list.
1181    rt_list : list
1182        List of retention times.
1183    tic_list : list
1184        List of total ion currents.
1185    """
1186
1187    def __init__(
1188        self, chromatogram_parent, mass_spectra, apex_index, rt_list, tic_list
1189    ):
1190        self._ri = None
1191        self._rt_list = list(rt_list)
1192        self._tic_list = list(tic_list)
1193        self.mass_spectra = list(mass_spectra)
1194        super().__init__(
1195            chromatogram_parent,
1196            self.mass_spectra[apex_index],
1197            (0, apex_index, len(self.mass_spectra) - 1),
1198        )
1199
1200    @property
1201    def rt_list(self):
1202        """Get the list of retention times.
1203
1204        Returns
1205        -------
1206        list
1207            The list of retention times.
1208        """
1209        return self._rt_list
1210
1211    @property
1212    def tic_list(self):
1213        """Get the list of total ion currents.
1214
1215        Returns
1216        -------
1217        list
1218            The list of total ion currents.
1219        """
1220        return self._tic_list