corems.molecular_id.search.database_interfaces

   1import os
   2import re
   3import warnings
   4from abc import ABC
   5from io import StringIO
   6from pathlib import Path
   7import sqlite3
   8
   9import numpy as np
  10import requests
  11import pandas as pd
  12from ms_entropy import FlashEntropySearch
  13
  14from corems.molecular_id.factory.EI_SQL import (
  15    EI_LowRes_SQLite,
  16    Metadatar,
  17    MetaboliteMetadata,
  18)
  19from corems.molecular_id.factory.lipid_molecular_metadata import LipidMetadata
  20from corems.mass_spectra.calc.lc_calc import find_closest
  21
  22
  23class SpectralDatabaseInterface(ABC):
  24    """
  25    Base class that facilitates connection to spectral reference databases,
  26    such as EMSL's Metabolomics Reference Database (MetabRef).
  27
  28    """
  29
  30    def __init__(self, key=None):
  31        """
  32        Initialize instance.
  33
  34        Parameters
  35        ----------
  36        key : str
  37            Token key.
  38
  39        """
  40
  41        self.key = key
  42
  43    def set_token(self, path):
  44        """
  45        Set environment variable for MetabRef database token.
  46
  47        Parameters
  48        ----------
  49        path : str
  50            Path to token.
  51
  52        """
  53
  54        # Read token from file
  55        with open(path, "r", encoding="utf-8") as f:
  56            token = f.readline().strip()
  57
  58        # Set environment variable
  59        os.environ[self.key] = token
  60
  61    def get_token(self):
  62        """
  63        Get environment variable for database token.
  64
  65        Returns
  66        -------
  67        str
  68            Token string.
  69
  70        """
  71
  72        # Check for token
  73        if self.key not in os.environ:
  74            raise ValueError("Must set {} environment variable.".format(self.key))
  75
  76        # Get token from environment variables
  77        return os.environ.get(self.key)
  78
  79    def get_header(self):
  80        """
  81        Access stored database token and prepare as header.
  82
  83        Returns
  84        -------
  85        str
  86            Header string.
  87
  88        """
  89
  90        # Get token
  91        token = self.get_token()
  92
  93        # Pad header information
  94        header = {"Authorization": f"Bearer {token}", "Content-Type": "text/plain"}
  95
  96        return header
  97
  98    def get_query(self, url, use_header=True):
  99        """
 100        Request payload from URL according to `get` protocol.
 101
 102        Parameters
 103        ----------
 104        url : str
 105            URL for request.
 106        use_header: bool
 107            Whether or not the query should include the header
 108
 109        Returns
 110        -------
 111        dict
 112            Response as JSON.
 113
 114        """
 115
 116        # Query URL via `get`
 117        if use_header:
 118            response = requests.get(url, headers=self.get_header())
 119        else:
 120            response = requests.get(url)
 121
 122        # Check response
 123        response.raise_for_status()
 124
 125        # Return as JSON
 126        return response.json()
 127
 128    def post_query(self, url, variable, values, tolerance):
 129        """
 130        Request payload from URL according to `post` protocol.
 131
 132        Parameters
 133        ----------
 134        url : str
 135            URL for request.
 136        variable : str
 137            Variable to query.
 138        values : str
 139            Specific values of `variable` to query.
 140        tolerance : str
 141            Query tolerance relative to `values`.
 142
 143        Returns
 144        -------
 145        dict
 146            Response as JSON.
 147
 148        """
 149
 150        # Coerce to string
 151        if not isinstance(variable, str):
 152            variable = str(variable).replace(" ", "")
 153
 154        if not isinstance(values, str):
 155            values = str(values).replace(" ", "")
 156
 157        if not isinstance(tolerance, str):
 158            tolerance = str(tolerance).replace(" ", "")
 159
 160        # Query URL via `post`
 161        response = requests.post(
 162            os.path.join(url, variable, tolerance),
 163            data=values,
 164            headers=self.get_header(),
 165        )
 166
 167        # Check response
 168        response.raise_for_status()
 169
 170        # Return as JSON
 171        return response.json()
 172
 173    def _check_flash_entropy_kwargs(self, fe_kwargs):
 174        """
 175        Check FlashEntropy keyword arguments.
 176
 177        Parameters
 178        ----------
 179        fe_kwargs : dict
 180            Keyword arguments for FlashEntropy search.
 181
 182
 183        Raises
 184        ------
 185        ValueError
 186            If "min_ms2_difference_in_da" or "max_ms2_tolerance_in_da" are present in `fe_kwargs` and they
 187            are not equal.
 188
 189        """
 190        # If "min_ms2_difference_in_da" in fe_kwargs, check that "max_ms2_tolerance_in_da" is also present and that min_ms2_difference_in_da = 2xmax_ms2_tolerance_in_da
 191        if (
 192            "min_ms2_difference_in_da" in fe_kwargs
 193            or "max_ms2_tolerance_in_da" in fe_kwargs
 194        ):
 195            if (
 196                "min_ms2_difference_in_da" not in fe_kwargs
 197                or "max_ms2_tolerance_in_da" not in fe_kwargs
 198            ):
 199                raise ValueError(
 200                    "Both 'min_ms2_difference_in_da' and 'max_ms2_tolerance_in_da' must be specified."
 201                )
 202            if (
 203                fe_kwargs["min_ms2_difference_in_da"]
 204                != 2 * fe_kwargs["max_ms2_tolerance_in_da"]
 205            ):
 206                raise ValueError(
 207                    "The values of 'min_ms2_difference_in_da' must be exactly 2x 'max_ms2_tolerance_in_da'."
 208                )
 209
 210    def _get_format_func(self, format):
 211        """
 212        Obtain format function by key.
 213
 214        Returns
 215        -------
 216        func
 217            Formatting function.
 218        """
 219
 220        if format.lower() in self.format_map.keys():
 221            return self.format_map[format.lower()]
 222
 223        raise ValueError(("{} not a supported format.").format(format))
 224
 225    def _dict_to_dataclass(self, source_dict, data_class):
 226        """
 227        Convert dictionary to dataclass.
 228
 229        Notes
 230        -----
 231        This function will pull the attributes a dataclass and its parent class
 232        and convert the dictionary to a dataclass instance with the appropriate
 233        attributes.
 234
 235        Parameters
 236        ----------
 237        data_class : :obj:`~dataclasses.dataclass`
 238            Dataclass to convert to.
 239        source_dict : dict
 240            Dictionary object to convert to dataclass.
 241
 242        Returns
 243        -------
 244        :obj:`~dataclasses.dataclass`
 245            Dataclass instance.
 246
 247        """
 248
 249        # Get list of expected attributes of data_class
 250        data_class_keys = list(data_class.__annotations__.keys())
 251
 252        # Does the data_class inherit from another class, if so, get the attributes of the parent class as well
 253        if len(data_class.__mro__) > 2:
 254            parent_class_keys = list(data_class.__bases__[0].__annotations__.keys())
 255            data_class_keys = list(set(data_class_keys + parent_class_keys))
 256
 257        # Remove keys that are not in the data_class from the input dictionary
 258        input_dict = {k: v for k, v in source_dict.items() if k in data_class_keys}
 259
 260        # Add keys that are in the data class but not in the input dictionary as None
 261        for key in data_class_keys:
 262            if key not in input_dict.keys():
 263                input_dict[key] = None
 264        return data_class(**input_dict)
 265
 266    def _spectrum_to_array(self, spectrum, normalize=True):
 267        """
 268        Convert a parenthesis-delimited spectrum string to array.
 269
 270        Parameters
 271        ----------
 272        spectrum : str
 273            Spectrum string, i.e. list of (m/z,abundance) pairs.
 274        normalize : bool
 275            Normalize the spectrum by its magnitude.
 276
 277        Returns
 278        -------
 279        :obj:`~numpy.array`
 280            Array of shape (N, 2), with m/z in the first column and abundance in
 281            the second.
 282        """
 283
 284        arr = np.array(
 285            re.findall(r"\(([^,]+),([^)]+)\)", spectrum), dtype=float
 286        ).reshape(-1, 2)
 287
 288        if normalize:
 289            arr = self.normalize_peaks(arr)
 290
 291        return arr
 292
 293    @staticmethod
 294    def normalize_peaks(arr):
 295        """
 296        Normalize peaks in an array.
 297
 298        Parameters
 299        ----------
 300        arr : :obj:`~numpy.array`
 301            Array of shape (N, 2), with m/z in the first column and abundance in
 302            the second.
 303
 304        Returns
 305        -------
 306        :obj:`~numpy.array`
 307            Normalized array of shape (N, 2), with m/z in the first column and
 308            normalized abundance in the second.
 309        """
 310        # Normalize the array
 311        arr[:, -1] = arr[:, -1] / arr[:, -1].sum()
 312
 313        return arr
 314
 315    @staticmethod
 316    def _build_flash_entropy_index(fe_lib, fe_kwargs={}, clean_spectra=True):
 317        """
 318        Build FlashEntropy index.
 319
 320        Parameters
 321        ----------
 322        fe_lib : list
 323            List of spectra to build index from. Can be a list of dictionaries or
 324            a FlashEntropy search instance.
 325        fe_kwargs : dict, optional
 326            Keyword arguments for FlashEntropy search.
 327        clean_spectra : bool, optional
 328            Clean spectra before building index. Default is True.
 329
 330        Returns
 331        -------
 332        :obj:`~ms_entropy.FlashEntropySearch`
 333            FlashEntropy search instance.
 334
 335        """
 336        # Initialize FlashEntropy
 337        fe_init_kws = [
 338            "max_ms2_tolerance_in_da",
 339            "mz_index_step",
 340            "low_memory",
 341            "path_data",
 342        ]
 343        fe_init_kws = {k: v for k, v in fe_kwargs.items() if k in fe_init_kws}
 344        fes = FlashEntropySearch(**fe_init_kws)
 345
 346        # Build FlashEntropy index
 347        fe_index_kws = [
 348            "max_indexed_mz",
 349            "precursor_ions_removal_da",
 350            "noise_threshold",
 351            "min_ms2_difference_in_da",
 352            "max_peak_num",
 353        ]
 354        fe_index_kws = {k: v for k, v in fe_kwargs.items() if k in fe_index_kws}
 355        fes.build_index(fe_lib, **fe_index_kws, clean_spectra=clean_spectra)
 356
 357        return fes
 358
 359
 360class MetabRefInterface(SpectralDatabaseInterface):
 361    """
 362    DEPRECATED interface retained for backward compatibility only.
 363    """
 364
 365    def __init__(self):
 366        """
 367        Initialize instance with deprecation warning.
 368
 369        """
 370
 371        super().__init__(key=None)
 372
 373        if self.__class__ is MetabRefInterface:
 374            warnings.warn(
 375                "MetabRefInterface is deprecated. Instantiate a concrete interface "
 376                "such as GCMSLibraryInterface or LCLipidLibraryInterface instead.",
 377                DeprecationWarning,
 378                stacklevel=2,
 379            )
 380
 381
 382class GCMSLibraryInterface(SpectralDatabaseInterface):
 383    """
 384    Interface to bundled GCMS spectral libraries in MSP format.
 385    
 386    Loads GCMS compound library and FAMES calibration library from local MSP files.
 387    Default files are bundled with CoreMS, but can be overridden via environment variables.
 388    """
 389
 390    def __init__(self):
 391        """
 392        Initialize instance.
 393        """
 394        super().__init__(key=None)
 395        
 396        # Local data file paths
 397        from pathlib import Path
 398        
 399        # Default to bundled data files
 400        data_dir = Path(__file__).parent.parent / "data"
 401        self.gcms_library_file = os.getenv(
 402            "GCMS_LIBRARY_PATH", 
 403            str(data_dir / "PNNLMetV20191015.msp")
 404        )
 405        self.fames_library_file = os.getenv(
 406            "FAMES_LIBRARY_PATH",
 407            str(data_dir / "FAMES_REF.msp")
 408        )
 409
 410        self.__init_format_map__()
 411
 412    def __init_format_map__(self):
 413        """
 414        Initialize database format mapper, enabling multiple format requests.
 415
 416        """
 417
 418        # Define format workflows
 419        self.format_map = {
 420            "json": lambda x, normalize, fe_kwargs: x,
 421            "dict": lambda x,
 422            normalize,
 423            fe_kwargs: self._to_LowResolutionEICompound_dict(x, normalize),
 424            "sql": lambda x,
 425            normalize,
 426            fe_kwargs: self._LowResolutionEICompound_dict_to_sqlite(
 427                self._to_LowResolutionEICompound_dict(x, normalize)
 428            ),
 429        }
 430
 431        # Add aliases
 432        self.format_map["metabref"] = self.format_map["json"]
 433        self.format_map["datadict"] = self.format_map["dict"]
 434        self.format_map["data-dict"] = self.format_map["dict"]
 435        self.format_map["lowreseicompound"] = self.format_map["dict"]
 436        self.format_map["lowres"] = self.format_map["dict"]
 437        self.format_map["lowresgc"] = self.format_map["dict"]
 438        self.format_map["sqlite"] = self.format_map["sql"]
 439
 440    def available_formats(self):
 441        """
 442        View list of available formats.
 443
 444        Returns
 445        -------
 446        list
 447            Format map keys.
 448        """
 449
 450        return list(self.format_map.keys())
 451
 452    def get_library(self, format="json", normalize=False):
 453        """
 454        Load GC/MS library from local MSP file.
 455
 456        Parameters
 457        ----------
 458        format : str
 459            Format of requested library, i.e. "json", "sql", "dict".
 460            See `available_formats` method for aliases.
 461        normalize : bool
 462            Normalize the spectrum by its magnitude.
 463
 464        Returns
 465        -------
 466        Library in requested format.
 467
 468        """
 469        # Load from local MSP file
 470        library_data = self._load_msp_file(self.gcms_library_file, normalize)
 471        
 472        # Init format function
 473        format_func = self._get_format_func(format)
 474        
 475        # Apply format conversion
 476        return format_func(library_data, normalize, {})
 477
 478    def get_fames(self, format="json", normalize=False):
 479        """
 480        Load GC/MS FAMEs library from local MSP file.
 481
 482        Parameters
 483        ----------
 484        format : str
 485            Format of requested library, i.e. "json", "sql", "dict".
 486            See `available_formats` method for aliases.
 487        normalize : bool
 488            Normalize the spectrum by its magnitude.
 489
 490        Returns
 491        -------
 492        Library in requested format.
 493
 494        """
 495        # Load from local MSP file
 496        library_data = self._load_msp_file(self.fames_library_file, normalize)
 497        
 498        # Init format function
 499        format_func = self._get_format_func(format)
 500        
 501        # Apply format conversion
 502        return format_func(library_data, normalize, {})
 503    
 504    def _load_msp_file(self, file_path, normalize=False):
 505        """
 506        Load and parse MSP file into format compatible with existing pipeline.
 507        
 508        Parameters
 509        ----------
 510        file_path : str
 511            Path to MSP file
 512        normalize : bool
 513            Normalize spectra
 514            
 515        Returns
 516        -------
 517        list of dict
 518            Library data in format compatible with _to_LowResolutionEICompound_dict
 519        """
 520        from pathlib import Path
 521        
 522        file_path = Path(file_path)
 523        if not file_path.exists():
 524            raise FileNotFoundError(
 525                f"Library file not found: {file_path}. "
 526                f"Set GCMS_LIBRARY_PATH or FAMES_LIBRARY_PATH environment variable to specify location."
 527            )
 528        
 529        # Parse MSP file
 530        spectra = []
 531        spectrum = {}
 532        peaks = []
 533        
 534        with open(file_path, 'r') as f:
 535            for line in f:
 536                line = line.strip()
 537                
 538                # Empty line marks end of spectrum
 539                if not line:
 540                    if spectrum and peaks:
 541                        # Convert peaks to the format expected by downstream code
 542                        # Format: "(mz,abundance)(mz,abundance)..."
 543                        peak_str = "".join([f"({int(mz)},{int(abun)})" for mz, abun in peaks])
 544                        spectrum['mz'] = peak_str
 545                        spectra.append(spectrum)
 546                    spectrum = {}
 547                    peaks = []
 548                    continue
 549                
 550                # Check if line contains peak data (starts with digit)
 551                if line and line[0].isdigit():
 552                    parts = line.split()
 553                    if len(parts) >= 2:
 554                        peaks.append((float(parts[0]), float(parts[1])))
 555                    continue
 556                
 557                # Handle metadata fields
 558                if ":" in line:
 559                    key, value = line.split(":", 1)
 560                    key = key.strip().lower()
 561                    value = value.strip()
 562                    
 563                    # Map MSP fields to expected format
 564                    field_mapping = {
 565                        "name": "molecule_name",
 566                        "formula": "formula",
 567                        "cas": "casno",
 568                        "retentiontime": "retention_time",
 569                        "ri": "ri",
 570                        "comment": "comments",
 571                        "num peaks": "peak_count",
 572                        "derivative": "derivative"
 573                    }
 574                    
 575                    # Metadata fields that go into the metadata dict
 576                    metadata_fields = {
 577                        "inchikey": "inchikey",
 578                        "inchi": "inchi",
 579                        "smiles": "smiles",
 580                        "pubchem": "pubchem",
 581                        "chebi": "chebi",
 582                        "kegg": "kegg",
 583                        "refmet": "refmet",
 584                        "iupac_name": "iupac_name"
 585                    }
 586                    
 587                    if key in field_mapping:
 588                        mapped_key = field_mapping[key]
 589                        # Convert numeric fields
 590                        if key in ["retentiontime", "ri"]:
 591                            try:
 592                                value = float(value)
 593                            except:
 594                                pass
 595                        elif key == "num peaks":
 596                            try:
 597                                value = int(value)
 598                            except:
 599                                pass
 600                        spectrum[mapped_key] = value
 601                    elif key in metadata_fields:
 602                        # Store in nested metadata dict
 603                        if "metadata" not in spectrum:
 604                            spectrum["metadata"] = {}
 605                        spectrum["metadata"][metadata_fields[key]] = value
 606                    else:
 607                        # Keep unmapped fields
 608                        spectrum[key] = value
 609        
 610        # Add last spectrum if file doesn't end with blank line
 611        if spectrum and peaks:
 612            peak_str = "".join([f"({int(mz)},{int(abun)})" for mz, abun in peaks])
 613            spectrum['mz'] = peak_str
 614            spectra.append(spectrum)
 615        
 616        return spectra
 617
 618    def _to_LowResolutionEICompound_dict(self, metabref_lib, normalize=False):
 619        """
 620        Convert MetabRef-formatted library to CoreMS LowResolutionEICompound-formatted
 621        dictionary for local ingestion.
 622
 623        Parameters
 624        ----------
 625        metabref_lib : dict
 626            MetabRef GC-MS library in JSON format.
 627        normalize : bool
 628            Normalize each spectrum by its magnitude.
 629
 630        Returns
 631        -------
 632        list of dict
 633            List of each spectrum contained in dictionary.
 634
 635        """
 636
 637        # All below key:value lookups are based on CoreMS class definitions
 638        # NOT MetabRef content. For example, MetabRef has keys for PubChem,
 639        # USI, etc. that are not considered below.
 640
 641        # Dictionary to map metabref keys to corems keys
 642        metadatar_cols = {
 643            "casno": "cas",
 644            "inchikey": "inchikey",
 645            "inchi": "inchi",
 646            "chebi": "chebi",
 647            "smiles": "smiles",
 648            "kegg": "kegg",
 649            "iupac_name": "iupac_name",
 650            "traditional_name": "traditional_name",  # Not present in metabref
 651            "common_name": "common_name",  # Not present in metabref
 652        }
 653
 654        # Dictionary to map metabref keys to corems keys
 655        lowres_ei_compound_cols = {
 656            "id": "metabref_id",
 657            "molecule_name": "name",  # Is this correct?
 658            "classify": "classify",  # Not present in metabref
 659            "formula": "formula",
 660            "ri": "ri",
 661            "rt": "retention_time",
 662            "source": "source",  # Not present in metabref
 663            "casno": "casno",
 664            "comments": "comment",
 665            "source_temp_c": "source_temp_c",  # Not present in metabref
 666            "ev": "ev",  # Not present in metabref
 667            "peak_count": "peaks_count",
 668            "mz": "mz",
 669            "abundance": "abundance",
 670        }
 671
 672        # Local result container
 673        corems_lib = []
 674
 675        # Enumerate spectra
 676        for i, source_ in enumerate(metabref_lib):
 677            # Copy source to prevent modification
 678            source = source_.copy()
 679
 680            # Parse target data
 681            target = {
 682                lowres_ei_compound_cols[k]: v
 683                for k, v in source.items()
 684                if k in lowres_ei_compound_cols
 685            }
 686
 687            # Explicitly add this to connect with LowResCompoundRef later
 688            if "retention_time" in source:
 689                target["rt"] = source["retention_time"]
 690            elif "rt" in source:
 691                target["rt"] = source["rt"]
 692
 693            # Parse (mz, abundance)
 694            arr = self._spectrum_to_array(target["mz"], normalize=normalize)
 695            target["mz"] = arr[:, 0]
 696            target["abundance"] = arr[:, 1]
 697
 698            # Parse meta data
 699            target["metadata"] = {
 700                metadatar_cols[k]: v for k, v in source.items() if k in metadatar_cols
 701            }
 702
 703            # Add anything else
 704            for k in source:
 705                if k not in lowres_ei_compound_cols:
 706                    target[k] = source[k]
 707
 708            # Add to CoreMS list
 709            corems_lib.append(target)
 710
 711        return corems_lib
 712
 713    def _LowResolutionEICompound_dict_to_sqlite(
 714        self, lowres_ei_compound_dict, url="sqlite://"
 715    ):
 716        """
 717        Convert CoreMS LowResolutionEICompound-formatted dictionary to SQLite
 718        database for local ingestion.
 719
 720        Parameters
 721        ----------
 722        lowres_ei_compound_dict : dict
 723            CoreMS GC-MS library formatted for LowResolutionEICompound.
 724        url : str
 725            URL to SQLite prefix.
 726
 727        Returns
 728        -------
 729        sqlite database
 730            Spectra contained in SQLite database.
 731
 732        """
 733
 734        # Dictionary to map corems keys to all-caps keys
 735        capped_cols = {
 736            "name": "NAME",
 737            "formula": "FORM",
 738            "ri": "RI",
 739            "retention_time": "RT",
 740            "source": "SOURCE",
 741            "casno": "CASNO",
 742            "comment": "COMMENT",
 743            "peaks_count": "NUM PEAKS",
 744        }
 745
 746        # Initialize SQLite object
 747        sqlite_obj = EI_LowRes_SQLite(url=url)
 748
 749        # Iterate spectra
 750        for _data_dict in lowres_ei_compound_dict:
 751            # Copy source to prevent modification
 752            data_dict = _data_dict.copy()
 753
 754            # Add missing capped values
 755            for k, v in capped_cols.items():
 756                # Key exists
 757                if k in data_dict:
 758                    # # This will replace the key
 759                    # data_dict[v] = data_dict.pop(k)
 760
 761                    # This will keep both keys
 762                    data_dict[v] = data_dict[k]
 763
 764            # Parse number of peaks
 765            if not data_dict.get("NUM PEAKS"):
 766                data_dict["NUM PEAKS"] = len(data_dict.get("mz"))
 767
 768            # Parse CAS number
 769            if not data_dict.get("CASNO"):
 770                data_dict["CASNO"] = data_dict.get("CAS")
 771
 772            if not data_dict["CASNO"]:
 773                data_dict["CASNO"] = 0
 774
 775            # Build linked metadata table
 776            if "metadata" in data_dict:
 777                metadata = data_dict.pop("metadata")
 778                # Only create metadata entry if we have required fields and valid data
 779                # Filter to only include fields that Metadatar model supports
 780                supported_metadata_fields = [
 781                    'cas', 'inchikey', 'inchi', 'chebi', 'smiles', 
 782                    'kegg', 'iupac_name', 'traditional_name', 'common_name'
 783                ]
 784                filtered_metadata = {
 785                    k: v for k, v in metadata.items() 
 786                    if k in supported_metadata_fields and v
 787                }
 788                # Inchikey is required by the database model
 789                if filtered_metadata and filtered_metadata.get("inchikey"):
 790                    data_dict["metadatar"] = Metadatar(**filtered_metadata)
 791
 792            # Attempt addition to sqlite
 793            try:
 794                sqlite_obj.add_compound(data_dict)
 795            except:
 796                print(data_dict["NAME"])
 797
 798        return sqlite_obj
 799
 800
 801class MetabRefGCInterface(GCMSLibraryInterface):
 802    """
 803    DEPRECATED: Use GCMSLibraryInterface instead.
 804    
 805    This interface is maintained for backward compatibility only.
 806    MetabRef API has been discontinued as of 2026.
 807    """
 808
 809    def __init__(self):
 810        """
 811        Initialize instance with deprecation warning.
 812        """
 813        warnings.warn(
 814            "MetabRefGCInterface is deprecated. Use GCMSLibraryInterface instead. "
 815            "MetabRef API has been discontinued; all data now loads from bundled local MSP files.",
 816            DeprecationWarning,
 817            stacklevel=2
 818        )
 819        super().__init__()
 820
 821
 822class LCLipidLibraryInterface(SpectralDatabaseInterface):
 823    """
 824    Interface to a local sqlite lipid library for LC-MS spectral searches.
 825    """
 826
 827    DEFAULT_DOWNLOAD_URL = (
 828        "https://nmdcdemo.emsl.pnnl.gov/minio/lipidomics/parameter_files/"
 829        "202412_lipid_ref.sqlite"
 830    )
 831
 832    def __init__(self, db_location=None):
 833        """
 834        Initialize instance.
 835
 836        Parameters
 837        ----------
 838        db_location : str | Path, optional
 839            Local path to the sqlite lipid library. If omitted, the
 840            COREMS_LIPIDOMICS_SQLITE_PATH environment variable is used.
 841        """
 842
 843        super().__init__(key=None)
 844        self.db_location = db_location
 845        self.__init_format_map__()
 846
 847    def _to_flashentropy(self, spectral_library, normalize=True, fe_kwargs={}):
 848        """
 849        Convert a spectral library to FlashEntropy format.
 850
 851        Parameters
 852        ----------
 853        spectral_library : dict
 854            MS2 library in JSON format or FlashEntropy search instance
 855            (for reformatting at different MS2 separation).
 856        normalize : bool
 857            Normalize each spectrum by its magnitude.
 858        fe_kwargs : dict, optional
 859            Keyword arguments for instantiation of FlashEntropy search and building index for FlashEntropy search;
 860            any keys not recognized will be ignored. By default, all parameters set to defaults.
 861
 862        Returns
 863        -------
 864        :obj:`~ms_entropy.FlashEntropySearch`
 865            MS2 library as FlashEntropy search instance.
 866
 867        Raises
 868        ------
 869        ValueError
 870            If "min_ms2_difference_in_da" or "max_ms2_tolerance_in_da" are present in `fe_kwargs` and they are not equal.
 871
 872        """
 873        self._check_flash_entropy_kwargs(fe_kwargs)
 874
 875        # Initialize empty library
 876        fe_lib = []
 877
 878        # Enumerate spectra
 879        for i, source in enumerate(spectral_library):
 880            if "spectrum_data" in source.keys():
 881                spectrum = source["spectrum_data"]
 882            else:
 883                spectrum = source
 884
 885            if "precursor_mz" not in spectrum.keys():
 886                spectrum["precursor_mz"] = spectrum.pop("precursor_ion")
 887
 888            spectrum["peaks"] = self._spectrum_to_array(
 889                spectrum["mz"], normalize=normalize
 890            )
 891            fe_lib.append(spectrum)
 892
 893        fe_search = self._build_flash_entropy_index(fe_lib, fe_kwargs=fe_kwargs)
 894
 895        return fe_search
 896
 897    def __init_format_map__(self):
 898        """
 899        Initialize database format mapper, enabling multiple format requests.
 900        """
 901
 902        self.format_map = {
 903            "json": lambda x, normalize, fe_kwargs: x,
 904            "flashentropy": lambda x, normalize, fe_kwargs: self._to_flashentropy(
 905                x, normalize, fe_kwargs
 906            ),
 907            "dataframe": lambda x, normalize, fe_kwargs: pd.DataFrame(x),
 908        }
 909
 910        self.format_map["fe"] = self.format_map["flashentropy"]
 911        self.format_map["flash-entropy"] = self.format_map["flashentropy"]
 912        self.format_map["df"] = self.format_map["dataframe"]
 913
 914    def available_formats(self):
 915        """
 916        View list of available formats.
 917
 918        Returns
 919        -------
 920        list
 921            Format map keys.
 922        """
 923
 924        return list(self.format_map.keys())
 925
 926    def _resolve_db_location(self):
 927        """
 928        Resolve and validate sqlite database location.
 929
 930        Returns
 931        -------
 932        Path
 933            Existing sqlite database file path.
 934        """
 935
 936        db_location = self.db_location or os.getenv("COREMS_LIPIDOMICS_SQLITE_PATH")
 937        if not db_location:
 938            raise ValueError(
 939                "A local lipid sqlite library path is required. "
 940                "Set COREMS_LIPIDOMICS_SQLITE_PATH or pass db_location."
 941            )
 942
 943        db_path = Path(db_location).expanduser()
 944        if not db_path.exists():
 945            raise FileNotFoundError(
 946                f"Lipid sqlite library not found at {db_path}. "
 947                f"Download it from {self.DEFAULT_DOWNLOAD_URL} "
 948                "and set COREMS_LIPIDOMICS_SQLITE_PATH."
 949            )
 950
 951        return db_path
 952
 953    def _get_candidate_spectra(self, connection, mz_list, polarity, mz_tol_ppm):
 954        """
 955        Fetch candidate spectra rows by precursor m/z and polarity.
 956
 957        Returns
 958        -------
 959        pandas.DataFrame
 960            Filtered rows from lipidMassSpectrumObject.
 961        """
 962
 963        mz_observed = np.sort(np.asarray(mz_list, dtype=float))
 964        if mz_observed.size == 0:
 965            return pd.DataFrame()
 966
 967        mz_all = pd.read_sql_query(
 968            "SELECT id, polarity, precursor_mz FROM lipidMassSpectrumObject", connection
 969        )
 970        mz_all = mz_all[mz_all["polarity"] == polarity].copy()
 971        if mz_all.empty:
 972            return pd.DataFrame()
 973
 974        mz_all = mz_all.sort_values(by="precursor_mz").reset_index(drop=True)
 975
 976        if mz_observed.size == 1:
 977            mz_all["closest_mz_obs"] = mz_observed[0]
 978        else:
 979            mz_all["closest_mz_obs"] = mz_observed[
 980                find_closest(mz_observed, mz_all.precursor_mz.values)
 981            ]
 982
 983        mz_all["ppm_error"] = (
 984            (mz_all["precursor_mz"] - mz_all["closest_mz_obs"])
 985            / mz_all["precursor_mz"]
 986            * 1e6
 987        )
 988
 989        mz_all = mz_all[np.abs(mz_all["ppm_error"]) <= mz_tol_ppm]
 990        if mz_all.empty:
 991            return pd.DataFrame()
 992
 993        mz_ids = tuple(mz_all["id"].tolist())
 994        return pd.read_sql_query(
 995            f"SELECT * FROM lipidMassSpectrumObject WHERE id IN {mz_ids}",
 996            connection,
 997        )
 998
 999    def get_lipid_library(
1000        self,
1001        mz_list,
1002        polarity,
1003        mz_tol_ppm,
1004        mz_tol_da_api=None,
1005        format="json",
1006        normalize=True,
1007        fe_kwargs={},
1008        api_delay=5,
1009        api_attempts=10,
1010    ):
1011        """
1012        Retrieve lipid spectra and metadata from a local sqlite library.
1013
1014        Parameters
1015        ----------
1016        mz_list : list
1017            List of precursor m/z values.
1018        polarity : str
1019            Ionization polarity, either "positive" or "negative".
1020        mz_tol_ppm : float
1021            Tolerance in ppm for precursor matching.
1022        mz_tol_da_api : float, optional
1023            Unused, kept for backward compatibility.
1024        format : str, optional
1025            Format of requested library, e.g. "json" or "flashentropy".
1026        normalize : bool, optional
1027            Normalize spectrum intensities.
1028        fe_kwargs : dict, optional
1029            Keyword arguments for FlashEntropy search.
1030        api_delay : int, optional
1031            Unused, kept for backward compatibility.
1032        api_attempts : int, optional
1033            Unused, kept for backward compatibility.
1034
1035        Returns
1036        -------
1037        tuple
1038            Library in requested format and lipid metadata dictionary.
1039        """
1040
1041        if not isinstance(mz_list, (list, np.ndarray)):
1042            raise ValueError("mz_list must be a list or numpy array")
1043        if not all(isinstance(mz, (float, int)) for mz in mz_list):
1044            raise ValueError("All elements in mz_list must be float or int")
1045        if polarity not in {"positive", "negative"}:
1046            raise ValueError("polarity must be either 'positive' or 'negative'")
1047        if not isinstance(mz_tol_ppm, (float, int)):
1048            raise ValueError("mz_tol_ppm must be a float or int")
1049
1050        db_path = self._resolve_db_location()
1051        connection = sqlite3.connect(str(db_path))
1052        try:
1053            # Step 1: Get candidate spectra records based on m/z and polarity
1054            spectra_df = self._get_candidate_spectra(
1055                connection=connection,
1056                mz_list=mz_list,
1057                polarity=polarity,
1058                mz_tol_ppm=float(mz_tol_ppm),
1059            )
1060
1061            if spectra_df.empty:
1062                format_func = self._get_format_func(format)
1063                return format_func([], normalize=normalize, fe_kwargs=fe_kwargs), {}
1064
1065            # Step 2: Get corresponding lipid metadata for candidate spectra from lipidTree view
1066            mol_ids = tuple(spectra_df["molecular_data_id"].tolist())
1067            mol_df = pd.read_sql_query(
1068                f"SELECT * FROM lipidTree WHERE id IN {mol_ids}",
1069                connection,
1070            )
1071        finally:
1072            connection.close()
1073
1074        mol_df["id_index"] = mol_df["id"]
1075        mol_df = mol_df.set_index("id_index")
1076        mol_records = mol_df.to_dict(orient="index")
1077        lipid_metadata = {
1078            int(k): self._dict_to_dataclass(v, LipidMetadata)
1079            for k, v in mol_records.items()
1080        }
1081
1082        spectra_records = spectra_df.to_dict(orient="records")
1083        format_func = self._get_format_func(format)
1084        library = format_func(spectra_records, normalize=normalize, fe_kwargs=fe_kwargs)
1085        return library, lipid_metadata
1086
1087
1088class MSPInterface(SpectralDatabaseInterface):
1089    """
1090    Interface to parse NIST MSP files
1091    """
1092
1093    def __init__(self, file_path):
1094        """
1095        Initialize instance.
1096
1097        Parameters
1098        ----------
1099        file_path : str
1100            Path to a local MSP file.
1101
1102        Attributes
1103        ----------
1104        file_path : str
1105            Path to the MSP file.
1106        _file_content : str
1107            Content of the MSP file.
1108        _data_frame : :obj:`~pandas.DataFrame`
1109            DataFrame of spectra from the MSP file with unaltered content.
1110        """
1111        super().__init__(key=None)
1112
1113        self.file_path = file_path
1114        if not os.path.exists(self.file_path):
1115            raise FileNotFoundError(
1116                f"File {self.file_path} does not exist. Please check the file path."
1117            )
1118        with open(self.file_path, "r") as f:
1119            self._file_content = f.read()
1120
1121        self._data_frame = self._read_msp_file()
1122        self.__init_format_map__()
1123
1124    def __init_format_map__(self):
1125        """
1126        Initialize database format mapper, enabling multiple format requests.
1127
1128        """
1129
1130        # x is a pandas dataframe similar to self._data_frame format
1131        # Define format workflows
1132        self.format_map = {
1133            "msp": lambda x, normalize, fe_kwargs: self._to_msp(x, normalize),
1134            "flashentropy": lambda x, normalize, fe_kwargs: self._to_flashentropy(
1135                x, normalize, fe_kwargs
1136            ),
1137            "df": lambda x, normalize, fe_kwargs: self._to_df(x, normalize),
1138        }
1139
1140        # Add aliases
1141        self.format_map["fe"] = self.format_map["flashentropy"]
1142        self.format_map["flash-entropy"] = self.format_map["flashentropy"]
1143        self.format_map["dataframe"] = self.format_map["df"]
1144        self.format_map["data-frame"] = self.format_map["df"]
1145
1146    def _read_msp_file(self):
1147        """
1148        Reads the MSP files into the pandas dataframe, and sort/remove zero intensity ions in MS/MS spectra.
1149
1150        Returns
1151        -------
1152        :obj:`~pandas.DataFrame`
1153            DataFrame of spectra from the MSP file, exacly as it is in the file (no sorting, filtering etc)
1154        """
1155        # If input_dataframe is provided, return it it
1156        spectra = []
1157        spectrum = {}
1158
1159        f = StringIO(self._file_content)
1160        for line in f:
1161            line = line.strip()
1162            if not line:
1163                continue  # Skip empty lines
1164
1165            # Handle metadata
1166            if ":" in line:
1167                key, value = line.split(":", 1)
1168                key = key.strip().lower()
1169                value = value.strip()
1170
1171                if key == "name":
1172                    # Save current spectrum and start a new one
1173                    if spectrum:
1174                        spectra.append(spectrum)
1175                    spectrum = {"name": value, "peaks": []}
1176                else:
1177                    spectrum[key] = value
1178
1179            # Handle peak data (assumed to start with a number)
1180            elif line[0].isdigit():
1181                peaks = line.split()
1182                m_z = float(peaks[0])
1183                intensity = float(peaks[1])
1184                spectrum["peaks"].append(([m_z, intensity]))
1185        # Save the last spectrum
1186        if spectrum:
1187            spectra.append(spectrum)
1188
1189        df = pd.DataFrame(spectra)
1190        for column in df.columns:
1191            if column != "peaks":  # Skip 'peaks' column
1192                try:
1193                    df[column] = pd.to_numeric(df[column], errors="raise")
1194                except:
1195                    pass
1196        return df
1197
1198    def _to_df(self, input_dataframe, normalize=True):
1199        """
1200        Convert MSP-derived library to FlashEntropy library. 
1201
1202        Parameters
1203        ----------
1204        input_dataframe : :obj:`~pandas.DataFrame`
1205            Input DataFrame containing MSP-formatted spectra.
1206        normalize : bool, optional
1207            Normalize each spectrum by its magnitude.
1208            Default is True.
1209
1210        Returns
1211        -------
1212        :obj:`~pandas.DataFrame`
1213            DataFrame of with desired normalization
1214        """
1215        if not normalize:
1216            return input_dataframe
1217        else:
1218            # Convert to dictionary
1219            db_dict = input_dataframe.to_dict(orient="records")
1220
1221            # Initialize empty library
1222            lib = []
1223
1224            # Enumerate spectra
1225            for i, source in enumerate(db_dict):
1226                spectrum = source
1227                # Check that spectrum["peaks"] exists
1228                if "peaks" not in spectrum.keys():
1229                    raise KeyError(
1230                        "MSP not interpretted correctly, 'peaks' key not found in spectrum, check _dataframe attribute."
1231                    )
1232
1233                # Convert spectrum["peaks"] to numpy array
1234                if not isinstance(spectrum["peaks"], np.ndarray):
1235                    spectrum["peaks"] = np.array(spectrum["peaks"])
1236
1237                # Normalize peaks, if requested
1238                if normalize:
1239                    spectrum["peaks"] = self.normalize_peaks(spectrum["peaks"])
1240                    spectrum["num peaks"] = len(spectrum["peaks"])
1241
1242                # Add spectrum to library
1243                lib.append(spectrum)
1244            
1245            # Convert to DataFrame
1246            df = pd.DataFrame(lib)
1247            return df
1248    
1249    def _to_flashentropy(self, input_dataframe, normalize=True, fe_kwargs={}):
1250        """
1251        Convert MSP-derived library to FlashEntropy library.
1252
1253        Parameters
1254        ----------
1255        input_dataframe : :obj:`~pandas.DataFrame`
1256            Input DataFrame containing MSP-formatted spectra.
1257        normalize : bool
1258            Normalize each spectrum by its magnitude.
1259        fe_kwargs : dict, optional
1260            Keyword arguments for instantiation of FlashEntropy search and building index for FlashEntropy search;
1261            any keys not recognized will be ignored. By default, all parameters set to defaults.
1262
1263        Returns
1264        -------
1265        :obj:`~ms_entropy.FlashEntropySearch`
1266            MS2 library as FlashEntropy search instance.
1267
1268        Raises
1269        ------
1270        ValueError
1271            If "min_ms2_difference_in_da" or "max_ms2_tolerance_in_da" are present in `fe_kwargs` and they
1272        """
1273        self._check_flash_entropy_kwargs(fe_kwargs)
1274
1275        db_df = input_dataframe
1276
1277        # Convert to dictionary
1278        db_dict = db_df.to_dict(orient="records")
1279
1280        # Initialize empty library
1281        fe_lib = []
1282
1283        # Enumerate spectra
1284        for i, source in enumerate(db_dict):
1285            # Reorganize source dict, if necessary
1286            if "spectrum_data" in source.keys():
1287                spectrum = source["spectrum_data"]
1288            else:
1289                spectrum = source
1290
1291            # Rename precursor_mz key for FlashEntropy
1292            if "precursor_mz" not in spectrum.keys():
1293                if "precursormz" in spectrum:
1294                    spectrum["precursor_mz"] = spectrum.pop("precursormz")
1295                elif "precursor_ion" in spectrum:
1296                    spectrum["precursor_mz"] = spectrum.pop("precursor_ion")
1297                else:
1298                    raise KeyError(
1299                        "MSP must have either 'precursormz' or 'precursor_ion' key to be converted to FlashEntropy format."
1300                    )
1301
1302            # Check that spectrum["peaks"] exists
1303            if "peaks" not in spectrum.keys():
1304                raise KeyError(
1305                    "MSP not interpretted correctly, 'peaks' key not found in spectrum, check _dataframe attribute."
1306                )
1307
1308            # Convert spectrum["peaks"] to numpy array
1309            if not isinstance(spectrum["peaks"], np.ndarray):
1310                spectrum["peaks"] = np.array(spectrum["peaks"])
1311
1312            # Normalize peaks, if requested
1313            if normalize:
1314                spectrum["peaks"] = self.normalize_peaks(spectrum["peaks"])
1315
1316            # Add spectrum to library
1317            fe_lib.append(spectrum)
1318
1319        # Build FlashEntropy index
1320        fe_search = self._build_flash_entropy_index(fe_lib, fe_kwargs=fe_kwargs)
1321
1322        return fe_search
1323    
1324    def _check_msp_compatibility(self):
1325        """
1326        Check if the MSP file is compatible with the get_metabolomics_spectra_library method and provide feedback if it is not.
1327        """
1328        # Check polarity
1329        if (
1330            "polarity" not in self._data_frame.columns
1331            and "ionmode" not in self._data_frame.columns
1332        ):
1333            raise ValueError(
1334                "Neither 'polarity' nor 'ionmode' columns found in the input MSP metadata. Please check the file."
1335            )
1336        polarity_column = (
1337            "polarity" if "polarity" in self._data_frame.columns else "ionmode"
1338        )
1339
1340        # Check if polarity_column contents is either "positive" or "negative"
1341        if not all(self._data_frame[polarity_column].isin(["positive", "negative"])):
1342            raise ValueError(
1343                f"Input field on MSP '{polarity_column}' must contain only 'positive' or 'negative' values."
1344            )
1345
1346        # Check if the MSP file contains the required columns for metabolite metadata
1347        # inchikey, by name, not null
1348        # either formula or molecular_formula, not null
1349        if not all(self._data_frame["inchikey"].notnull()):
1350            raise ValueError(
1351                "Input field on MSP 'inchikey' must contain only non-null values."
1352            )
1353        if (
1354            "formula" not in self._data_frame.columns
1355            and "molecular_formula" not in self._data_frame.columns
1356        ):
1357            raise ValueError(
1358                "Input field on MSP must contain either 'formula' or 'molecular_formula' columns."
1359            )
1360        molecular_formula_column = (
1361            "formula" if "formula" in self._data_frame.columns else "molecular_formula"
1362        )
1363        if not all(self._data_frame[molecular_formula_column].notnull()):
1364            raise ValueError(
1365                f"Input field on MSP '{molecular_formula_column}' must contain only non-null values."
1366            )
1367
1368    def get_metabolomics_spectra_library(
1369        self,
1370        polarity,
1371        metabolite_metadata_mapping={},
1372        format="fe",
1373        normalize=True,
1374        fe_kwargs={},
1375    ):
1376        """
1377        Prepare metabolomics spectra library and associated metabolite metadata
1378
1379        Note: this uses the inchikey as the index for the metabolite metadata dataframe and for connecting to the spectra, so it must be in the input
1380
1381        """
1382        # Check if the MSP file is compatible with the get_metabolomics_spectra_library method
1383        self._check_msp_compatibility()
1384
1385        # Check if the polarity parameter is valid and if a polarity column exists in the dataframe
1386        if polarity not in ["positive", "negative"]:
1387            raise ValueError("Polarity must be 'positive' or 'negative'")
1388        polarity_column = (
1389            "polarity" if "polarity" in self._data_frame.columns else "ionmode"
1390        )
1391
1392        # Get a subset of the initial dataframea by polarity
1393        db_df = self._data_frame[self._data_frame[polarity_column] == polarity].copy()
1394
1395        # Rename the columns of the db_df to match the MetaboliteMetadata dataclass using the metabolite_metadata_mapping
1396        # If the mapping is not provided, use the default mapping
1397        if not metabolite_metadata_mapping:
1398            metabolite_metadata_mapping = {
1399                "chebi_id": "chebi",
1400                "kegg_id": "kegg",
1401                "refmet_name": "common_name",
1402                "molecular_formula": "formula",
1403                "gnps_spectra_id":"id",
1404                "precursormz": "precursor_mz",
1405                "precursortype":"ion_type"
1406            }
1407        db_df.rename(columns=metabolite_metadata_mapping, inplace=True)
1408        db_df["molecular_data_id"] = db_df["inchikey"]
1409
1410
1411
1412        # Check if the resulting dataframe has the required columns for the flash entropy search
1413        required_columns = ["molecular_data_id", "precursor_mz", "ion_type", "id"]
1414        for col in required_columns:
1415            if col not in db_df.columns:
1416                raise ValueError(
1417                    f"Input field on MSP must contain '{col}' column for FlashEntropy search."
1418                )
1419
1420        # Pull out the metabolite metadata from the dataframe and put it into a different dataframe
1421        # First get a list of the possible attributes of the MetaboliteMetadata dataclass
1422        metabolite_metadata_keys = list(MetaboliteMetadata.__annotations__.keys())
1423        # Replace id with molecular_data_id in metabolite_metadata_keys
1424        metabolite_metadata_keys = [
1425            "molecular_data_id" if x == "id" else x for x in metabolite_metadata_keys
1426        ]
1427        metabolite_metadata_df = db_df[
1428            db_df.columns[db_df.columns.isin(metabolite_metadata_keys)]
1429        ].copy()
1430
1431        # Make unique and recast the id column for metabolite metadata
1432        metabolite_metadata_df.drop_duplicates(subset=["molecular_data_id"], inplace=True)
1433        metabolite_metadata_df["id"] = metabolite_metadata_df["molecular_data_id"]
1434
1435        # Convert to a dictionary using the inchikey as the key
1436        metabolite_metadata_dict = metabolite_metadata_df.to_dict(
1437            orient="records"
1438        )
1439        metabolite_metadata_dict = {
1440            v["id"]: self._dict_to_dataclass(v, MetaboliteMetadata)
1441            for v in metabolite_metadata_dict
1442        }
1443
1444        # Remove the metabolite metadata columns from the original dataframe
1445        for key in metabolite_metadata_keys:
1446            if key != "molecular_data_id":
1447                if key in db_df.columns:
1448                    db_df.drop(columns=key, inplace=True)
1449
1450        # Format the spectral library
1451        format_func = self._get_format_func(format)
1452        lib = format_func(db_df, normalize=normalize, fe_kwargs=fe_kwargs)
1453        return (lib, metabolite_metadata_dict)