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