Integration_with_mzMatch.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "f74772f3",
   "metadata": {},
   "source": [
    "# Integration with mzMatch and PeakMLViewerPy\n",
    "\n",
    "The new ipaPy2 package is fully integrated with [mzMatch](https://github.com/UoMMIB/mzmatch.R).\n",
    "mzMatch is a data processing pipeline for metabolomics LC/MS data. This pipeline revolves around the PeakML file format, which is designed to contain all the relevant information obtained from the processing of an LC/MS-based untargeted metabolomics experiment.  \n",
    "\n",
    "As shown below, the new IPA implementation includes the necessary functions to read a .peakml file and extracts the information required to run the IPA method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "a201bd87",
   "metadata": {},
   "outputs": [],
   "source": [
    "from ipaPy2 import PeakMLIO\n",
    "from ipaPy2 import ipa\n",
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "529ad22a",
   "metadata": {},
   "source": [
    "All the information about the ReadPeakML() function can be found in the help:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "9f100106",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on function ReadPeakML in module ipaPy2.PeakMLIO:\n",
      "\n",
      "ReadPeakML(filename)\n",
      "    Loading data from a PeakML file\n",
      "    \n",
      "    Parameters\n",
      "    ----------\n",
      "    filename : string with the name (including path) of the .peakml file.\n",
      "                It is important to use the allpeaks file generated by the\n",
      "                mzmatch anlysis, for the proper use of the IPA method.\n",
      "    \n",
      "    Returns\n",
      "    -------\n",
      "    df : pandas dataframe (necessary)\n",
      "         A dataframe containing the MS1 data including the following columns:\n",
      "            -ids: an unique id for each feature\n",
      "            -rel.ids:   relation ids. In a previous step of the data processing\n",
      "                        pipeline, features are clustered based on peak shape\n",
      "                        similarity/retention time. Features in the same\n",
      "                        cluster are likely to come from the same metabolite.\n",
      "                        All isotope patterns must be in the same rel.id\n",
      "                        cluster.\n",
      "            -mzs: mass-to-charge ratios, usually the average across\n",
      "                  different samples.\n",
      "            -RTs: retention times in seconds, usually the average across\n",
      "                  different samples.\n",
      "            -Ints: representative (e.g., maximum or average) intensity detected\n",
      "                   for each feature across samples (either peak area or peak\n",
      "                   intensity)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(PeakMLIO.ReadPeakML)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00edd3e7",
   "metadata": {},
   "source": [
    "Using this function, one can read a .peakml file and obtain the dataframe needed to run the IPA pipeline.\n",
    "For this example, we used a rather small example .peakml file that can be downloaded from [here](https://drive.google.com/file/d/123bDs8kMTlDbjd1gETSAEe6uCW_fgsSW/view?usp=share_link)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "9a3e032d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "loading Example_allpeaks.peakml...\n",
      "parsing 959 peaks...\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>ids</th>\n",
       "      <th>rel.ids</th>\n",
       "      <th>mzs</th>\n",
       "      <th>RTs</th>\n",
       "      <th>Int</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "      <td>116.070547</td>\n",
       "      <td>45.770423</td>\n",
       "      <td>2.170017e+09</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>88</td>\n",
       "      <td>0</td>\n",
       "      <td>117.073691</td>\n",
       "      <td>45.787586</td>\n",
       "      <td>1.256520e+08</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>372</td>\n",
       "      <td>0</td>\n",
       "      <td>70.065197</td>\n",
       "      <td>45.795901</td>\n",
       "      <td>3.487656e+07</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>501</td>\n",
       "      <td>0</td>\n",
       "      <td>231.133686</td>\n",
       "      <td>46.183948</td>\n",
       "      <td>2.519223e+07</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>104.106842</td>\n",
       "      <td>40.843309</td>\n",
       "      <td>1.889172e+09</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   ids  rel.ids         mzs        RTs           Int\n",
       "0    1        0  116.070547  45.770423  2.170017e+09\n",
       "1   88        0  117.073691  45.787586  1.256520e+08\n",
       "2  372        0   70.065197  45.795901  3.487656e+07\n",
       "3  501        0  231.133686  46.183948  2.519223e+07\n",
       "4    2        1  104.106842  40.843309  1.889172e+09"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = PeakMLIO.ReadPeakML('Example_allpeaks.peakml')\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "19a29601",
   "metadata": {},
   "source": [
    "In order to run the IPA method on this dataset, it is necessary to load the necessary databases:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "17c0d799",
   "metadata": {},
   "outputs": [],
   "source": [
    "DB=pd.read_csv('DB/IPA_MS1.csv')\n",
    "adducts = pd.read_csv('DB/adducts.csv')\n",
    "Bio = pd.read_csv('DB/allBIO_reactions.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d9d424d3",
   "metadata": {},
   "source": [
    "Finally, we can run the whole pipeline with the simpleIPA() function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "bb74007d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mapping isotope patterns ....\n",
      "1.1 seconds elapsed\n",
      "computing all adducts - Parallelized ....\n",
      "23.7 seconds elapsed\n",
      "annotating based on MS1 information - Parallelized ...\n",
      "23.9 seconds elapsed\n",
      "computing posterior probabilities including biochemical and adducts connections\n",
      "initialising sampler ...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Gibbs Sampler Progress Bar: 100%|█████████| 3000/3000 [1:09:33<00:00,  1.39s/it]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "parsing results ...\n",
      "Done -  4175.6 seconds elapsed\n"
     ]
    }
   ],
   "source": [
    "annotations = ipa.simpleIPA(df=df,ionisation=1,DB=DB,adductsAll=adducts,ppm=3,ppmthr=5,Bio=Bio,\n",
    "                            delta_add=0.1,delta_bio=0.1,noits=3000,ncores=70)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b7b23490",
   "metadata": {},
   "source": [
    "As an example, below are shown the results for the annotation of the mass spectrometry feature associated with id=9"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "a76cf990",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>id</th>\n",
       "      <th>name</th>\n",
       "      <th>formula</th>\n",
       "      <th>adduct</th>\n",
       "      <th>m/z</th>\n",
       "      <th>charge</th>\n",
       "      <th>RT range</th>\n",
       "      <th>ppm</th>\n",
       "      <th>isotope pattern score</th>\n",
       "      <th>fragmentation pattern score</th>\n",
       "      <th>prior</th>\n",
       "      <th>post</th>\n",
       "      <th>post Gibbs</th>\n",
       "      <th>chi-square pval</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>C08307</td>\n",
       "      <td>Hordatine A</td>\n",
       "      <td>C28H40N8O4</td>\n",
       "      <td>M+2H</td>\n",
       "      <td>276.158077</td>\n",
       "      <td>2</td>\n",
       "      <td>None</td>\n",
       "      <td>-0.584944</td>\n",
       "      <td>4.877535e-01</td>\n",
       "      <td>None</td>\n",
       "      <td>0.486508</td>\n",
       "      <td>6.540751e-01</td>\n",
       "      <td>0.791852</td>\n",
       "      <td>1.262275e-50</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>NPA028941</td>\n",
       "      <td>Microginin 550</td>\n",
       "      <td>C27H44N4O8</td>\n",
       "      <td>M+2H</td>\n",
       "      <td>276.157408</td>\n",
       "      <td>2</td>\n",
       "      <td>None</td>\n",
       "      <td>1.836519</td>\n",
       "      <td>5.033153e-01</td>\n",
       "      <td>None</td>\n",
       "      <td>0.248086</td>\n",
       "      <td>3.441760e-01</td>\n",
       "      <td>0.208148</td>\n",
       "      <td>1.262275e-50</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Unknown</td>\n",
       "      <td>Unknown</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>3.000000</td>\n",
       "      <td>8.931157e-03</td>\n",
       "      <td>None</td>\n",
       "      <td>0.071043</td>\n",
       "      <td>1.748911e-03</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.262275e-50</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>C21776</td>\n",
       "      <td>2-({[(4-Methoxyphenyl)methyl](methyl)amino}met...</td>\n",
       "      <td>C14H23NNaO3</td>\n",
       "      <td>M+Na</td>\n",
       "      <td>276.157012</td>\n",
       "      <td>1</td>\n",
       "      <td>None</td>\n",
       "      <td>3.272514</td>\n",
       "      <td>8.931157e-09</td>\n",
       "      <td>None</td>\n",
       "      <td>0.048591</td>\n",
       "      <td>1.196186e-09</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.262275e-50</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>NPA001615</td>\n",
       "      <td>(Z)-N-(4-decenoyl)-L-homoserine lactone</td>\n",
       "      <td>C14H23NNaO3</td>\n",
       "      <td>M+Na</td>\n",
       "      <td>276.157012</td>\n",
       "      <td>1</td>\n",
       "      <td>None</td>\n",
       "      <td>3.272514</td>\n",
       "      <td>8.931157e-09</td>\n",
       "      <td>None</td>\n",
       "      <td>0.048591</td>\n",
       "      <td>1.196186e-09</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.262275e-50</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>NPA031672</td>\n",
       "      <td>Laricinin A</td>\n",
       "      <td>C14H23NNaO3</td>\n",
       "      <td>M+Na</td>\n",
       "      <td>276.157012</td>\n",
       "      <td>1</td>\n",
       "      <td>None</td>\n",
       "      <td>3.272514</td>\n",
       "      <td>8.931157e-09</td>\n",
       "      <td>None</td>\n",
       "      <td>0.048591</td>\n",
       "      <td>1.196186e-09</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.262275e-50</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>ET28010x_1</td>\n",
       "      <td>MPL-dm</td>\n",
       "      <td>C14H23NNaO3</td>\n",
       "      <td>M+Na</td>\n",
       "      <td>276.157012</td>\n",
       "      <td>1</td>\n",
       "      <td>None</td>\n",
       "      <td>3.272514</td>\n",
       "      <td>8.931157e-09</td>\n",
       "      <td>None</td>\n",
       "      <td>0.048591</td>\n",
       "      <td>1.196186e-09</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>1.262275e-50</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "           id                                               name      formula  \\\n",
       "0      C08307                                        Hordatine A   C28H40N8O4   \n",
       "1   NPA028941                                     Microginin 550   C27H44N4O8   \n",
       "2     Unknown                                            Unknown         None   \n",
       "3      C21776  2-({[(4-Methoxyphenyl)methyl](methyl)amino}met...  C14H23NNaO3   \n",
       "4   NPA001615            (Z)-N-(4-decenoyl)-L-homoserine lactone  C14H23NNaO3   \n",
       "5   NPA031672                                        Laricinin A  C14H23NNaO3   \n",
       "6  ET28010x_1                                             MPL-dm  C14H23NNaO3   \n",
       "\n",
       "  adduct         m/z charge RT range       ppm  isotope pattern score  \\\n",
       "0   M+2H  276.158077      2     None -0.584944           4.877535e-01   \n",
       "1   M+2H  276.157408      2     None  1.836519           5.033153e-01   \n",
       "2   None        None   None     None  3.000000           8.931157e-03   \n",
       "3   M+Na  276.157012      1     None  3.272514           8.931157e-09   \n",
       "4   M+Na  276.157012      1     None  3.272514           8.931157e-09   \n",
       "5   M+Na  276.157012      1     None  3.272514           8.931157e-09   \n",
       "6   M+Na  276.157012      1     None  3.272514           8.931157e-09   \n",
       "\n",
       "  fragmentation pattern score     prior          post  post Gibbs  \\\n",
       "0                        None  0.486508  6.540751e-01    0.791852   \n",
       "1                        None  0.248086  3.441760e-01    0.208148   \n",
       "2                        None  0.071043  1.748911e-03    0.000000   \n",
       "3                        None  0.048591  1.196186e-09    0.000000   \n",
       "4                        None  0.048591  1.196186e-09    0.000000   \n",
       "5                        None  0.048591  1.196186e-09    0.000000   \n",
       "6                        None  0.048591  1.196186e-09    0.000000   \n",
       "\n",
       "   chi-square pval  \n",
       "0     1.262275e-50  \n",
       "1     1.262275e-50  \n",
       "2     1.262275e-50  \n",
       "3     1.262275e-50  \n",
       "4     1.262275e-50  \n",
       "5     1.262275e-50  \n",
       "6     1.262275e-50  "
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "annotations[9]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "acd382de",
   "metadata": {},
   "source": [
    "The IPA annotations obtained can be added directly to a .peakml file via the add_IPA_to_PeakML() function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "b623c176",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on function add_IPA_to_PeakML in module ipaPy2.PeakMLIO:\n",
      "\n",
      "add_IPA_to_PeakML(file, IPA_Data, out_File)\n",
      "    Adding IPA annotation to PeakML file\n",
      "    \n",
      "    Parameters\n",
      "    ----------\n",
      "    file : string with the name (including path) of the .peakml file.\n",
      "    IPA_Data : Dictionary containing the IPA annotation\n",
      "    out_File : string with the name (including path) where the annotated\n",
      "                .peakml file will be saved.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(PeakMLIO.add_IPA_to_PeakML)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "c2d2760a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Adding IPA annotation to Example_allpeaks.peakml and saving it as Example_allpeaks_IPA_annotated.peakml\n"
     ]
    }
   ],
   "source": [
    "PeakMLIO.add_IPA_to_PeakML(\"Example_allpeaks.peakml\",annotations,\"Example_allpeaks_IPA_annotated.peakml\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f73dd795",
   "metadata": {},
   "source": [
    "The dictionary containing all annotations can be also saved in a .pickle file."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d14e93e8",
   "metadata": {},
   "source": [
    "This file can be read by the [PeakMLViewerpy](https://github.com/UoMMIB/PeakMLViewerPy), allowing the exploration of the .peakml file together with the IPA annotation.\n",
    "PeakMLViewerpy is availabe on GitHub [here](https://github.com/UoMMIB/PeakMLViewerPy) together with a detailed installation guide."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4946e9f5",
   "metadata": {},
   "source": [
    "![screenshot](PeakMLViewerpy_example.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "727e2637",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.10.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}