documentation_builder/sampling.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Flux sampling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basic usage"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The easiest way to get started with flux sampling is using the `sample` function in the `flux_analysis` submodule. `sample` takes at least two arguments: a cobra model and the number of samples you want to generate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>ACALD</th>\n",
       "      <th>ACALDt</th>\n",
       "      <th>ACKr</th>\n",
       "      <th>ACONTa</th>\n",
       "      <th>ACONTb</th>\n",
       "      <th>ACt2r</th>\n",
       "      <th>ADK1</th>\n",
       "      <th>AKGDH</th>\n",
       "      <th>AKGt2r</th>\n",
       "      <th>ALCD2x</th>\n",
       "      <th>...</th>\n",
       "      <th>RPI</th>\n",
       "      <th>SUCCt2_2</th>\n",
       "      <th>SUCCt3</th>\n",
       "      <th>SUCDi</th>\n",
       "      <th>SUCOAS</th>\n",
       "      <th>TALA</th>\n",
       "      <th>THD2</th>\n",
       "      <th>TKT1</th>\n",
       "      <th>TKT2</th>\n",
       "      <th>TPI</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>-3.706944</td>\n",
       "      <td>-0.163964</td>\n",
       "      <td>-0.295823</td>\n",
       "      <td>8.975852</td>\n",
       "      <td>8.975852</td>\n",
       "      <td>-0.295823</td>\n",
       "      <td>4.847986</td>\n",
       "      <td>6.406533</td>\n",
       "      <td>-0.081797</td>\n",
       "      <td>-3.542980</td>\n",
       "      <td>...</td>\n",
       "      <td>-1.649393</td>\n",
       "      <td>20.917568</td>\n",
       "      <td>20.977290</td>\n",
       "      <td>744.206008</td>\n",
       "      <td>-6.406533</td>\n",
       "      <td>1.639515</td>\n",
       "      <td>1.670533</td>\n",
       "      <td>1.639515</td>\n",
       "      <td>1.635542</td>\n",
       "      <td>6.256787</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>-1.340710</td>\n",
       "      <td>-0.175665</td>\n",
       "      <td>-0.429169</td>\n",
       "      <td>11.047827</td>\n",
       "      <td>11.047827</td>\n",
       "      <td>-0.429169</td>\n",
       "      <td>2.901598</td>\n",
       "      <td>7.992916</td>\n",
       "      <td>-0.230564</td>\n",
       "      <td>-1.165045</td>\n",
       "      <td>...</td>\n",
       "      <td>-0.066975</td>\n",
       "      <td>24.735567</td>\n",
       "      <td>24.850041</td>\n",
       "      <td>710.481004</td>\n",
       "      <td>-7.992916</td>\n",
       "      <td>0.056442</td>\n",
       "      <td>9.680476</td>\n",
       "      <td>0.056442</td>\n",
       "      <td>0.052207</td>\n",
       "      <td>7.184752</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>-1.964087</td>\n",
       "      <td>-0.160334</td>\n",
       "      <td>-0.618029</td>\n",
       "      <td>9.811474</td>\n",
       "      <td>9.811474</td>\n",
       "      <td>-0.618029</td>\n",
       "      <td>17.513791</td>\n",
       "      <td>8.635576</td>\n",
       "      <td>-0.284992</td>\n",
       "      <td>-1.803753</td>\n",
       "      <td>...</td>\n",
       "      <td>-4.075515</td>\n",
       "      <td>23.425719</td>\n",
       "      <td>23.470968</td>\n",
       "      <td>696.114154</td>\n",
       "      <td>-8.635576</td>\n",
       "      <td>4.063291</td>\n",
       "      <td>52.316496</td>\n",
       "      <td>4.063291</td>\n",
       "      <td>4.058376</td>\n",
       "      <td>5.122237</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>-0.838442</td>\n",
       "      <td>-0.123865</td>\n",
       "      <td>-0.376067</td>\n",
       "      <td>11.869552</td>\n",
       "      <td>11.869552</td>\n",
       "      <td>-0.376067</td>\n",
       "      <td>7.769872</td>\n",
       "      <td>9.765178</td>\n",
       "      <td>-0.325219</td>\n",
       "      <td>-0.714577</td>\n",
       "      <td>...</td>\n",
       "      <td>-0.838094</td>\n",
       "      <td>23.446704</td>\n",
       "      <td>23.913036</td>\n",
       "      <td>595.787313</td>\n",
       "      <td>-9.765178</td>\n",
       "      <td>0.822987</td>\n",
       "      <td>36.019720</td>\n",
       "      <td>0.822987</td>\n",
       "      <td>0.816912</td>\n",
       "      <td>8.364314</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>-0.232088</td>\n",
       "      <td>-0.034346</td>\n",
       "      <td>-1.067684</td>\n",
       "      <td>7.972039</td>\n",
       "      <td>7.972039</td>\n",
       "      <td>-1.067684</td>\n",
       "      <td>5.114975</td>\n",
       "      <td>5.438125</td>\n",
       "      <td>-0.787864</td>\n",
       "      <td>-0.197742</td>\n",
       "      <td>...</td>\n",
       "      <td>-3.109205</td>\n",
       "      <td>8.902309</td>\n",
       "      <td>9.888083</td>\n",
       "      <td>584.552692</td>\n",
       "      <td>-5.438125</td>\n",
       "      <td>3.088152</td>\n",
       "      <td>12.621811</td>\n",
       "      <td>3.088152</td>\n",
       "      <td>3.079686</td>\n",
       "      <td>6.185089</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 95 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "      ACALD    ACALDt      ACKr     ACONTa     ACONTb     ACt2r       ADK1  \\\n",
       "0 -3.706944 -0.163964 -0.295823   8.975852   8.975852 -0.295823   4.847986   \n",
       "1 -1.340710 -0.175665 -0.429169  11.047827  11.047827 -0.429169   2.901598   \n",
       "2 -1.964087 -0.160334 -0.618029   9.811474   9.811474 -0.618029  17.513791   \n",
       "3 -0.838442 -0.123865 -0.376067  11.869552  11.869552 -0.376067   7.769872   \n",
       "4 -0.232088 -0.034346 -1.067684   7.972039   7.972039 -1.067684   5.114975   \n",
       "\n",
       "      AKGDH    AKGt2r    ALCD2x    ...          RPI   SUCCt2_2     SUCCt3  \\\n",
       "0  6.406533 -0.081797 -3.542980    ...    -1.649393  20.917568  20.977290   \n",
       "1  7.992916 -0.230564 -1.165045    ...    -0.066975  24.735567  24.850041   \n",
       "2  8.635576 -0.284992 -1.803753    ...    -4.075515  23.425719  23.470968   \n",
       "3  9.765178 -0.325219 -0.714577    ...    -0.838094  23.446704  23.913036   \n",
       "4  5.438125 -0.787864 -0.197742    ...    -3.109205   8.902309   9.888083   \n",
       "\n",
       "        SUCDi    SUCOAS      TALA       THD2      TKT1      TKT2       TPI  \n",
       "0  744.206008 -6.406533  1.639515   1.670533  1.639515  1.635542  6.256787  \n",
       "1  710.481004 -7.992916  0.056442   9.680476  0.056442  0.052207  7.184752  \n",
       "2  696.114154 -8.635576  4.063291  52.316496  4.063291  4.058376  5.122237  \n",
       "3  595.787313 -9.765178  0.822987  36.019720  0.822987  0.816912  8.364314  \n",
       "4  584.552692 -5.438125  3.088152  12.621811  3.088152  3.079686  6.185089  \n",
       "\n",
       "[5 rows x 95 columns]"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from cobra.test import create_test_model\n",
    "from cobra.flux_analysis import sample\n",
    "\n",
    "model = create_test_model(\"textbook\")\n",
    "s = sample(model, 100)\n",
    "s.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By default sample uses the `optgp` method based on the [method presented here](http://dx.doi.org/10.1371/journal.pone.0086587) as it is suited for larger models and can run in parallel. By default the sampler uses a single process. This can be changed by using the `processes` argument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "One process:\n",
      "CPU times: user 5.31 s, sys: 433 ms, total: 5.74 s\n",
      "Wall time: 5.27 s\n",
      "Two processes:\n",
      "CPU times: user 217 ms, sys: 488 ms, total: 705 ms\n",
      "Wall time: 2.8 s\n"
     ]
    }
   ],
   "source": [
    "print(\"One process:\")\n",
    "%time s = sample(model, 1000)\n",
    "print(\"Two processes:\")\n",
    "%time s = sample(model, 1000, processes=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Alternatively you can also user Artificial Centering Hit-and-Run for sampling by setting the method to `achr`.  `achr` does not support parallel execution but has good convergence and is almost Markovian."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "s = sample(model, 100, method=\"achr\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In general setting up the sampler is expensive since initial search directions are generated by solving many linear programming problems. Thus, we recommend to generate as many samples as possible in one go. However, this might require finer control over the sampling procedure as described in the following section."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Advanced usage"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sampler objects"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The sampling process can be controlled on a lower level by using the sampler classes directly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from cobra.flux_analysis.sampling import OptGPSampler, ACHRSampler"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Both sampler classes have standardized interfaces and take some additional argument. For instance the `thinning` factor. \"Thinning\" means only recording samples every n iterations. A higher thinning factors mean less correlated samples but also larger computation times. By default the samplers use a thinning factor of 100 which creates roughly uncorrelated samples. If you want less samples but better mixing feel free to increase this parameter. If you want to study convergence for your own model you might want to set it to 1 to obtain all iterates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "achr = ACHRSampler(model, thinning=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`OptGPSampler` has an additional `processes` argument specifying how many processes are used to create parallel sampling chains. This should be in the order of your CPU cores for maximum efficiency. As noted before class initialization can take up to a few minutes due to generation of initial search directions. Sampling on the other hand is quick."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "optgp = OptGPSampler(model, processes=4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sampling and validation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Both samplers have a sample function that generates samples from the initialized object and act like the `sample` function described above, only that this time it will only accept a single argument, the number of samples. For `OptGPSampler` the number of samples should be a multiple of the number of processes, otherwise it will be increased to the nearest multiple automatically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "s1 = achr.sample(100)\n",
    "\n",
    "s2 = optgp.sample(100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can call `sample` repeatedly and both samplers are optimized to generate large amount of samples without falling into \"numerical traps\". All sampler objects have a `validate` function in order to check if a set of points are feasible and give detailed information about feasibility violations in a form of a short code denoting feasibility. Here the short code is a combination of any of the following letters:\n",
    "\n",
    "- \"v\" - valid point\n",
    "- \"l\" - lower bound violation\n",
    "- \"u\" - upper bound violation\n",
    "- \"e\" - equality violation (meaning the point is not a steady state)\n",
    "\n",
    "For instance for a random flux distribution (should not be feasible):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['le'], \n",
       "      dtype='<U3')"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "bad = np.random.uniform(-1000, 1000, size=len(model.reactions))\n",
    "achr.validate(np.atleast_2d(bad))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And for our generated samples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(['v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v',\n",
       "       'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v',\n",
       "       'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v',\n",
       "       'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v',\n",
       "       'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v',\n",
       "       'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v',\n",
       "       'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v',\n",
       "       'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v', 'v'], \n",
       "      dtype='<U3')"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "achr.validate(s1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Batch sampling"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sampler objects are made for generating billions of samples, however using the `sample` function might quickly fill up your RAM when working with genome-scale models. Here, the `batch` method of the sampler objects might come in handy. `batch` takes two arguments, the number of samples in each batch and the number of batches. This will make sense with a small example. \n",
    "\n",
    "Let's assume we want to quantify what proportion of our samples will grow. For that we might want to generate 10 batches of 50 samples each and measure what percentage of the individual 100 samples show a growth rate larger than 0.1. Finally, we want to calculate the mean and standard deviation of those individual percentages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Usually 8.70% +- 2.72% grow...\n"
     ]
    }
   ],
   "source": [
    "counts = [np.mean(s.Biomass_Ecoli_core > 0.1) for s in optgp.batch(100, 10)]\n",
    "print(\"Usually {:.2f}% +- {:.2f}% grow...\".format(\n",
    "    np.mean(counts) * 100.0, np.std(counts) * 100.0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adding constraints"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Flux sampling will respect additional contraints defined in the model. For instance we can add a constraint enforcing growth in asimilar manner as the section before."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "co = model.problem.Constraint(model.reactions.Biomass_Ecoli_core.flux_expression, lb=0.1)\n",
    "model.add_cons_vars([co])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Note that this is only for demonstration purposes. usually you could set the lower bound of the reaction directly instead of creating a new constraint.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0    0.175547\n",
      "1    0.111499\n",
      "2    0.123073\n",
      "3    0.151874\n",
      "4    0.122541\n",
      "5    0.121878\n",
      "6    0.147333\n",
      "7    0.106499\n",
      "8    0.174448\n",
      "9    0.143273\n",
      "Name: Biomass_Ecoli_core, dtype: float64\n"
     ]
    }
   ],
   "source": [
    "s = sample(model, 10)\n",
    "print(s.Biomass_Ecoli_core)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see our new constraint was respected."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.5.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}