corems.mass_spectra.input.corems_hdf5

   1__author__ = "Yuri E. Corilo"
   2__date__ = "Oct 29, 2019"
   3
   4
   5from threading import Thread
   6import h5py
   7import tomlkit
   8import json
   9import multiprocessing
  10from pathlib import Path
  11from typing import Union, Tuple, List, Optional
  12
  13import numpy as np
  14import pandas as pd
  15import warnings
  16
  17from corems.chroma_peak.factory.chroma_peak_classes import LCMSMassFeature
  18from corems.encapsulation.spawn_pool import spawn_pool
  19from corems.encapsulation.input.parameter_from_json import (
  20    load_and_set_json_parameters_lcms,
  21    load_and_set_toml_parameters_lcms,
  22)
  23from corems.mass_spectra.factory.lc_class import (
  24    LCMSBase,
  25    MassSpectraBase,
  26    LCMSCollection,
  27)
  28from corems.mass_spectra.factory.chromat_data import EIC_Data
  29from corems.mass_spectra.input.parserbase import SpectraParserInterface
  30from corems.mass_spectrum.input.coremsHDF5 import ReadCoreMSHDF_MassSpectrum
  31from corems.molecular_id.factory.spectrum_search_results import SpectrumSearchResults
  32from corems.mass_spectra.input.mzml import MZMLSpectraParser
  33
  34
  35def create_manifest_from_folder(
  36    folder_path: Path,
  37    output_path: Path = None,
  38    batch_time_threshold_hours: float = 12.0,
  39    center_name: str = None,
  40    overwrite: bool = False,
  41) -> Path:
  42    """
  43    Create a manifest CSV file for ReadCoreMSHDFMassSpectraCollection from CoreMS HDF5 files.
  44
  45    Scans a folder for .corems subdirectories and generates a manifest with columns:
  46    sample_name, batch, order, center, time. Files are batched by creation time, and
  47    one sample is designated as the retention time alignment center.
  48
  49    Parameters
  50    ----------
  51    folder_path : Path
  52        Path to folder containing .corems subdirectories with HDF5 files.
  53    output_path : Path, optional
  54        Output manifest CSV path. Default: folder_path/manifest.csv.
  55    batch_time_threshold_hours : float, optional
  56        Time gap in hours for batch separation. Default: 12.0.
  57    center_name : str, optional
  58        Sample name to designate as RT alignment center (must exist in samples).
  59        If None, the middle sample (by creation time) is used.
  60    overwrite : bool, optional
  61        Whether to overwrite existing manifest. Default: False.
  62
  63    Returns
  64    -------
  65    Path
  66        Path to created manifest file.
  67
  68    Raises
  69    ------
  70    FileNotFoundError
  71        If folder_path doesn't exist or contains no .corems subdirectories.
  72    FileExistsError
  73        If output file exists and overwrite is False.
  74    ValueError
  75        If no HDF5 files found, or center_name doesn't match any sample.
  76    """
  77    if not folder_path.exists():
  78        raise FileNotFoundError(f"Folder {folder_path} does not exist.")
  79
  80    # Set default output path if not provided
  81    if output_path is None:
  82        output_path = folder_path / "manifest.csv"
  83
  84    # Check if output file exists
  85    if output_path.exists() and not overwrite:
  86        raise FileExistsError(
  87            f"Manifest file {output_path} already exists. "
  88            "Set overwrite=True to replace it."
  89        )
  90
  91    # Find all .corems subdirectories
  92    corems_dirs = sorted(
  93        [d for d in folder_path.iterdir() if d.is_dir() and d.suffix == ".corems"]
  94    )
  95
  96    if not corems_dirs:
  97        raise FileNotFoundError(
  98            f"No .corems subdirectories found in {folder_path}. "
  99            "Ensure the folder contains processed CoreMS data."
 100        )
 101
 102    # Collect sample information
 103    sample_data = []
 104
 105    for corems_dir in corems_dirs:
 106        sample_name = corems_dir.stem  # Remove .corems extension
 107        hdf5_file = corems_dir / f"{sample_name}.hdf5"
 108
 109        if not hdf5_file.exists():
 110            print(f"Warning: HDF5 file not found for {sample_name}, skipping.")
 111            continue
 112
 113        # Get creation time using the ReadCoreMSHDFMassSpectra method
 114        try:
 115            # Use context manager to ensure file is properly closed
 116            with ReadCoreMSHDFMassSpectra(str(hdf5_file)) as parser:
 117                # Use the get_original_creation_time() method which checks HDF5 attrs first,
 118                # then falls back to original parser if needed
 119                creation_time = parser.get_original_creation_time()
 120
 121                # Skip sample if creation time unavailable
 122                if creation_time is None:
 123                    print(
 124                        f"Warning: Could not get original creation time for {sample_name}, skipping."
 125                    )
 126                    continue
 127
 128        except Exception as e:
 129            print(
 130                f"Warning: Error getting creation time for {sample_name}: {e}, skipping."
 131            )
 132            continue
 133
 134        sample_data.append(
 135            {
 136                "sample_name": sample_name,
 137                "creation_time": creation_time,
 138                "hdf5_path": hdf5_file,
 139            }
 140        )
 141
 142    if not sample_data:
 143        raise ValueError(
 144            f"No valid HDF5 files found in {folder_path}. "
 145            "Ensure .corems subdirectories contain .hdf5 files."
 146        )
 147
 148    # Sort by creation time
 149    sample_data.sort(key=lambda x: x["creation_time"])
 150
 151    # Assign batches based on time threshold
 152    batch_assignments = []
 153    current_batch = 1
 154
 155    for i, sample in enumerate(sample_data):
 156        if i == 0:
 157            batch_assignments.append(current_batch)
 158        else:
 159            time_diff = sample["creation_time"] - sample_data[i - 1]["creation_time"]
 160            time_diff_hours = time_diff.total_seconds() / 3600
 161
 162            if time_diff_hours > batch_time_threshold_hours:
 163                current_batch += 1
 164
 165            batch_assignments.append(current_batch)
 166
 167    # Determine which sample should be the center for retention time alignment
 168    sample_names = [s["sample_name"] for s in sample_data]
 169
 170    if center_name is not None:
 171        # Validate that center_name is in the discovered samples
 172        if center_name not in sample_names:
 173            raise ValueError(
 174                f"Specified center_name '{center_name}' not found in discovered samples. "
 175                f"Available samples: {', '.join(sample_names)}"
 176            )
 177        center_sample = center_name
 178    else:
 179        # Use the middle sample (by creation time) as center
 180        middle_idx = len(sample_data) // 2
 181        center_sample = sample_data[middle_idx]["sample_name"]
 182        print(
 183            f"Auto-selected center sample: {center_sample} (index {middle_idx} of {len(sample_data)}, middle by creation time)"
 184        )
 185
 186    # Create manifest dataframe with center column as TRUE/FALSE
 187    manifest_df = pd.DataFrame(
 188        {
 189            "sample_name": sample_names,
 190            "batch": batch_assignments,
 191            "order": list(range(1, len(sample_data) + 1)),
 192            "center": [
 193                "TRUE" if name == center_sample else "FALSE" for name in sample_names
 194            ],
 195            "time": [
 196                s["creation_time"].strftime("%Y-%m-%dT%H:%M:%SZ") for s in sample_data
 197            ],
 198        }
 199    )
 200
 201    # Sort manifest by time before saving to ensure proper order
 202    manifest_df = manifest_df.sort_values("time").reset_index(drop=True)
 203    # Update order column to reflect sorted order
 204    manifest_df["order"] = list(range(1, len(manifest_df) + 1))
 205
 206    # Save manifest
 207    manifest_df.to_csv(output_path, index=False)
 208
 209    print(f"Manifest created successfully at {output_path}")
 210    print(f"Total samples: {len(sample_data)}")
 211    print(f"Number of batches: {current_batch}")
 212    print(
 213        f"Batch assignments: {dict(zip(range(1, current_batch + 1), [batch_assignments.count(b) for b in range(1, current_batch + 1)]))}"
 214    )
 215
 216    return output_path
 217
 218
 219class ReadCoreMSHDFMassSpectra(
 220    SpectraParserInterface, ReadCoreMSHDF_MassSpectrum, Thread
 221):
 222    """Class to read CoreMS HDF5 files and populate a LCMS or MassSpectraBase object.
 223
 224    Parameters
 225    ----------
 226    file_location : str
 227        The location of the HDF5 file to read, including the suffix.
 228
 229    Attributes
 230    ----------
 231    file_location : str
 232        The location of the HDF5 file to read.
 233    h5pydata : h5py.File
 234        The HDF5 file object.
 235    scans : list
 236        A list of the location of individual mass spectra within the HDF5 file.
 237    scan_number_list : list
 238        A list of the scan numbers of the mass spectra within the HDF5 file.
 239    parameters_location : str
 240        The location of the parameters file (json or toml).
 241
 242    Methods
 243    -------
 244    * import_mass_spectra(mass_spectra).
 245        Imports all mass spectra from the HDF5 file onto the LCMS or MassSpectraBase object.
 246    * get_mass_spectrum_from_scan(scan_number).
 247        Return mass spectrum data object from scan number.
 248    * load().
 249        Placeholder method to meet the requirements of the SpectraParserInterface.
 250    * run(mass_spectra).
 251        Runs the importer functions to populate a LCMS or MassSpectraBase object.
 252    * import_scan_info(mass_spectra).
 253        Imports the scan info from the HDF5 file to populate the _scan_info attribute
 254        on the LCMS or MassSpectraBase object
 255    * import_ms_unprocessed(mass_spectra).
 256        Imports the unprocessed mass spectra from the HDF5 file to populate the
 257        _ms_unprocessed attribute on the LCMS or MassSpectraBase object
 258    * import_parameters(mass_spectra).
 259        Imports the parameters from the HDF5 file to populate the parameters
 260        attribute on the LCMS or MassSpectraBase object
 261    * import_mass_features(mass_spectra).
 262        Imports the mass features from the HDF5 file to populate the mass_features
 263        attribute on the LCMS or MassSpectraBase object
 264    * import_eics(mass_spectra).
 265        Imports the extracted ion chromatograms from the HDF5 file to populate the
 266        eics attribute on the LCMS or MassSpectraBase object
 267    * import_spectral_search_results(mass_spectra).
 268        Imports the spectral search results from the HDF5 file to populate the
 269        spectral_search_results attribute on the LCMS or MassSpectraBase object
 270    * get_mass_spectra_obj().
 271        Return mass spectra data object, populating the _ms list on the LCMS or
 272        MassSpectraBase object from the HDF5 file
 273    * get_lcms_obj().
 274        Return LCMSBase object, populating the majority of the attributes on the
 275        LCMS object from the HDF5 file
 276
 277    """
 278
 279    def __init__(self, file_location: str):
 280        Thread.__init__(self)
 281        ReadCoreMSHDF_MassSpectrum.__init__(self, file_location)
 282
 283        # override the scans attribute on ReadCoreMSHDF_MassSpectrum class to expect a nested location within the HDF5 file
 284        self.scans = [
 285            "mass_spectra/" + x for x in list(self.h5pydata["mass_spectra"].keys())
 286        ]
 287        self.scan_number_list = sorted(
 288            [int(float(i)) for i in list(self.h5pydata["mass_spectra"].keys())]
 289        )
 290
 291        # set the location of the parameters file (json or toml)
 292        add_files = [
 293            x
 294            for x in self.file_location.parent.glob(
 295                self.file_location.name.replace(".hdf5", ".*")
 296            )
 297            if x.suffix != ".hdf5"
 298        ]
 299        if len([x for x in add_files if x.suffix == ".json"]) > 0:
 300            self.parameters_location = [x for x in add_files if x.suffix == ".json"][0]
 301        elif len([x for x in add_files if x.suffix == ".toml"]) > 0:
 302            self.parameters_location = [x for x in add_files if x.suffix == ".toml"][0]
 303        else:
 304            self.parameters_location = None
 305
 306    def __enter__(self):
 307        """Context manager entry."""
 308        return self
 309
 310    def __exit__(self, exc_type, exc_val, exc_tb):
 311        """Context manager exit - closes the HDF5 file."""
 312        if hasattr(self, "h5pydata") and self.h5pydata is not None:
 313            self.h5pydata.close()
 314        return False
 315
 316    def close(self):
 317        """Explicitly close the HDF5 file."""
 318        if hasattr(self, "h5pydata") and self.h5pydata is not None:
 319            self.h5pydata.close()
 320
 321    def get_mass_spectrum_from_scan(self, scan_number):
 322        """Return mass spectrum data object from scan number."""
 323        if scan_number in self.scan_number_list:
 324            mass_spec = self.get_mass_spectrum(scan_number)
 325            return mass_spec
 326        else:
 327            raise Exception("Scan number not found in HDF5 file.")
 328
 329    def get_mass_spectra_from_scan_list(
 330        self, scan_list, spectrum_mode, auto_process=True
 331    ):
 332        """Return a list of mass spectrum data objects from a list of scan numbers.
 333
 334        Parameters
 335        ----------
 336        scan_list : list
 337            A list of scan numbers to retrieve mass spectra for.
 338        spectrum_mode : str
 339            The spectrum mode to use when retrieving the mass spectra.
 340            Note that this parameter is not used for CoreMS HDF5 files, as the spectra are already processed and only
 341            centroided spectra are saved.
 342        auto_process : bool
 343            If True, automatically process the mass spectra when retrieving them.
 344            Note that this parameter is not used for CoreMS HDF5 files, as the spectra are already processed and only
 345            centroided spectra are saved.
 346
 347        Returns
 348        -------
 349        list
 350            A list of mass spectrum data objects corresponding to the provided scan numbers.
 351        """
 352        mass_spectra_list = []
 353        for scan_number in scan_list:
 354            if scan_number in self.scan_number_list:
 355                mass_spec = self.get_mass_spectrum_from_scan(scan_number)
 356                mass_spectra_list.append(mass_spec)
 357            else:
 358                warnings.warn(f"Scan number {scan_number} not found in HDF5 file.")
 359        return mass_spectra_list
 360
 361    def load(self) -> None:
 362        """ """
 363        pass
 364
 365    def get_ms_raw(
 366        self,
 367        spectra=None,
 368        scan_df=None,
 369        time_range: Optional[
 370            Union[Tuple[float, float], List[Tuple[float, float]]]
 371        ] = None,
 372    ) -> dict:
 373        """ """
 374        # Warn if spectra or scan_df are not None that they are not used for CoreMS HDF5 files and should be rerun after instantiation
 375        if spectra is not None or scan_df is not None:
 376            SyntaxWarning(
 377                "get_ms_raw method for CoreMS HDF5 files can only access saved data, consider rerunning after instantiation."
 378            )
 379        ms_unprocessed = {}
 380        dict_group_load = self.h5pydata["ms_unprocessed"]
 381        dict_group_keys = dict_group_load.keys()
 382        for k in dict_group_keys:
 383            dataset = dict_group_load[k]
 384            ms_up_int = dataset[:]
 385            column_labels_json = dataset.attrs.get("column_labels")
 386            column_labels = (
 387                json.loads(column_labels_json)
 388                if column_labels_json
 389                else ["scan", "mz", "intensity"]
 390            )
 391            ms_unprocessed[int(k)] = pd.DataFrame(ms_up_int, columns=column_labels)
 392        return ms_unprocessed
 393
 394    def get_scan_df(
 395        self,
 396        time_range: Optional[
 397            Union[Tuple[float, float], List[Tuple[float, float]]]
 398        ] = None,
 399    ) -> pd.DataFrame:
 400        scan_info = {}
 401        dict_group_load = self.h5pydata["scan_info"]
 402        dict_group_keys = dict_group_load.keys()
 403        for k in dict_group_keys:
 404            scan_info[k] = dict_group_load[k][:]
 405        scan_df = pd.DataFrame(scan_info)
 406        scan_df.set_index("scan", inplace=True, drop=False)
 407        # HDF5 stores text as fixed-width bytes (numpy "S"). Decode only those
 408        # columns. Do not use select_dtypes([object]) + str.decode: under pandas
 409        # 3 / infer_string, already-decoded string columns can be selected and
 410        # corrupted by a blanket decode.
 411        for col in scan_df.columns:
 412            series = scan_df[col]
 413            sample = series.dropna()
 414            if sample.empty:
 415                continue
 416            first = sample.iloc[0]
 417            if isinstance(first, (bytes, bytearray, np.bytes_)):
 418                scan_df[col] = series.map(
 419                    lambda x: (
 420                        x.decode("utf-8")
 421                        if isinstance(x, (bytes, bytearray, np.bytes_))
 422                        else x
 423                    )
 424                )
 425
 426        # Apply time range filtering if specified
 427        if time_range is not None:
 428            time_ranges = self._normalize_time_range(time_range)
 429            mask = pd.Series([False] * len(scan_df), index=scan_df.index)
 430            for start_time, end_time in time_ranges:
 431                mask |= (scan_df.scan_time >= start_time) & (
 432                    scan_df.scan_time <= end_time
 433                )
 434            scan_df = scan_df[mask]
 435
 436        return scan_df
 437
 438    def run(self, mass_spectra, load_raw=True, load_light=False) -> None:
 439        """Runs the importer functions to populate a LCMS or MassSpectraBase object.
 440
 441        Notes
 442        -----
 443        The following functions are run in order, if the HDF5 file contains the necessary data:
 444        1. import_parameters(), which populates the parameters attribute on the LCMS or MassSpectraBase object.
 445        2. import_mass_spectra(), which populates the _ms list on the LCMS or MassSpectraBase object.
 446        3. import_scan_info(), which populates the _scan_info on the LCMS or MassSpectraBase object.
 447        4. import_ms_unprocessed(), which populates the _ms_unprocessed attribute on the LCMS or MassSpectraBase object.
 448        5. import_mass_features(), which populates the mass_features attribute on the LCMS or MassSpectraBase object.
 449        6. import_eics(), which populates the eics attribute on the LCMS or MassSpectraBase object.
 450        7. import_spectral_search_results(), which populates the spectral_search_results attribute on the LCMS or MassSpectraBase object.
 451
 452        Parameters
 453        ----------
 454        mass_spectra : LCMSBase or MassSpectraBase
 455            The LCMS or MassSpectraBase object to populate with mass spectra, generally instantiated with only the file_location, analyzer, and instrument_label attributes.
 456        load_raw : bool
 457            If True, load raw data (unprocessed) from HDF5 files for overall lcms object and individual mass spectra. Default is True.
 458        load_light : bool
 459            If True, only load the parameters, mass features, and scan info. Default is False.
 460
 461        Returns
 462        -------
 463        None, but populates several attributes on the LCMS or MassSpectraBase object.
 464
 465        """
 466        if self.parameters_location is not None:
 467            # Populate the parameters attribute on the LCMS object
 468            self.import_parameters(mass_spectra)
 469
 470        if "mass_spectra" in self.h5pydata and not load_light:
 471            # Populate the _ms list on the LCMS object
 472            self.import_mass_spectra(mass_spectra, load_raw=load_raw)
 473
 474        if "scan_info" in self.h5pydata:
 475            # Populate the _scan_info attribute on the LCMS object
 476            self.import_scan_info(mass_spectra)
 477
 478        if "ms_unprocessed" in self.h5pydata and load_raw and not load_light:
 479            # Populate the _ms_unprocessed attribute on the LCMS object
 480            self.import_ms_unprocessed(mass_spectra)
 481
 482        if "mass_features" in self.h5pydata:
 483            # Populate the mass_features attribute on the LCMS object
 484            self.import_mass_features(mass_spectra)
 485
 486        if "eics" in self.h5pydata and not load_light:
 487            # Populate the eics attribute on the LCMS object
 488            self.import_eics(mass_spectra)
 489
 490        if "spectral_search_results" in self.h5pydata and not load_light:
 491            # Populate the spectral_search_results attribute on the LCMS object
 492            self.import_spectral_search_results(mass_spectra)
 493
 494    def import_mass_spectra(self, mass_spectra, load_raw=True) -> None:
 495        """Imports all mass spectra from the HDF5 file.
 496
 497        Parameters
 498        ----------
 499        mass_spectra : LCMSBase | MassSpectraBase
 500            The MassSpectraBase or LCMSBase object to populate with mass spectra.
 501        load_raw : bool
 502            If True, load raw data (unprocessed) from HDF5 files for overall lcms object and individual mass spectra. Default
 503
 504        Returns
 505        -------
 506        None, but populates the '_ms' list on the LCMSBase or MassSpectraBase
 507        object with mass spectra from the HDF5 file.
 508        """
 509        for scan_number in self.scan_number_list:
 510            mass_spec = self.get_mass_spectrum(scan_number, load_raw=load_raw)
 511            mass_spec.scan_number = scan_number
 512            mass_spectra.add_mass_spectrum(mass_spec)
 513
 514    def import_scan_info(self, mass_spectra) -> None:
 515        """Imports the scan info from the HDF5 file.
 516
 517        Parameters
 518        ----------
 519        lcms : LCMSBase | MassSpectraBase
 520            The MassSpectraBase or LCMSBase objects
 521
 522        Returns
 523        -------
 524        None, but populates the 'scan_df' attribute on the LCMSBase or MassSpectraBase
 525        object with a pandas DataFrame of the 'scan_info' from the HDF5 file.
 526
 527        """
 528        scan_df = self.get_scan_df()
 529        mass_spectra.scan_df = scan_df
 530
 531    def import_ms_unprocessed(self, mass_spectra) -> None:
 532        """Imports the unprocessed mass spectra from the HDF5 file.
 533
 534        Parameters
 535        ----------
 536        lcms : LCMSBase | MassSpectraBase
 537            The MassSpectraBase or LCMSBase objects
 538
 539        Returns
 540        -------
 541        None, but populates the '_ms_unprocessed' attribute on the LCMSBase or MassSpectraBase
 542        object with a dictionary of the 'ms_unprocessed' from the HDF5 file.
 543
 544        """
 545        ms_unprocessed = self.get_ms_raw()
 546        mass_spectra._ms_unprocessed = ms_unprocessed
 547
 548    def import_parameters(self, mass_spectra) -> None:
 549        """Imports the parameters from the HDF5 file.
 550
 551        Parameters
 552        ----------
 553        mass_spectra : LCMSBase | MassSpectraBase
 554            The MassSpectraBase or LCMSBase object to populate with parameters.
 555
 556        Returns
 557        -------
 558        None, but populates the 'parameters' attribute on the LCMS or MassSpectraBase
 559        object with a dictionary of the 'parameters' from the HDF5 file.
 560
 561        """
 562        if ".json" == self.parameters_location.suffix:
 563            load_and_set_json_parameters_lcms(mass_spectra, self.parameters_location)
 564        if ".toml" == self.parameters_location.suffix:
 565            load_and_set_toml_parameters_lcms(mass_spectra, self.parameters_location)
 566        else:
 567            raise Exception(
 568                "Parameters file must be in JSON format, TOML format is not yet supported."
 569            )
 570
 571    def import_mass_features(self, mass_spectra, mf_ids=None) -> None:
 572        """Imports the mass features from the HDF5 file.
 573
 574        Parameters
 575        ----------
 576        mass_spectra : LCMSBase | MassSpectraBase
 577            The MassSpectraBase or LCMSBase object to populate with mass features.
 578        mf_ids : list, optional
 579            A list of mass feature IDs to import. If None, all mass features are imported.
 580
 581        Returns
 582        -------
 583        None, but populates the 'mass_features' attribute on the LCMSBase or MassSpectraBase
 584        object with a dictionary of the 'mass_features' from the HDF5 file.
 585
 586        """
 587        dict_group_load = self.h5pydata["mass_features"]
 588        dict_group_keys = dict_group_load.keys()
 589        for k in dict_group_keys:
 590            if mf_ids is not None and int(k) not in mf_ids:
 591                continue
 592            # Instantiate the MassFeature object
 593            mass_feature = LCMSMassFeature(
 594                mass_spectra,
 595                mz=dict_group_load[k].attrs["_mz_exp"],
 596                retention_time=dict_group_load[k].attrs["_retention_time"],
 597                intensity=dict_group_load[k].attrs["_intensity"],
 598                apex_scan=dict_group_load[k].attrs["_apex_scan"],
 599                persistence=dict_group_load[k].attrs["_persistence"],
 600                id=int(k),
 601            )
 602
 603            # Populate additional attributes on the MassFeature object
 604            for key in dict_group_load[k].attrs.keys() - {
 605                "_mz_exp",
 606                "_mz_cal",
 607                "_retention_time",
 608                "_intensity",
 609                "_apex_scan",
 610                "_persistence",
 611            }:
 612                setattr(mass_feature, key, dict_group_load[k].attrs[key])
 613
 614            # Populate attributes on MassFeature object that are lists
 615            for key in dict_group_load[k].keys():
 616                setattr(mass_feature, key, dict_group_load[k][key][:])
 617                # Convert _noise_score from array to tuple
 618                if key == "_noise_score":
 619                    mass_feature._noise_score = tuple(mass_feature._noise_score)
 620            mass_spectra.mass_features[int(k)] = mass_feature
 621
 622        # Associate mass features with ms1 and ms2 spectra, if available
 623        for mf_id in mass_spectra.mass_features.keys():
 624            if mass_spectra.mass_features[mf_id].apex_scan in mass_spectra._ms.keys():
 625                mass_spectra.mass_features[mf_id].mass_spectrum = mass_spectra._ms[
 626                    mass_spectra.mass_features[mf_id].apex_scan
 627                ]
 628            if mass_spectra.mass_features[mf_id].ms2_scan_numbers is not None:
 629                for ms2_scan in mass_spectra.mass_features[mf_id].ms2_scan_numbers:
 630                    if ms2_scan in mass_spectra._ms.keys():
 631                        mass_spectra.mass_features[mf_id].ms2_mass_spectra[ms2_scan] = (
 632                            mass_spectra._ms[ms2_scan]
 633                        )
 634
 635    def import_eics(self, mass_spectra, mz_list=None, mz_tolerance=0.0001):
 636        """Imports the extracted ion chromatograms from the HDF5 file.
 637
 638        Parameters
 639        ----------
 640        mass_spectra : LCMSBase | MassSpectraBase
 641            The MassSpectraBase or LCMSBase object to populate with extracted ion chromatograms.
 642        mz_list : list of float, optional
 643            List of m/z values to load EICs for. If None, loads all EICs. Default is None.
 644        mz_tolerance : float, optional
 645            Tolerance in Daltons for matching m/z values when mz_list is provided.
 646            Default is 0.0001 Da.
 647
 648        Returns
 649        -------
 650        None, but populates the 'eics' attribute on the LCMSBase or MassSpectraBase
 651        object with a dictionary of the 'eics' from the HDF5 file.
 652
 653        """
 654        if "eics" not in self.h5pydata:
 655            # Samples with no exported EICs (e.g. empty / fully gap-filled) have no group
 656            return
 657
 658        dict_group_load = self.h5pydata["eics"]
 659        dict_group_keys = dict_group_load.keys()
 660
 661        # Prefilter dict_group_keys if mz_list is provided to EICs within tolerance
 662        if mz_list is not None:
 663            if len(mz_list) == 0:
 664                return
 665            target_mz_array = np.array(sorted(mz_list), dtype=float)
 666            mzs = [
 667                float(k)
 668                for k in dict_group_keys
 669                if np.abs(float(k) - target_mz_array).min() < mz_tolerance
 670            ]
 671            dict_group_keys = [str(mz) for mz in mzs]
 672
 673        for k in dict_group_keys:
 674            # Check if we should load this EIC (filter by m/z if list provided)
 675            eic_mz = dict_group_load[k].attrs["mz"]
 676
 677            my_eic = EIC_Data(
 678                scans=dict_group_load[k]["scans"][:],
 679                time=dict_group_load[k]["time"][:],
 680                eic=dict_group_load[k]["eic"][:],
 681            )
 682            for key in dict_group_load[k].keys():
 683                if key not in ["scans", "time", "eic"]:
 684                    setattr(my_eic, key, dict_group_load[k][key][:])
 685                    # if key is apexes, convert to a tuple of a list
 686                    if key == "apexes" and len(my_eic.apexes) > 0:
 687                        my_eic.apexes = [tuple(x) for x in my_eic.apexes]
 688            # Add to mass_spectra object
 689            mass_spectra.eics[eic_mz] = my_eic
 690
 691        # Associate EICs with mass features using tolerance-based matching
 692        mass_spectra.associate_eics_with_mass_features()
 693
 694    @staticmethod
 695    def _load_eics_from_hdf5_group(eics_group, lcms_obj, mz_filter=None):
 696        """Load EICs from an HDF5 group.
 697
 698        This is a static helper method that can be reused to load EIC data
 699        from any HDF5 group in a consistent format.
 700
 701        Parameters
 702        ----------
 703        eics_group : h5py.Group
 704            The HDF5 group containing EIC data.
 705        lcms_obj : LCMSBase
 706            The LCMS object to associate EICs with (for reference, not modified).
 707        mz_filter : list, optional
 708            List of m/z values to load. If None, loads all EICs. Default is None.
 709            Uses tolerance-based matching (0.0001).
 710
 711        Returns
 712        -------
 713        dict
 714            Dictionary of EIC_Data objects keyed by m/z value.
 715        """
 716        from corems.mass_spectra.factory.chromat_data import EIC_Data
 717
 718        loaded_eics = {}
 719        tolerance = 0.0001
 720
 721        for eic_key_str in eics_group.keys():
 722            eic_mz = (
 723                float(eic_key_str)
 724                if eic_key_str.replace(".", "", 1).replace("-", "", 1).isdigit()
 725                else eics_group[eic_key_str].attrs.get("mz")
 726            )
 727
 728            # If mz_filter is provided, check if this EIC matches any requested m/z
 729            if mz_filter is not None:
 730                if not any(abs(eic_mz - mz) < tolerance for mz in mz_filter):
 731                    continue
 732
 733            eic_data = eics_group[eic_key_str]
 734
 735            # Create EIC_Data object from datasets
 736            eic = EIC_Data(
 737                scans=list(eic_data["scans"][:]) if "scans" in eic_data else [],
 738                time=list(eic_data["time"][:]) if "time" in eic_data else [],
 739                eic=list(eic_data["eic"][:]) if "eic" in eic_data else [],
 740                apexes=list(eic_data["apexes"][:]) if "apexes" in eic_data else [],
 741            )
 742
 743            # Load any additional datasets
 744            for key in eic_data.keys():
 745                if key not in ["scans", "time", "eic", "apexes"]:
 746                    setattr(eic, key, eic_data[key][:])
 747
 748            loaded_eics[eic_mz] = eic
 749
 750        return loaded_eics
 751
 752    def import_spectral_search_results(self, mass_spectra):
 753        """Imports the spectral search results from the HDF5 file.
 754
 755        Parameters
 756        ----------
 757        mass_spectra : LCMSBase | MassSpectraBase
 758            The MassSpectraBase or LCMSBase object to populate with spectral search results.
 759
 760        Returns
 761        -------
 762        None, but populates the 'spectral_search_results' attribute on the LCMSBase or MassSpectraBase
 763        object with a dictionary of the 'spectral_search_results' from the HDF5 file.
 764
 765        """
 766        overall_results_dict = {}
 767        ms2_results_load = self.h5pydata["spectral_search_results"]
 768        for k in ms2_results_load.keys():
 769            overall_results_dict[int(k)] = {}
 770            for k2 in ms2_results_load[k].keys():
 771                ms2_search_res = SpectrumSearchResults(
 772                    query_spectrum=mass_spectra._ms[int(k)],
 773                    precursor_mz=ms2_results_load[k][k2].attrs["precursor_mz"],
 774                    spectral_similarity_search_results={},
 775                )
 776
 777                for key in ms2_results_load[k][k2].keys() - {"precursor_mz"}:
 778                    data = list(ms2_results_load[k][k2][key][:])
 779                    if data and isinstance(data[0], bytes):
 780                        data = [x.decode("utf-8") for x in data]
 781                    setattr(ms2_search_res, key, data)
 782
 783                overall_results_dict[int(k)][
 784                    ms2_results_load[k][k2].attrs["precursor_mz"]
 785                ] = ms2_search_res
 786
 787        # add to mass_spectra
 788        mass_spectra.spectral_search_results.update(overall_results_dict)
 789
 790        # If there are mass features, associate the results with each mass feature
 791        if len(mass_spectra.mass_features) > 0:
 792            for mass_feature_id, mass_feature in mass_spectra.mass_features.items():
 793                scan_ids = mass_feature.ms2_scan_numbers
 794                for ms2_scan_id in scan_ids:
 795                    precursor_mz = mass_feature.mz
 796                    try:
 797                        mass_spectra.spectral_search_results[ms2_scan_id][precursor_mz]
 798                    except KeyError:
 799                        pass
 800                    else:
 801                        mass_spectra.mass_features[
 802                            mass_feature_id
 803                        ].ms2_similarity_results.append(
 804                            mass_spectra.spectral_search_results[ms2_scan_id][
 805                                precursor_mz
 806                            ]
 807                        )
 808
 809    def get_mass_spectra_obj(
 810        self,
 811        load_raw=True,
 812        load_light=False,
 813        time_range: Optional[
 814            Union[Tuple[float, float], List[Tuple[float, float]]]
 815        ] = None,
 816    ) -> MassSpectraBase:
 817        """
 818        Return mass spectra data object, populating the _ms list on MassSpectraBase object from the HDF5 file.
 819
 820        Parameters
 821        ----------
 822        load_raw : bool
 823            If True, load raw data (unprocessed) from HDF5 files for overall spectra object and individual mass spectra. Default is True.
 824        load_light : bool
 825            If True, only load the parameters, mass features, and scan info. Default is False.
 826        time_range : tuple or list of tuples, optional
 827            Retention time range(s) to load. Can be:
 828            - Single range: (start_time, end_time) in minutes
 829            - Multiple ranges: [(start1, end1), (start2, end2), ...] in minutes
 830            If None, loads all scans. Note: For HDF5 files, this parameter is accepted for
 831            interface consistency but not currently used in filtering.
 832
 833        """
 834        # Instantiate the LCMS object
 835        spectra_obj = MassSpectraBase(
 836            file_location=self.file_location,
 837            analyzer=self.analyzer,
 838            instrument_label=self.instrument_label,
 839            sample_name=self.sample_name,
 840        )
 841
 842        # This will populate the _ms list on the LCMS or MassSpectraBase object
 843        self.run(spectra_obj, load_raw=load_raw, load_light=load_light)
 844
 845        return spectra_obj
 846
 847    def get_lcms_obj(
 848        self,
 849        load_raw=True,
 850        load_light=False,
 851        use_original_parser=True,
 852        raw_file_path=None,
 853        time_range: Optional[
 854            Union[Tuple[float, float], List[Tuple[float, float]]]
 855        ] = None,
 856    ) -> LCMSBase:
 857        """
 858        Return LCMSBase object, populating attributes on the LCMSBase object from the HDF5 file.
 859
 860        Parameters
 861        ----------
 862        load_raw : bool
 863            If True, load raw data (unprocessed) from HDF5 files for overall lcms object and individual mass spectra. Default is True.
 864        load_light : bool
 865            If True, only load the parameters, mass features, and scan info. Default is False.
 866        use_original_parser : bool
 867            If True, use the original parser to populate the LCMS object. Default is True.
 868        raw_file_path : str
 869            The location of the raw file to parse if attempting to use original parser.
 870            Default is None, which attempts to get the raw file path from the HDF5 file.
 871            If the original file path has moved, this parameter can be used to specify the new location.
 872        time_range : tuple or list of tuples, optional
 873            Retention time range(s) to load. Can be:
 874            - Single range: (start_time, end_time) in minutes
 875            - Multiple ranges: [(start1, end1), (start2, end2), ...] in minutes
 876            If None, loads all scans. Note: For HDF5 files, this parameter is accepted for
 877            interface consistency. If use_original_parser=True, time_range can be passed to the
 878            original parser for filtering.
 879        """
 880        # Instantiate the LCMS object
 881        lcms_obj = LCMSBase(
 882            file_location=self.file_location,
 883            analyzer=self.analyzer,
 884            instrument_label=self.instrument_label,
 885            sample_name=self.sample_name,
 886        )
 887
 888        # This will populate the majority of the attributes on the LCMS object
 889        self.run(lcms_obj, load_raw=load_raw, load_light=load_light)
 890
 891        # Set final attributes of the LCMS object
 892        lcms_obj.polarity = self.h5pydata.attrs["polarity"]
 893        lcms_obj._scans_number_list = list(lcms_obj.scan_df.scan)
 894        lcms_obj._retention_time_list = list(lcms_obj.scan_df.scan_time)
 895        lcms_obj._tic_list = list(lcms_obj.scan_df.tic)
 896
 897        # If use_original_parser is True, instantiate the original parser and populate the LCMS object
 898        if use_original_parser:
 899            lcms_obj = self.add_original_parser(lcms_obj, raw_file_path=raw_file_path)
 900        else:
 901            lcms_obj.spectra_parser_class = self.__class__
 902
 903        return lcms_obj
 904
 905    def get_raw_file_location(self):
 906        """
 907        Get the raw file location from the HDF5 file attributes.
 908
 909        Returns
 910        -------
 911        str
 912            The raw file location.
 913        """
 914        if "original_file_location" in self.h5pydata.attrs:
 915            return self.h5pydata.attrs["original_file_location"]
 916        else:
 917            return None
 918
 919    def add_original_parser(self, mass_spectra, raw_file_path=None):
 920        """
 921        Add the original parser to the mass spectra object.
 922
 923        Parameters
 924        ----------
 925        mass_spectra : MassSpectraBase | LCMSBase
 926            The MassSpectraBase or LCMSBase object to add the original parser to.
 927        raw_file_path : str
 928            The location of the raw file to parse. Default is None, which attempts to get the raw file path from the HDF5 file.
 929        """
 930        # Get the original parser type
 931        og_parser_type = self.h5pydata.attrs["parser_type"]
 932
 933        # If raw_file_path is None, get it from the HDF5 file attributes
 934        if raw_file_path is None:
 935            raw_file_path = self.get_raw_file_location()
 936            if raw_file_path is None:
 937                raise ValueError(
 938                    "Raw file path not found in HDF5 file attributes, cannot instantiate original parser."
 939                )
 940
 941        # Set the raw file path on the mass_spectra object so the parser knows where to find the raw file
 942        mass_spectra.raw_file_location = raw_file_path
 943
 944        if og_parser_type == "ImportMassSpectraThermoMSFileReader":
 945            # Imported lazily: rawFileReader pulls in pythonnet (`clr`), which
 946            # needs a .NET/Mono runtime. Only reopening a Thermo .raw file
 947            # needs it; reading an HDF5 file of any other origin must not.
 948            from corems.mass_spectra.input.rawFileReader import (
 949                ImportMassSpectraThermoMSFileReader,
 950            )
 951
 952            # Check that the parser can be instantiated with the raw file path
 953            parser_class = ImportMassSpectraThermoMSFileReader
 954        elif og_parser_type == "MZMLSpectraParser":
 955            # Check that the parser can be instantiated with the raw file path
 956            parser_class = MZMLSpectraParser
 957
 958        # Set the spectra parser class on the mass_spectra object so the spectra_parser property can be used with the original parser
 959        mass_spectra.spectra_parser_class = parser_class
 960
 961        return mass_spectra
 962
 963    def get_original_creation_time(self):
 964        """
 965        Get the creation time of the original raw data file.
 966
 967        First checks if creation_time is saved in the HDF5 file attributes.
 968        If not found, attempts to instantiate the original parser and get the creation time.
 969
 970        Returns
 971        -------
 972        datetime
 973            The creation time of the original raw data file, or None if not available.
 974        """
 975        # Check if creation_time is saved in HDF5 attributes
 976        if "creation_time" in self.h5pydata.attrs:
 977            from datetime import datetime
 978
 979            return datetime.fromisoformat(self.h5pydata.attrs["creation_time"])
 980
 981        # Fall back to using original parser to get creation time
 982        try:
 983            # Get the original parser type and raw file path
 984            og_parser_type = self.h5pydata.attrs.get("parser_type")
 985            raw_file_path = self.get_raw_file_location()
 986
 987            if og_parser_type is None or raw_file_path is None:
 988                warnings.warn(
 989                    "Cannot retrieve creation time: parser_type or original_file_location not found in HDF5 attributes."
 990                )
 991                return None
 992
 993            # Check if raw file exists
 994            from pathlib import Path
 995
 996            if not Path(raw_file_path).exists():
 997                warnings.warn(
 998                    f"Cannot retrieve creation time: original raw file not found at {raw_file_path}"
 999                )
