OpenMS_SOP.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pyopenms\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1>Introduction</h1>\n",
    "PyOpenMS is a set of Python bindings for the OpenMS C++ library. The included functions are powerful but can be daunting compared to native Python libraries. This SOP will serve as a guide to PyOpenMS concepts and include examples of how to do some routine data processing. \n",
    "\n",
    "<h1>File Processing</h1>\n",
    "\n",
    "We'll start by opening an MzML file and viewing its contents. The first step is to import the MzML module as well as the classes that describe experiments and spectra."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyopenms import MzMLFile, MSExperiment, MSSpectrum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll instantiate an MzMLFile object and an MSExperiment obejct, and load our file into the experiment. In this example we'll open the file P3_5+_3.mzML."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "file = MzMLFile()  # Instantiate an MzMLFile object.\n",
    "exp = MSExperiment()  # Instantiate an empty MSExperiment object.\n",
    "\n",
    "file.load(\"P3_5+_3.mzML\", exp)  # Load the data from P3_5+_3.mzML into the newly instantiated experiment."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In OpenMS terminology, an \"experiment\" is a collection of spectra, such as an LC-MS dataset or an infusion experiment with a lot of averages. The experiment, called exp here, now contains the spectral data from P3_5+\\_3.mzml. Individual spectra can be accessed by their index within the experiment like this: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<pyopenms.pyopenms.MSSpectrum object at 0x7f3571ce1b70>\n"
     ]
    }
   ],
   "source": [
    "spectrum = exp[0]  # Access the first element (i.e. scan) of the experiment.\n",
    "\n",
    "print(spectrum)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see, accessing the first element of our experiment returns a MSSpectrum object. This object contains a lot of information about the spectrum. If you're planning on doing any sort of spectral processing, the most important method is <code>get_peaks()</code>, which returns a tuple of arrays containing the raw mass and intensity values for each point in the spectrum:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[  204.99993902   205.00009595   205.00025288 ...,  2002.37720191\n",
      "  2002.38199264  2002.38678339]\n",
      "[ 0.  0.  0. ...,  0.  0.  0.]\n"
     ]
    }
   ],
   "source": [
    "peaks = spectrum.get_peaks()\n",
    "\n",
    "print(peaks[0])  # This is the array of masses.\n",
    "print(peaks[1])  # This is the array of intensity values\n",
    "\n",
    "q = spectrum.getInstrumentSettings()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In most cases, you'll probably be more interested in working with a centroided spectrum rather than the raw instrument data. OpenMS contains some powerful algorithms for signal processing, including a function specifically designed to pick peaks in FTMS data. In this next section we'll generate a profile-mode spectrum from our spectrum object and highlight how OpenMS method parameters can be configured. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{b'SignalToNoise:auto_max_percentile': 95,\n",
       " b'SignalToNoise:auto_max_stdev_factor': 3.0,\n",
       " b'SignalToNoise:auto_mode': 0,\n",
       " b'SignalToNoise:bin_count': 30,\n",
       " b'SignalToNoise:max_intensity': -1,\n",
       " b'SignalToNoise:min_required_elements': 10,\n",
       " b'SignalToNoise:noise_for_empty_window': 1e+20,\n",
       " b'SignalToNoise:win_len': 200.0,\n",
       " b'SignalToNoise:write_log_messages': b'true',\n",
       " b'missing': 1,\n",
       " b'ms_levels': [1],\n",
       " b'report_FWHM': b'false',\n",
       " b'report_FWHM_unit': b'relative',\n",
       " b'signal_to_noise': 0.0,\n",
       " b'spacing_difference': 1.5,\n",
       " b'spacing_difference_gap': 4.0}"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from pyopenms import PeakPickerHiRes\n",
    "\n",
    "peakPicker = PeakPickerHiRes()\n",
    "\n",
    "peakPicker.getParameters().asDict()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We've instantiated a PeakPickerHiRes object and printed out its parameters as a dictionary. It looks like right now the peak picker is set to only pick MS1 spectra, since the <code>ms_levels</code> parameter is set to <code>[1]</code>. Let's change it so that it will process both MS1 and MS2 spectra. The first step is going to be to create a new parameter object.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyopenms import Param\n",
    "\n",
    "peakPickerParams = Param()\n",
    "\n",
    "peakPickerParams.setValue(b'ms_levels', [1, 2], '')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notice that we've passed an empty string as the third argument in order to make the debugger happy. We also have to pass the parameter name, 'ms_levels,' as a bytestring by putting the <code>b</code> in front. Now that we have our param object, let's apply it to <code>peakPicker</code>:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1, 2]"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "peakPicker.setParameters(peakPickerParams)\n",
    "\n",
    "peakPicker.getParameters().asDict()[b'ms_levels']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now our peak picker is properly configured to process MS1 and MS2 spectra. This basic pattern of instantiating a parameter object and applying it to a method is how settings are changed for any method in the pyOpenMS library. \n",
    "\n",
    "The next step is to actually apply our peak picker in order to process the spectrum. There are two ways to do this: we could either process individual spectra, or apply the algorithm to every scan in the experiment. First let's process the <code>spectrum</code> we defined earlier. The peak picker works by processing one <code>MSSpectrum</code> object and writing to a second <code>MSSpectrum</code>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The profile spectrum contains 198120 data points.\n",
      "Our peak picker found 15500 peaks.\n"
     ]
    }
   ],
   "source": [
    "centroidSpectrum = MSSpectrum()\n",
    "\n",
    "peakPicker.pick(spectrum, centroidSpectrum) \n",
    "\n",
    "peaksinProfileSpectrum = spectrum.get_peaks()[0].size\n",
    "peaksInCentroidSpectrum = centroidSpectrum.get_peaks()[0].size\n",
    "\n",
    "print(\"The profile spectrum contains %d data points.\"%peaksinProfileSpectrum)\n",
    "\n",
    "print(\"Our peak picker found %d peaks.\"%peaksInCentroidSpectrum)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If your experiment contains multiple spectra, you can also apply the peakPicker to the entire experiment. Again, we create a new <code>MSExperiment</code> object to write our peaks into:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The profile spectrum contains 198120 data points.\n",
      "Our peak picker found 15500 peaks.\n"
     ]
    }
   ],
   "source": [
    "centroidExperiment = MSExperiment()\n",
    "\n",
    "peakPicker.pickExperiment(exp, centroidExperiment)\n",
    "\n",
    "# Here we're looking at the number of data poins in the first spectrum to make sure our peak picking worked. \n",
    "peaksinProfileExperiment = exp[0].get_peaks()[0].size\n",
    "peaksInCentroidExperiment = centroidExperiment[0].get_peaks()[0].size\n",
    "\n",
    "print(\"The profile spectrum contains %d data points.\"%peaksinProfileExperiment)\n",
    "\n",
    "print(\"Our peak picker found %d peaks.\"%peaksInCentroidExperiment)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The last basic operation on spectrum objects that I'm going to cover is how to write a processed spectra to a file. This is useful if you want to open up your spectrum in a MzML viewer such as TOPP or mMass. To write a file, we instantiate a new <code>MzMLFile</code> object and pass an <code>MSExperiment</code> to its <code>store</code> method. Alternatively, if we want to write an individual spectrum, we need to create a new <code>MSExperiment</code> and add our spectrum to it before writing to the file. The following chunk of code demonstrates both methods.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Method 1: Storing an experiment to a file\n",
    "fileToExport = MzMLFile()\n",
    "fileToExport.store('ProcessedExperiment.mzML', centroidExperiment)\n",
    "\n",
    "# Method 2: Adding a spectrum to a new experiment, then writing to a file. \n",
    "experimentToWrite = MSExperiment()\n",
    "experimentToWrite.addSpectrum(centroidSpectrum) \n",
    "fileToExport2 = MzMLFile()\n",
    "fileToExport2.store('ProcessedSpectrum.mzML', experimentToWrite)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you've read this far, you now know how to open a mass spectral file using pyOpenMS and perform a peak-picking analysis on both individual spectra and the full experiment. The library contains a lot of other useful processing methods, but the workflow for using them is going to be very similar: instantiate the method, provide it with a parameter object if necessary, and run the analysis on your data.\n",
    "\n",
    "<h1>The Chemistry Module</h1>\n",
    "\n",
    "The next section will cover some useful features of PyOpenMS's <code>Chemistry</code> module. The chemistry module contains a lot of neat methods that let you manipulate amino acid sequences to compute their masses, theoretical fragmentation, and isotope distributions. The following chunk of code shows how to create a sequence object for the peptide MASSSPECISNEAT and get some useful information out of it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyopenms import AASequence\n",
    "\n",
    "# generate a new sequence object from a string containing the sequence. \n",
    "myPeptideSequence = AASequence().fromString(\"MASSSPECISNEAT\", True)\n",
    "\n",
    "# Get the monoisotopic weight (for average weight, use getAverageWeight instead of getMonoWeight)\n",
    "monoIsotopicWeightUncharged = myPeptideSequence.getMonoWeight(0, 0)\n",
    "monoIsotopicWeight_3plus = myPeptideSequence.getMonoWeight(0, 3)\n",
    "\n",
    "# Compute the molecular composition\n",
    "formula_unprotonated = myPeptideSequence.getFormula(0, 0).toString()\n",
    "\n",
    "# get the first 5 residues- positional arguments are start and endpoints\n",
    "first5Residues = myPeptideSequence.getSubsequence(0, 5).toString()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You might have noticed that the functions to compute the weight and composition of the peptide take two positional arguments- e.g. <code>getFormula(0, 0)</code>. The first argument is a ResidueType code and determines whether the sequence should be treated as an intact peptide or as a fragment, according to the following lookup table:\n",
    "<table>\n",
    "    <tr><th>Key</th><th>Label</th><th>Description (From OpenMS documentation)</th></tr>\n",
    "    <tr><th>0</th><th>Full</th><th>With N-terminus and C-terminus</th></tr>\n",
    "    <tr><th>1</th><th>Internal</th><th>Without Termini</th></tr>\n",
    "    <tr><th>2</th><th>N-terminal</th><th>Only N-terminus</th></tr>\n",
    "    <tr><th>3</th><th>C-terminal</th><th>Only C-terminus</th></tr>\n",
    "    <tr><th>4</th><th>A-Ion</th><th>N-terminus up to C-alpha/carbonyl carbon bond</th></tr>\n",
    "    <tr><th>5</th><th>B-Ion</th><th>N-terminus up to the peptide bond</th></tr>\n",
    "    <tr><th>6</th><th>C-Ion</th><th>N-terminus up to the amide/c-alpha bond</th></tr>\n",
    "    <tr><th>7</th><th>X-Ion</th><th>Amide/c-alpha bond up to the C-terminus</tr>\n",
    "    <tr><th>8</th><th>Y-Ion</th><th>Peptide bond up to the C-terminus</tr>\n",
    "    <tr><th>9</th><th>Z-Ion</th><th>C-alpha/carbonyl carbon bond up to C-terminus</tr>\n",
    "</table>\n",
    "\n",
    "The second argument is the charge, and determines how many protons should be added to the sequence before computation. As a concrete example, calling <code>myPeptideaSequence.getMonoWeight(4, 3)</code> would return the monoisotopic mass of the a<sub>14</sub><sup>3+</sup> ion for a peptide with MASSSPECISNEAT as the first 14 residues.\n",
    "\n",
    "<h2>Modification handling</h2>\n",
    "\n",
    "PyOpenMS has Unimod integration, and you can add modifications using the SetModification command. Unfortunately modifications that aren't in Unimod can't be added in this way."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The composition of the unoxidized peptide is b'C55H91N15O25S2'\n",
      "The composition of the oxidized peptide is b'C55H91N15O26S2'\n"
     ]
    }
   ],
   "source": [
    "myPeptideSequence.setModification(0, \"Oxidation\")\n",
    "\n",
    "formula_oxidized = myPeptideSequence.getFormula(0, 0).toString()\n",
    "\n",
    "print(\"The composition of the unoxidized peptide is %s\"%formula_unprotonated)\n",
    "\n",
    "print(\"The composition of the oxidized peptide is %s\"%formula_oxidized)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can see that setting the odixation caused the formula to change to reflect the additional oxygen. \n",
    "\n",
    "PyOpenMS contains functions to compute the isotopic distributins of sequences and compositions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(1441, 0.4383429723511864),\n",
       " (1442, 0.30098864515812396),\n",
       " (1443, 0.164891991080421),\n",
       " (1444, 0.06560323788654057),\n",
       " (1445, 0.021944548182978953),\n",
       " (1446, 0.006238062479179644),\n",
       " (1447, 0.0015597236074226704),\n",
       " (1448, 0.0003477242543181491),\n",
       " (1449, 7.01571171361993e-05),\n",
       " (1450, 1.293788269227358e-05)]"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from pyopenms import EmpiricalFormula\n",
    "\n",
    "# compute the first 10 isotopes\n",
    "myPeptide_isotopes = myPeptideSequence.getFormula(0, 0).getIsotopeDistribution(10)\n",
    "\n",
    "# this getContainer function returns the isotope distribution as tuples of formula mass and abundance. \n",
    "myPeptide_isotopes.getContainer()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The EmpiricalFormula class can let you calculate isotopes distributions for arbitrary formulas. In the next example we'll add a bromine to the formula and compute the new isotope distribution. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(1520, 0.2222051573510236),\n",
       " (1521, 0.1525773959589096),\n",
       " (1522, 0.29974295165602416),\n",
       " (1523, 0.18167922636301667),\n",
       " (1524, 0.09243571975023783),\n",
       " (1525, 0.035512485566261465),\n",
       " (1526, 0.011611956395607774),\n",
       " (1527, 0.0032523832966202416),\n",
       " (1528, 0.0008046953251989233),\n",
       " (1529, 0.00017802833709978724)]"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "myComposition_brominated = myPeptideSequence.getFormula(0, 0).toString() + b'Br'\n",
    "\n",
    "formulaObject = EmpiricalFormula(myComposition_brominated)\n",
    "\n",
    "brominatedDistribution = formulaObject.getIsotopeDistribution(10)\n",
    "\n",
    "brominatedDistribution.getContainer()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "I hope that I've given you a good idea of what basic features OpenMS is capable of and that you're starting to have some ideas about how to harness this power for your own research. OpenMS has hundreds of classes and functions and I'll try to keep this SOP updated with the most useful ones. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.6.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}