usALEX - Corrections - Gamma factor fit.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fit Gamma factor\n",
    "\n",
    "> *This notebook estimates the gamma factor from a set of 5 μs-ALEX smFRET measurements.*\n",
    "\n",
    "## What this notebook does?\n",
    "\n",
    "According to [Lee 2005](http://dx.doi.org/10.1529/biophysj.104.054114) ([PDF](http://www.chem.ucla.edu/~michalet/papers/BJ2005.pdf), [SI PDF](http://www.chem.ucla.edu/~michalet/papers/BJ2005SI.pdf)), we estimate the $\\gamma$-factor \n",
    "from Proximity Ratio (PR) and S values (with background, leakage and direct excitation correction) \n",
    "for a set of 5 μs-ALEX measurements.\n",
    "\n",
    "The PR and S values are computed by the notebook\n",
    "\n",
    "- [usALEX-5samples-PR-leakage-dir-ex-all-ph](usALEX-5samples-PR-leakage-dir-ex-all-ph.ipynb)\n",
    "\n",
    "which is executed by [8-spots paper analysis](8-spots paper analysis.ipynb).\n",
    "\n",
    "From [Lee 2005](http://dx.doi.org/10.1529/biophysj.104.054114) (equation 20), the following linear relation holds:\n",
    "\n",
    "$$\\frac{1}{S} = \\Omega + \\Sigma \\cdot E_{PR}$$\n",
    "\n",
    "Once $\\Omega$ and $\\Sigma$ are fitted, we can compute the $\\gamma$-factor as ([equation 22](http://www.sciencedirect.com/science/article/pii/S0006349505733464#eq22)):\n",
    "\n",
    "$$\\gamma = (\\Omega-1)/(\\Omega + \\Sigma-1)$$\n",
    "\n",
    "$$\\beta = \\Omega + \\Sigma - 1$$\n",
    "\n",
    "The definition of $\\beta$ based on physical parameters is:\n",
    "\n",
    "$$ \\beta = \\frac{I_{A_{ex}}\\sigma_{A_{ex}}^A}{I_{D_{ex}}\\sigma_{D_{ex}}^D}$$\n",
    "\n",
    "Note that, calling $S_\\gamma$ the corrected S, the following relation holds:\n",
    "\n",
    "$$ S_\\gamma = (1 + \\beta)^{-1}$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from __future__ import division\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import lmfit\n",
    "from scipy.stats import linregress"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Computation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook read data from the file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data_file = 'results/usALEX-5samples-PR-leakage-dir-ex-all-ph.csv'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data = pd.read_csv(data_file).set_index('sample')\n",
    "data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data[['E_gauss_w', 'E_kde_w', 'S_gauss']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "E_ref, S_ref = data.E_gauss_w, data.S_gauss"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "res = linregress(E_ref, 1/S_ref)\n",
    "slope, intercept, r_val, p_val, stderr = res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For more info see [`scipy.stats.linearregress`](http://docs.scipy.org/doc/scipy-0.13.0/reference/generated/scipy.stats.linregress.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Sigma = slope \n",
    "Sigma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Omega = intercept\n",
    "Omega"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Pearson correlation coefficient](http://en.wikipedia.org/wiki/Pearson_product-moment_correlation_coefficient):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "r_val"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[Coefficient of determination](http://en.wikipedia.org/wiki/Coefficient_of_determination) $R^2$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "r_val**2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "P-value (to test the null hypothesis that the slope is zero):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "p_val"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Gamma computed from the previous fitted values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "gamma = (Omega - 1)/(Omega + Sigma - 1)\n",
    "'%.6f' % gamma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "with open('results/usALEX - gamma factor - all-ph.csv', 'w') as f:\n",
    "    f.write('%.6f' % gamma)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "beta = Omega + Sigma - 1\n",
    "'%.6f' % beta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "with open('results/usALEX - beta factor - all-ph.csv', 'w') as f:\n",
    "    f.write('%.6f' % beta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fit plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "%matplotlib inline\n",
    "%config InlineBackend.figure_format='retina'  # for hi-dpi displays\n",
    "\n",
    "sns.set_style('whitegrid')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "x = np.arange(0, 1, 0.01)\n",
    "plt.plot(E_ref, 1/S_ref, 's', label='dsDNA samples')\n",
    "plt.plot(x, intercept + slope*x, 'k', label='fit (slope = %.2f)' % slope)\n",
    "plt.legend(loc=4)\n",
    "plt.ylim(1, 2)\n",
    "plt.xlabel('PR')\n",
    "plt.ylabel('1/SR');"
   ]
  }
 ],
 "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}