06_Spectral_Matching.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import glob\n",
    "import shutil\n",
    "import os\n",
    "from pyteomics import mgf, auxiliary\n",
    "from pyopenms import *\n",
    "import pyteomics\n",
    "from pyteomics import mztab\n",
    "from src.export_feature_matrix import export"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "path = os.path.join(\"results\", \"interim\", \"annotations\")\n",
    "\n",
    "if os.path.exists(path):\n",
    "    shutil.rmtree(path)\n",
    "\n",
    "os.mkdir(path)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First,import the feature table:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DF_features = pd.read_csv(os.path.join(\"results\", \"interim\", \"FeatureMatrix.tsv\"), sep=\"\\t\")\n",
    "DF_features.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### `2) Spectral matching`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the MGF file to an MSExperiment format:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mgf_file= os.path.join(\"results\", \"GNPSexport\", \"MSMS.mgf\")\n",
    "exp = MSExperiment()\n",
    "MascotGenericFile().load(mgf_file, exp)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Perform spectral matching with a library in MGF format that is located under \"resources\":\n",
    "(here we are using the one freely available from GNPS (all libraries))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "database = os.path.join(\"resources\", \"GNPS-LIBRARY.mgf\")\n",
    "speclib = MSExperiment()\n",
    "MascotGenericFile().load(database, speclib)\n",
    "mztab = MzTab()\n",
    "output_mztab = os.path.join(\"results\", \"interim\", \"annotations\", \"MSMS.mzTab\")\n",
    "out_merged = \"\"\n",
    "MSMS_match = MetaboliteSpectralMatching()\n",
    "MSMS_match_par = MSMS_match.getDefaults()\n",
    "MSMS_match_par.setValue('merge_spectra', 'false')\n",
    "MSMS_match.setParameters(MSMS_match_par)\n",
    "MSMS_match.run(exp, speclib, mztab,  String(out_merged))\n",
    "MzTabFile().store(output_mztab, mztab)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clean up the mzTab to a dataframe:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spectralmatch =  pyteomics.mztab.MzTab(output_mztab, encoding=\"UTF8\", table_format=\"df\")\n",
    "spectralmatch.metadata\n",
    "df = spectralmatch.small_molecule_table\n",
    "spectralmatch_DF = df.drop(columns= [\"identifier\", \"inchi_key\", \"modifications\", \"calc_mass_to_charge\", \"opt_adduct_ion\", \"taxid\", \"species\", \"database\", \"spectra_ref\", \"search_engine\", \"opt_sec_id\",\"smallmolecule_abundance_std_error_study_variable[1]\", \"smallmolecule_abundance_stdev_study_variable[1]\", \"smallmolecule_abundance_study_variable[1]\", \"chemical_formula\"])\n",
    "spectralmatch_DF =spectralmatch_DF[spectralmatch_DF[\"opt_ppm_error\"] <= 10] \n",
    "spectralmatch_DF =spectralmatch_DF[spectralmatch_DF[\"opt_ppm_error\"] >= -10]\n",
    "spectralmatch_DF =spectralmatch_DF[spectralmatch_DF[\"opt_match_score\"] >= 60]\n",
    "spectralmatch_DF[\"opt_spec_native_id\"]= spectralmatch_DF[\"opt_spec_native_id\"].str.replace(r\"index=\", \"\")\n",
    "spectralmatch_DF       "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Annotate the feature matrix with the spectral matches:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "path = os.path.join(\"results\", \"GNPSexport\", \"MSMS.mgf\")\n",
    "file= mgf.MGF(source=path, use_header=True, convert_arrays=2, read_charges=True, read_ions=False, dtype=None, encoding=None)\n",
    "parameters=[]\n",
    "for spectrum in file:\n",
    "    parameters.append(spectrum['params'])\n",
    "mgf_file= pd.DataFrame(parameters)\n",
    "mgf_file[\"feature_id\"] = mgf_file[\"feature_id\"].str.replace(r\"e_\", \"\")\n",
    "\n",
    "df = exp.get_df() #mzML loaded in MSExperiment()\n",
    "for spec in exp:\n",
    "    df[\"index\"] = [spec.getNativeID() for spec in exp]\n",
    "    df[\"SCANS\"] = [spec.getMetaValue(\"Scan_ID\") for spec in exp]\n",
    "df[\"index\"] = df[\"index\"].str.replace(r\"index=\", \"\")\n",
    "\n",
    "# Add Scan numbers to spectral match DF\n",
    "spectralmatch_DF[\"SCANS\"] = \"\"\n",
    "for i, idx in zip(spectralmatch_DF.index, spectralmatch_DF[\"opt_spec_native_id\"]):\n",
    "    hits = []\n",
    "    for index, scan_number, in zip(df[\"index\"], df[\"SCANS\"]):\n",
    "        if idx==index:\n",
    "            hit = f\"{scan_number}\"\n",
    "            if hit not in hits:\n",
    "                hits.append(hit)\n",
    "    spectralmatch_DF.loc[i, \"SCANS\"] = \" ## \".join(hits)\n",
    "\n",
    "# Add Scan numbers to feature DF\n",
    "scans = []\n",
    "for consensus_id in DF_features[\"id\"].astype(str):\n",
    "    hits = []\n",
    "    for scan, mgf_id in zip(mgf_file[\"scans\"], mgf_file[\"feature_id\"]):\n",
    "        if consensus_id == mgf_id:\n",
    "            hit = f\"{scan}\"\n",
    "            if hit not in hits:\n",
    "                hits.append(hit)\n",
    "    scans.append(\" ## \".join(hits))\n",
    "\n",
    "DF_features[\"SCANS\"] = scans\n",
    "\n",
    "DF_features[\"SpectralMatch\"] = \"\"\n",
    "DF_features[\"SpectralMatch_smiles\"] = \"\"\n",
    "\n",
    "for i, scan in zip(DF_features.index, DF_features[\"SCANS\"]):\n",
    "    hits1 = []\n",
    "    hits2=[]\n",
    "    for name, smiles, scan_number, in zip(spectralmatch_DF[\"description\"],spectralmatch_DF[\"smiles\"], spectralmatch_DF[\"SCANS\"]):\n",
    "        if scan==scan_number:\n",
    "            hit1 = f\"{name}\"\n",
    "            hit2 = f\"{smiles}\"\n",
    "            if hit1 not in hits1:\n",
    "                hits1.append(hit1)\n",
    "                hits2.append(hit2)\n",
    "    DF_features.loc[i, \"SpectralMatch\"] = \" ## \".join(hits1)\n",
    "    DF_features.loc[i, \"SpectralMatch_smiles\"] = \" ## \".join(hits2)\n",
    "\n",
    "export(DF_features)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.0"
  },
  "vscode": {
   "interpreter": {
    "hash": "edde62aa2661007f0756e9790e7a328c288a583bf6ce768a355147dac67c8db8"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}