1000                return None
1001
1002            # Instantiate the original parser
1003            if og_parser_type == "ImportMassSpectraThermoMSFileReader":
1004                # Lazy import; see add_original_parser for why.
1005                from corems.mass_spectra.input.rawFileReader import (
1006                    ImportMassSpectraThermoMSFileReader,
1007                )
1008
1009                parser = ImportMassSpectraThermoMSFileReader(raw_file_path)
1010            elif og_parser_type == "MZMLSpectraParser":
1011                parser = MZMLSpectraParser(raw_file_path)
1012            else:
1013                warnings.warn(
1014                    f"Unknown parser type: {og_parser_type}, cannot retrieve creation time."
1015                )
1016                return None
1017
1018            # Get creation time from parser
1019            return parser.get_creation_time()
1020
1021        except Exception as e:
1022            warnings.warn(f"Failed to retrieve creation time from original parser: {e}")
1023            return None
1024
1025    def get_creation_time(self):
1026        """
1027        Get the creation time of the original raw data file.
1028
1029        This is an alias for get_original_creation_time() for backward compatibility.
1030
1031        Returns
1032        -------
1033        datetime
1034            The creation time of the original raw data file, or None if not available.
1035        """
1036        return self.get_original_creation_time()
1037
1038    def get_instrument_info(self):
1039        """
1040        Raise a NotImplemented Warning, as instrument info is not available in CoreMS HDF5 files and returning None.
1041        """
1042        warnings.warn(
1043            "Instrument info is not available in CoreMS HDF5 files, returning None."
1044            "This should be accessed through the original parser.",
1045        )
1046        return None
1047
1048    def get_scans_in_time_range(
1049        self,
1050        time_range: Union[Tuple[float, float], List[Tuple[float, float]]],
1051        ms_level: Optional[int] = None,
1052    ) -> List[int]:
1053        """Return scan numbers within specified retention time range(s).
1054
1055        Parameters
1056        ----------
1057        time_range : tuple or list of tuples
1058            Retention time range(s) in minutes. Can be:
1059            - Single range: (start_time, end_time)
1060            - Multiple ranges: [(start1, end1), (start2, end2), ...]
1061        ms_level : int, optional
1062            If specified, only return scans of this MS level (e.g., 1 for MS1, 2 for MS2).
1063            If None, returns scans of all MS levels.
1064
1065        Returns
1066        -------
1067        list of int
1068            List of scan numbers within the specified time range(s) and MS level.
1069        """
1070        # Normalize time range to list of tuples
1071        time_ranges = self._normalize_time_range(time_range)
1072
1073        # Get all scan data
1074        scan_df = self.get_scan_df()
1075
1076        # Filter by time range
1077        mask = pd.Series([False] * len(scan_df), index=scan_df.index)
1078        for start_time, end_time in time_ranges:
1079            mask |= (scan_df.scan_time >= start_time) & (scan_df.scan_time <= end_time)
1080
1081        filtered_df = scan_df[mask]
1082
1083        # Filter by MS level if specified
1084        if ms_level is not None:
1085            filtered_df = filtered_df[filtered_df.ms_level == ms_level]
1086
1087        return filtered_df.scan.tolist()
1088
1089
1090class ReadCoreMSHDFMassSpectraCollection:
1091    """Read a collection of CoreMS HDF5 files and populate an LCMSCollection object.
1092
1093    Parameters
1094    ----------
1095    folder_location : Path
1096        Folder containing .corems subdirectories with HDF5 files.
1097    manifest_file : Path, optional
1098        Manifest CSV with columns: sample_name, order, batch, center, time.
1099        One sample must have center='TRUE' for RT alignment.
1100        If None, checks if auto-generated manifest_auto.csv exists in the folder. If not,
1101        auto-generates from folder contents. Default: None.
1102    cores : int, optional
1103        Number of cores for multiprocessing. Default: 1.
1104    auto_manifest_batch_threshold_hours : float, optional
1105        Time gap (hours) for auto-generated batch separation. Default: 12.0.
1106    auto_manifest_center_name : str, optional
1107        Sample name for RT alignment center when auto-generating.
1108        Must match a discovered sample. If None, uses middle sample. Default: None.
1109
1110    Attributes
1111    ----------
1112    folder_location : Path
1113        Folder containing CoreMS HDF5 files.
1114    manifest_filepath : Path
1115        Path to manifest file.
1116    manifest : dict
1117        Manifest data indexed by sample_name.
1118    """
1119
1120    def __init__(
1121        self,
1122        folder_location: Path,
1123        manifest_file: Path = None,
1124        cores: int = 1,
1125        auto_manifest_batch_threshold_hours: float = 12.0,
1126        auto_manifest_center_name: str = None,
1127    ):
1128        # Check for folder location
1129        folder_location = Path(folder_location)
1130        if not folder_location.exists():
1131            raise FileNotFoundError(f"Folder location {folder_location} not found.")
1132
1133        # Auto-generate manifest if not provided
1134        if manifest_file is None:
1135            # Check if manifest_auto.csv already exists
1136            auto_manifest_path = folder_location / "manifest_auto.csv"
1137            if auto_manifest_path.exists():
1138                print(
1139                    f"No manifest file provided. Using existing manifest_auto.csv from {folder_location}"
1140                )
1141                manifest_file = auto_manifest_path
1142            else:
1143                print(
1144                    f"No manifest file provided. Auto-generating manifest from {folder_location}"
1145                )
1146                manifest_file = create_manifest_from_folder(
1147                    folder_path=folder_location,
1148                    output_path=auto_manifest_path,
1149                    batch_time_threshold_hours=auto_manifest_batch_threshold_hours,
1150                    center_name=auto_manifest_center_name,
1151                    overwrite=True,
1152                )
1153        else:
1154            manifest_file = Path(manifest_file)
1155            if not manifest_file.exists():
1156                raise FileNotFoundError(f"Manifest file {manifest_file} not found.")
1157
1158            # Check if the manifest file is a CSV
1159            if manifest_file.suffix != ".csv":
1160                raise ValueError("Manifest file must be a CSV.")
1161
1162        self.folder_location = folder_location
1163        self._manifest_dict = None
1164        self._parse_manifest(manifest_file)
1165        self._validate_manifest()
1166        self._validate_parameters()
1167        self._validate_cores(cores)
1168
1169    def _validate_cores(self, cores):
1170        # Check if the cores parameter is an integer greater than 0 and less than the number of cores available
1171        if not isinstance(cores, int) or cores < 1:
1172            raise ValueError("Cores must be an integer greater than 0.")
1173        if cores > multiprocessing.cpu_count():
1174            raise ValueError(
1175                f"Cores must be less than or equal to the number of cores available ({multiprocessing.cpu_count()})."
1176            )
1177        self._cores = cores
1178
1179    def _parse_manifest(self, manifest_file):
1180        """Parse the manifest file and set the manifest dictionary."""
1181        self.manifest_filepath = manifest_file
1182        manifest = pd.read_csv(manifest_file)
1183        # Check if the following columns exisit in the manifest file
1184        if not all(
1185            col in manifest.columns for col in ["sample_name", "order", "batch"]
1186        ):
1187            raise ValueError(
1188                "Manifest file must contain the following columns: 'sample_name', 'order', 'batch'."
1189            )
1190        # Set index to the 'sample_name' column
1191        manifest.set_index("sample_name", inplace=True)
1192        self._manifest_dict = manifest.to_dict(orient="index")
1193
1194    def _validate_manifest(self):
1195        """Validate the manifest dictionary against the CoreMS folder location."""
1196        # Check if the folder location contains HDF5 files for each sample
1197        for sample_name in self._manifest_dict.keys():
1198            corems_dir = self.folder_location / f"{sample_name}.corems"
1199            if not corems_dir.exists():
1200                raise FileNotFoundError(f"CoreMS folder for {sample_name} not found.")
1201            hdf5_file = corems_dir / f"{sample_name}.hdf5"
1202            if not hdf5_file.exists():
1203                raise FileNotFoundError(f"HDF5 file for {sample_name} not found.")
1204
1205        # Check that at least one sample has center='TRUE' for retention time alignment
1206        center_values = [
1207            sample_data.get("center") for sample_data in self._manifest_dict.values()
1208        ]
1209        if not any(
1210            center_val == "TRUE" or center_val == True for center_val in center_values
1211        ):
1212            raise ValueError(
1213                "Manifest must contain at least one sample with center='TRUE' for retention time alignment. "
1214                "None of the samples in the manifest have center='TRUE'."
1215            )
1216
1217    def _validate_parameters(self):
1218        """Validate that the parameters used for all samples within a batch are the same."""
1219        # Check if parameters files are saved as JSON or TOML
1220        if self.parameters_files[0].suffix == ".json":
1221            importer = json
1222            suffix = ".json"
1223
1224        elif self.parameters_files[0].suffix == ".toml":
1225            importer = tomlkit
1226            suffix = ".toml"
1227
1228        manfiest_df = self.manifest_dataframe
1229
1230        # Split up samples by batch
1231        batches = manfiest_df["batch"].unique()
1232
1233        for batch in batches:
1234            samples = manfiest_df[manfiest_df["batch"] == batch].index
1235            # check if self.parameters_files end with .json or .toml
1236            batch_param_files = [
1237                self.folder_location / f"{sample_name}.corems/{sample_name}{suffix}"
1238                for sample_name in self._manifest_dict.keys()
1239                if sample_name in samples
1240            ]
1241            with open(
1242                batch_param_files[0],
1243                "r",
1244                encoding="utf8",
1245            ) as stream:
1246                first_parameters = importer.load(stream)
1247            for parameters_file in batch_param_files[1:]:
1248                with open(
1249                    parameters_file,
1250                    "r",
1251                    encoding="utf8",
1252                ) as stream:
1253                    parameters = importer.load(stream)
1254                if parameters != first_parameters:
1255                    raise ValueError(
1256                        f"Parameters files for samples in batch {batch} are not equal."
1257                    )
1258
1259    def get_lcms_obj(
1260        self,
1261        sample_name: str,
1262        load_raw=False,
1263        load_light=True,
1264        use_original_parser=True,
1265        raw_file_path=None,
1266    ) -> LCMSBase:
1267        """Return a LCMSBase object for a given sample name within the collection.
1268
1269        Parameters
1270        ----------
1271        sample_name : str
1272            The sample name to retrieve the LCMS object for.
1273        load_raw : bool
1274            If True, load raw data from HDF5 files. Default is False.
1275        load_light : bool
1276            If True, only load the parameters, mass features, and scan info are initially loaded for each lcms object. Default is True.
1277        """
1278        hdf5_file = self.folder_location / f"{sample_name}.corems/{sample_name}.hdf5"
1279        with ReadCoreMSHDFMassSpectra(hdf5_file) as parser:
1280            lcms_obj = parser.get_lcms_obj(
1281                load_raw=load_raw,
1282                load_light=load_light,
1283                use_original_parser=use_original_parser,
1284                raw_file_path=raw_file_path,
1285            )
1286            if load_light:
1287                # mass_features_to_df always includes _eic_mz (object value or mz fallback)
1288                mf_df = lcms_obj.mass_features_to_df()
1289                lcms_obj.mass_features = {}
1290                lcms_obj.light_mf_df = mf_df
1291        return lcms_obj
1292
1293    def get_lcms_collection(
1294        self, load_raw=False, load_light=True, use_original_parser=True
1295    ) -> LCMSCollection:
1296        """Return a LCMSCollection object
1297
1298        Parameters
1299        ----------
1300        load_raw : bool
1301            If True, load raw data from HDF5 files. Default is False.
1302        load_light : bool
1303            If True, only load the parameters, mass features, and scan info are initially loaded for each lcms object.
1304            After concatenating the mass_features, remove the mass_features attribute from the individual LCMS objects for memory efficiency. Default is True.
1305            Default is True.
1306        """
1307        # Instantiate the LCMSCollection object
1308        lcms_coll = LCMSCollection(
1309            collection_location=self.folder_location,
1310            manifest=self.manifest,
1311            collection_parser=self,
1312        )
1313
1314        # Set the number of cores on the LCMSCollection object from the ReadCoreMSHDFMassSpectraCollection object
1315        lcms_coll.parameters.lcms_collection.cores = self._cores
1316
1317        # Add LCMS objects to the collection
1318        samples = self._manifest_dict.keys()
1319
1320        # Initialize the LCMS object dictionary
1321        if self._cores > 1:
1322            if self._cores > len(samples):
1323                ncores = len(samples)
1324            else:
1325                ncores = self._cores
1326            # Create a pool of workers (one for each core or sample, whichever is smaller)
1327            args = [
1328                (sample, load_raw, load_light, use_original_parser)
1329                for sample in samples
1330            ]
1331            with spawn_pool(ncores) as pool:
1332                lcms_objs = pool.starmap(self.get_lcms_obj, args)
1333            for sample_name, lcms_obj in zip(samples, lcms_objs):
1334                lcms_coll._lcms[sample_name] = lcms_obj
1335
1336        elif self._cores == 1:
1337            # Load the LCMS objects sequentially
1338            for sample_name in samples:
1339                lcms_coll._lcms[sample_name] = self.get_lcms_obj(
1340                    sample_name,
1341                    load_raw=load_raw,
1342                    load_light=load_light,
1343                    use_original_parser=use_original_parser,
1344                )
1345
1346        else:
1347            raise ValueError(
1348                "Number of cores must be greater than 0 and set on the ReadCoreMSHDFMassSpectraCollection object."
1349            )
1350
1351        # Check that all LCMS objects have the same polarity
1352        if len(set([x.polarity for k, x in lcms_coll._lcms.items()])) != 1:
1353            raise ValueError("All samples must have the same polarity.")
1354
1355        # Set ids on the LCMS objects in the manifest
1356        i = 0
1357        for sample in lcms_coll.samples:
1358            lcms_coll._manifest_dict[sample]["collection_id"] = i
1359            i += 1
1360
1361        # Reorder the LCMS objects
1362        lcms_coll._reorder_lcms_objects()
1363
1364        # Collect the mass features from the LCMS objects and combine them into a single dataframe for the collection
1365        lcms_coll._combine_mass_features()
1366
1367        # If load_light, remove the mass_feature attribute from the individual LCMS objects
1368        if load_light:
1369            for sample_name in lcms_coll.samples:
1370                lcms_coll._lcms[sample_name].mass_features = {}
1371                # Remove the light_mf_df attribute from the individual LCMS objects
1372                del lcms_coll._lcms[sample_name].light_mf_df
1373
1374        return lcms_coll
1375
1376    @property
1377    def manifest(self):
1378        return self._manifest_dict
1379
1380    @property
1381    def manifest_dataframe(self):
1382        return pd.DataFrame(self._manifest_dict).T
1383
1384    @property
1385    def hdf5_files(self):
1386        return [
1387            self.folder_location / f"{sample_name}.corems/{sample_name}.hdf5"
1388            for sample_name in self._manifest_dict.keys()
1389        ]
1390
1391    @property
1392    def parameters_files(self):
1393        # Check if parameters files are saved as JSON or TOML
1394        json_files = [
1395            self.folder_location / f"{sample_name}.corems/{sample_name}.json"
1396            for sample_name in self._manifest_dict.keys()
1397        ]
1398        toml_files = [
1399            self.folder_location / f"{sample_name}.corems/{sample_name}.toml"
1400            for sample_name in self._manifest_dict.keys()
1401        ]
1402        if all([x.exists() for x in json_files]):
1403            return json_files
1404        elif all([x.exists() for x in toml_files]):
1405            return toml_files
1406        else:
1407            raise ValueError("Parameters files are not saved for all samples.")
1408
1409
1410class ReadSavedLCMSCollection(ReadCoreMSHDFMassSpectraCollection):
1411    """
1412    Subclass to read and re-instantiate a LCMSCollection from a saved HDF5 file.
1413
1414
1415    Parameters
1416    ----------
1417    collection_hdf5_path : str or Path
1418        Path to the saved LCMSCollection HDF5 file.
1419    cores : int, optional
1420        Number of cores for processing. Default is 1.
1421    """
1422
1423    def __init__(self, collection_hdf5_path: str, cores: int = 1):
1424        # Convert to Path objects
1425        self.collection_hdf5_path = Path(collection_hdf5_path)
1426
1427        # Validate the collection file exists
1428        if not self.collection_hdf5_path.exists():
1429            raise FileNotFoundError(
1430                f"Collection HDF5 file {self.collection_hdf5_path} not found."
1431            )
1432
1433        # Validate cores
1434        self._validate_cores(cores)
1435
1436        # Load metadata from saved collection
1437        self._load_collection_metadata()
1438
1439        if not self.folder_location.exists():
1440            raise FileNotFoundError(
1441                f"Folder location {self.folder_location} not found."
1442            )
1443
1444        # Load the mass spectra data
1445        self._validate_manifest()
1446
1447        # Set the parameters file location
1448        self.parameters_location = self._get_parameters_location()
1449
1450    def _get_parameters_location(self):
1451        """Find the parameters file (JSON or TOML) associated with the collection HDF5 file."""
1452        # Check for TOML file first (preferred)
1453        toml_path = self.collection_hdf5_path.with_suffix(".toml")
1454        if toml_path.exists():
1455            return toml_path
1456
1457        # Check for JSON file
1458        json_path = self.collection_hdf5_path.with_suffix(".json")
1459        if json_path.exists():
1460            return json_path
1461
1462        # No parameters file found
1463        return None
1464
1465    def _load_collection_metadata(self):
1466        """Load metadata and manifest from the saved collection HDF5 file."""
1467        with h5py.File(self.collection_hdf5_path, "r") as f:
1468            self.folder_location = Path(f.attrs.get("lcms_objects_folder", ""))
1469            self.missing_mass_features_searched = f.attrs.get(
1470                "missing_mass_features_searched", False
1471            )
1472
1473            # Call the _load_manifest function to process the manifest
1474            self._manifest_dict = self._load_manifest(f)
1475
1476    def _load_manifest(self, hdf_handle):
1477        """Load and clean the manifest from the HDF5 file."""
1478        manifest_json = hdf_handle.attrs.get("manifest", "{}")
1479        if isinstance(manifest_json, bytes):
1480            manifest_json = manifest_json.decode("utf-8")
1481        loaded_manifest = json.loads(manifest_json)
1482
1483        # Convert integer values for 'use_rt_alignment' back to booleans
1484        def convert_back_to_bool(data):
1485            if isinstance(data, dict):
1486                # Process each key-value pair recursively
1487                return {
1488                    k: (
1489                        bool(v)
1490                        if k == "use_rt_alignment" and isinstance(v, int)
1491                        else convert_back_to_bool(v)
1492                    )
1493                    for k, v in data.items()
1494                }
1495            elif isinstance(data, list):
1496                # Recursively process lists
1497                return [convert_back_to_bool(item) for item in data]
1498            else:
1499                # Return non-dict/list types unchanged
1500                return data
1501
1502        # Clean the loaded manifest
1503        return convert_back_to_bool(loaded_manifest)
1504
1505    def _load_rt_alignments(self, lcms_collection):
1506        """Load retention time alignments from the saved collection HDF5 file."""
1507        # First Set the rt_aligned flag from the collection-level attribute saved directly
1508        with h5py.File(self.collection_hdf5_path, "r") as f:
1509            lcms_collection.rt_aligned = f.attrs.get("rt_aligned", False)
1510            lcms_collection.rt_alignment_attempted = f.attrs.get(
1511                "rt_alignment_attempted", False
1512            )
1513
1514        if lcms_collection.rt_aligned:
1515            with h5py.File(self.collection_hdf5_path, "r") as f:
1516                if "rt_alignments" in f:
1517                    # Iterate over the group `rt_alignments` containing datasets and add to the corresponding lcms object
1518                    rt_alignments_group = f["rt_alignments"]
1519                    for sample_idx, lcms_obj in zip(
1520                        rt_alignments_group.keys(), lcms_collection
1521                    ):
1522                        alignment_data = rt_alignments_group[sample_idx][:]
1523                        scan_df = lcms_obj.scan_df
1524                        scan_df["scan_time_aligned"] = alignment_data
1525                        lcms_obj.scan_df = scan_df
1526        elif lcms_collection.rt_alignment_attempted:
1527            # This means it was attempted and not used, so we populate the "scan_time_aligned"
1528            for lcms_obj in lcms_collection:
1529                scan_df = lcms_obj.scan_df
1530                scan_df["scan_time_aligned"] = scan_df["scan_time"]
1531                lcms_obj.scan_df = scan_df
1532
1533    def _load_cluster_assignments(self, lcms_collection):
1534        """Load cluster assignments from the saved collection HDF5 file."""
1535        with h5py.File(self.collection_hdf5_path, "r") as f:
1536            if "cluster_assignments" in f:
1537                # Access the group containing cluster assignments
1538                cluster_grp = f["cluster_assignments"]
1539
1540                # Reload index and cluster data
1541                index = cluster_grp["index"][:]  # Extract index
1542                index = [
1543                    idx.decode("utf-8") for idx in index
1544                ]  # Convert byte strings back to regular strings
1545                cluster_data = cluster_grp["cluster"][:]  # Extract cluster column
1546
1547                # Reassemble the DataFrame
1548                cluster_df = pd.DataFrame({"cluster": cluster_data}, index=index)
1549
1550                # Assign cluster data back to lcms_collection.mass_features_dataframe
1551                lcms_collection.mass_features_dataframe = (
1552                    lcms_collection.mass_features_dataframe.join(cluster_df, how="left")
1553                )
1554
1555                # Drop rows with NaN cluster values
1556                lcms_collection.mass_features_dataframe.dropna(
1557                    subset=["cluster"], inplace=True
1558                )
1559
1560    def _load_feature_group_labels(self, lcms_collection):
1561        """Restore ``feature_group_dataframe`` and merge onto feature tables."""
1562        from corems.mass_spectra.calc.feature_grouping import (
1563            GROUP_COLUMNS,
1564            empty_group_labels,
1565            merge_feature_group_labels_into_frames,
1566        )
1567
1568        with h5py.File(self.collection_hdf5_path, "r") as f:
1569            if "feature_group_labels" not in f:
1570                return
1571            grp = f["feature_group_labels"]
1572            if "cluster" not in grp:
1573                raise ValueError(
1574                    "feature_group_labels HDF5 group is missing 'cluster' dataset"
1575                )
1576            clusters = np.asarray(grp["cluster"][:], dtype=np.int64)
1577            labels = empty_group_labels(clusters)
1578            for col in GROUP_COLUMNS:
1579                if col not in grp:
1580                    continue
1581                raw = grp[col][:]
1582                if col in ("feature_group_id", "mono_cluster_id"):
1583                    vals = []
1584                    for v in raw:
1585                        if v is None or (isinstance(v, float) and np.isnan(v)):
1586                            vals.append(pd.NA)
1587                        else:
1588                            vals.append(int(v))
1589                    labels[col] = pd.Series(vals, index=labels.index, dtype="Int64")
1590                else:
1591                    decoded = []
1592                    for v in raw:
1593                        if isinstance(v, (bytes, bytearray, np.bytes_)):
1594                            s = v.decode("utf-8")
1595                        else:
1596                            s = "" if v is None else str(v)
1597                        decoded.append(None if s == "" else s)
1598                    labels[col] = decoded
1599
1600        lcms_collection.feature_group_dataframe = labels
1601        mf, induced = merge_feature_group_labels_into_frames(
1602            labels,
1603            lcms_collection.mass_features_dataframe,
1604            getattr(lcms_collection, "induced_mass_features_dataframe", None),
1605        )
1606        lcms_collection.mass_features_dataframe = mf
1607        if (
1608            getattr(lcms_collection, "induced_mass_features_dataframe", None)
1609            is not None
1610        ):
1611            lcms_collection.induced_mass_features_dataframe = induced
1612
1613    def get_lcms_collection(
1614        self,
1615        load_raw=False,
1616        load_light=False,
1617        load_representatives=False,
1618        load_eics=False,
1619        load_ms1=False,
1620        load_ms2=False,
1621    ):
1622        """Get the LCMS collection from the saved HDF5 file.
1623
1624        Parameters
1625        ----------
1626        load_raw : bool, optional
1627            If True, load raw data. Default is False.
1628        load_light : bool, optional
1629            If True, load light data (minimal). Default is False.
1630        load_representatives : bool, optional
1631            If True, load representative mass features from clusters. Default is False.
1632        load_eics : bool, optional
1633            If True, load EIC data for clustered mass features. Default is False.
1634        load_ms1 : bool, optional
1635            If True, load MS1 spectra for loaded mass features. Default is False.
1636        load_ms2 : bool, optional
1637            If True, load MS2 spectra for loaded mass features. Default is False.
1638
1639        Returns
1640        -------
1641        LCMSCollection
1642            The loaded LCMS collection object.
1643        """
1644        # First load the LCMSCollection object exactly as in the parent class
1645        lcms_collection = super().get_lcms_collection(
1646            load_raw=load_raw, load_light=load_light
1647        )
1648
1649        # Set the missing_mass_features_searched flag from saved metadata
1650        lcms_collection.missing_mass_features_searched = (
1651            self.missing_mass_features_searched
1652        )
1653
1654        # Load parameters if a parameters file exists
1655        if self.parameters_location:
1656            self._load_parameters(lcms_collection)
1657
1658        # Add retention time alignments if they exist
1659        self._load_rt_alignments(lcms_collection)
1660
1661        # Add cluster assignments if they exist
1662        self._load_cluster_assignments(lcms_collection)
1663
1664        # Load induced mass features if they exist
1665        self._load_induced_mass_features(lcms_collection)
1666
1667        # Load EICs for induced mass features from collection HDF5
1668        if lcms_collection.missing_mass_features_searched and load_eics:
1669            self._load_induced_eics_from_collection(lcms_collection)
1670
1671        # Combine induced mass features into the collection-level dataframe if any were loaded
1672        if lcms_collection.missing_mass_features_searched:
1673            lcms_collection._combine_mass_features(induced_features=True)
1674
1675        # Restore feature-group labels after clusters (+ induced) are in place
1676        self._load_feature_group_labels(lcms_collection)
1677
1678        # Load representative mass features if requested
1679        if load_representatives:
1680            self._load_representative_mass_features(lcms_collection)
1681
1682        # Load MS1 and/or MS2 spectra for loaded mass features if requested
1683        if load_ms1 or load_ms2:
1684            # Reuse the existing ReloadFeaturesOperation from the pipeline system
1685            from corems.mass_spectra.calc.lc_calc_operations import (
1686                ReloadFeaturesOperation,
1687            )
1688
1689            operations = [
1690                ReloadFeaturesOperation(
1691                    "reload_spectra", add_ms1=load_ms1, add_ms2=load_ms2
1692                )
1693            ]
1694            lcms_collection.process_samples_pipeline(
1695                operations, keep_raw_data=False, show_progress=False
1696            )
1697
1698        # Load EICs for clustered features if requested
1699        if load_eics:
1700            # Reuse the existing LoadEICsOperation from the pipeline system
1701            from corems.mass_spectra.calc.lc_calc_operations import LoadEICsOperation
1702
1703            operations = [LoadEICsOperation("load_eics")]
1704            lcms_collection.process_samples_pipeline(
1705                operations, keep_raw_data=False, show_progress=False
1706            )
1707
1708            # Associate EICs with mass features (same as in process_consensus_features)
1709            for sample_id in range(len(lcms_collection.samples)):
1710                sample = lcms_collection[sample_id]
1711                if sample.eics:  # Only if EICs were loaded
1712                    # Associate EICs with regular mass features
1713                    sample.associate_eics_with_mass_features(induced=False)
1714                    # Associate EICs with induced mass features
1715                    sample.associate_eics_with_mass_features(induced=True)
1716
1717        return lcms_collection
1718
1719    def _load_parameters(self, lcms_collection):
1720        """Load collection-level parameters from the saved parameters file."""
1721        from corems.encapsulation.input.parameter_from_json import (
1722            load_and_set_json_parameters_lcms_collection,
1723            load_and_set_toml_parameters_lcms_collection,
1724        )
1725
1726        if self.parameters_location.suffix == ".json":
1727            load_and_set_json_parameters_lcms_collection(
1728                lcms_collection, self.parameters_location
1729            )
1730        elif self.parameters_location.suffix == ".toml":
1731            load_and_set_toml_parameters_lcms_collection(
1732                lcms_collection, self.parameters_location
1733            )
1734        else:
1735            warnings.warn(
1736                f"Unknown parameter file format: {self.parameters_location.suffix}. Skipping parameter loading."
1737            )
1738
1739    def _load_induced_mass_features(self, lcms_collection):
1740        """Load induced mass features from the saved collection HDF5 file.
1741
1742        Induced mass features are gap-filled features that exist at the collection level.
1743        This method loads them from the collection HDF5 file with all their attributes
1744        and datasets, and distributes them to individual LCMS objects.
1745
1746        Parameters
1747        ----------
1748        lcms_collection : LCMSCollection
1749            The LCMS collection object to populate with induced mass features.
1750        """
1751        with h5py.File(self.collection_hdf5_path, "r") as f:
1752            if "induced_mass_features" not in f:
1753                return
1754
1755            # Access the top-level induced mass features group
1756            imf_group = f["induced_mass_features"]
1757
1758            # Iterate through each sample's induced mass features
1759            for sample_idx in imf_group.keys():
1760                lcms_obj = lcms_collection[int(sample_idx)]
1761                sample_group = imf_group[sample_idx]
1762
1763                # Load each mass feature for this sample
1764                for mf_id_str in sample_group.keys():
1765                    mf_group = sample_group[mf_id_str]
1766
1767                    # Note: Induced mass feature IDs are strings like 'c2923_0_i', not integers
1768                    # Keep them as strings since that's how they're stored
1769                    mf_id = mf_id_str
1770
1771                    # Instantiate the LCMSMassFeature object with required attributes
1772                    mass_feature = LCMSMassFeature(
1773                        lcms_obj,
1774                        mz=mf_group.attrs["_mz_exp"],
1775                        retention_time=mf_group.attrs["_retention_time"],
1776                        intensity=mf_group.attrs["_intensity"],
1777                        apex_scan=mf_group.attrs["_apex_scan"],
1778                        persistence=mf_group.attrs.get("_persistence", 0),
1779                        id=mf_id,
1780                    )
1781
1782                    # Populate additional attributes from HDF5 attributes
1783                    for key in mf_group.attrs.keys() - {
1784                        "_mz_exp",
1785                        "_mz_cal",
1786                        "_retention_time",
1787                        "_intensity",
1788                        "_apex_scan",
1789                        "_persistence",
1790                    }:
1791                        setattr(mass_feature, key, mf_group.attrs[key])
1792
1793                    # Populate attributes from HDF5 datasets (arrays)
1794                    for key in mf_group.keys():
1795                        setattr(mass_feature, key, mf_group[key][:])
1796                        # Convert _noise_score from array to tuple
1797                        if key == "_noise_score":
1798                            mass_feature._noise_score = tuple(mass_feature._noise_score)
1799
1800                    # Add to the LCMS object's induced_mass_features dictionary
1801                    lcms_obj.induced_mass_features[mf_id] = mass_feature
1802
1803    def _load_induced_eics_from_collection(self, lcms_collection):
1804        """Load EICs for induced mass features from the collection HDF5 file.
1805
1806        Induced mass features are gap-filled features. Their EICs are saved at the
1807        collection level and need to be loaded and associated with the induced mass features.
1808
1809        Parameters
1810        ----------
1811        lcms_collection : LCMSCollection
1812            The LCMS collection object with induced mass features to associate EICs with.
1813        """
1814        with h5py.File(self.collection_hdf5_path, "r") as f:
1815            if "induced_eics" not in f:
1816                return
1817
1818            # Access the top-level induced EICs group
1819            induced_eics_group = f["induced_eics"]
1820
1821            # Iterate through each sample's induced EICs
1822            for sample_idx in induced_eics_group.keys():
1823                lcms_obj = lcms_collection[int(sample_idx)]
1824                sample_group = induced_eics_group[sample_idx]
1825
1826                # Use the static helper to load EICs
1827                loaded_eics = ReadCoreMSHDFMassSpectra._load_eics_from_hdf5_group(
1828                    sample_group, lcms_obj
1829                )
1830
1831                # Ensure eics dictionary exists (should already be initialized in __init__)
1832                if not hasattr(lcms_obj, "eics") or lcms_obj.eics is None:
1833                    lcms_obj.eics = {}
1834
1835                # Add to lcms_obj.eics dictionary
1836                for eic_mz, eic in loaded_eics.items():
1837                    lcms_obj.eics[eic_mz] = eic
1838
1839            # Associate EICs with induced mass features after all samples processed
1840            # This is done outside the loop to handle all samples at once
1841            for lcms_obj in lcms_collection:
1842                if len(lcms_obj.induced_mass_features) > 0:
1843                    lcms_obj.associate_eics_with_mass_features(induced=True)
1844
1845    def _load_representative_mass_features(self, lcms_collection):
1846        """Load representative mass features for all clusters from HDF5 files.
1847
1848        This method uses the same logic as process_consensus_features() when loading
1849        representatives, calling get_sample_mf_map_for_representatives() (DRY helper)
1850        to determine which features to load.
1851
1852        Parameters
1853        ----------
1854        lcms_collection : LCMSCollection
1855            The LCMS collection object to populate with representative mass features.
1856        """
1857        # Get cluster assignments from the mass_features_dataframe
1858        if "cluster" not in lcms_collection.mass_features_dataframe.columns:
1859            return
1860
1861        # Use DRY helper method to build sample_mf_map with cluster IDs
1862        sample_mf_map = lcms_collection.get_sample_mf_map_for_representatives(
1863            include_cluster_id=True
1864        )
1865
1866        # Load mass features for each sample
1867        for sample_id, mf_list in sample_mf_map.items():
1868            lcms_obj = lcms_collection[sample_id]
1869
1870            # Load each mass feature
1871            for mf_id, cluster_id in mf_list:
1872                self._load_single_mass_feature(lcms_obj, mf_id, cluster_id)
1873
1874    def _load_single_mass_feature(self, lcms_obj, feature_id, cluster_index=None):
1875        """Load a single mass feature from an LCMS object's HDF5 file.
1876
1877        Parameters
1878        ----------
1879        lcms_obj : LCMSBase
1880            The LCMS object to add the mass feature to.
1881        feature_id : int
1882            The ID of the mass feature to load.
1883        cluster_index : int, optional
1884            The cluster index to assign to the loaded mass feature.
1885        """
1886        hdf5_path = lcms_obj.file_location.with_suffix(".hdf5")
1887
1888        if not hdf5_path.exists():
1889            return
1890
1891        with h5py.File(hdf5_path, "r") as f:
1892            if "mass_features" not in f:
1893                return
1894
1895            mf_group = f["mass_features"]
1896            feature_id_str = str(feature_id)
1897
1898            if feature_id_str not in mf_group:
1899                return
1900
1901            mf_data = mf_group[feature_id_str]
1902
1903            # Create LCMSMassFeature object
1904            mass_feature = LCMSMassFeature(
1905                lcms_obj,
1906                mz=mf_data.attrs["_mz_exp"],
1907                retention_time=mf_data.attrs["_retention_time"],
1908                intensity=mf_data.attrs["_intensity"],
1909                apex_scan=mf_data.attrs["_apex_scan"],
1910                persistence=mf_data.attrs.get("_persistence", 0),
1911                id=feature_id,
1912            )
1913
1914            # Set cluster_index if provided
1915            if cluster_index is not None:
1916                mass_feature.cluster_index = cluster_index
1917
1918            # Populate additional attributes from HDF5 attributes
1919            for key in mf_data.attrs.keys() - {
1920                "_mz_exp",
1921                "_mz_cal",
1922                "_retention_time",
1923                "_intensity",
1924                "_apex_scan",
1925                "_persistence",
1926            }:
1927                setattr(mass_feature, key, mf_data.attrs[key])
1928
1929            # Populate attributes from HDF5 datasets (arrays)
1930            for key in mf_data.keys():
1931                setattr(mass_feature, key, mf_data[key][:])
1932                # Convert _noise_score from array to tuple
1933                if key == "_noise_score":
1934                    mass_feature._noise_score = tuple(mass_feature._noise_score)
1935
1936            # Add to the LCMS object's mass_features dictionary
1937            lcms_obj.mass_features[feature_id] = mass_feature