Multi-spot 5-Samples analysis - PR-raw.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Analysis of 5 smFRET samples\n",
    "\n",
    "In this notebook:\n",
    "\n",
    "* Per-sample FRET histograms\n",
    "* Collapsed (merged channels) FRET histograms\n",
    "* FRET vs distance plot (multi-spot vs usALEX)\n",
    "* D-only fraction plots\n",
    "* Burst size vs measurement (corrected and uncorrected)\n",
    "* Background (per CH) vs measuremets\n",
    "* Burst-per-second (per CH) vs measuremet\n",
    "* Burst accumulation vs time (per CH and all-CH)\n",
    "\n",
    "> **USAGE TIP**: to comment-out a code cell, select all (CTRL+a) and hit CTRL+/.\n",
    "\n",
    "## FRET fitting remarks\n",
    "\n",
    "In this notebook we fit the Proximity Ratio histogram using different models."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load FRETBursts software"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from fretbursts import *\n",
    "sns = init_notebook()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import pandas as pd\n",
    "from IPython.display import display\n",
    "from IPython.display import display, Math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import ipywidgets as widgets\n",
    "from ipywidgets import interact, interactive, fixed"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import lmfit\n",
    "print('lmfit version:', lmfit.__version__)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 8-spot paper plot style"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "PLOT_DIR = './figure/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "import matplotlib as mpl\n",
    "from cycler import cycler\n",
    "\n",
    "bmap = sns.color_palette(\"Set1\", 9)\n",
    "colors = np.array(bmap)[(1,0,2,3,4,8,6,7), :]\n",
    "mpl.rcParams['axes.prop_cycle'] = cycler('color', colors)\n",
    "colors_labels = ['blue', 'red', 'green', 'violet', 'orange', 'gray', 'brown', 'pink', ]\n",
    "for c, cl in zip(colors, colors_labels):\n",
    "    locals()[cl] = tuple(c) # assign variables with color names\n",
    "sns.palplot(colors)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Data folder:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data_dir = './data/multispot/'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check that the folder exists:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data_dir = os.path.abspath(data_dir) + '/'\n",
    "assert os.path.exists(data_dir), \"Path '%s' does not exist.\" % data_dir"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "List of data files in `data_dir`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "from glob import glob\n",
    "file_list = sorted(glob(data_dir + '*.hdf5'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "labels = ['7d', '12d', '17d', '22d', '27d', 'DO']\n",
    "files_dict = {lab: fname for lab, fname in zip(sorted(labels), file_list)}\n",
    "files_dict"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parameters from file"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the **leakage coefficient** from disk (computed in [Multi-spot 5-Samples analyis - Leakage coefficient fit](Multi-spot 5-Samples analyis - Leakage coefficient fit.ipynb)):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "leakage_coeff_fname = 'results/Multi-spot - leakage coefficient KDE wmean DexDem.csv'\n",
    "leakage = np.loadtxt(leakage_coeff_fname, ndmin=1)\n",
    "\n",
    "print('Leakage coefficient:', leakage)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load the **direct excitation coefficient** ($d_{dirT}$) from disk (computed in [usALEX - Corrections - Direct excitation physical parameter](usALEX - Corrections - Direct excitation physical parameter.ipynb)):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dir_ex_coeff_fname = 'results/usALEX - direct excitation coefficient dir_ex_t beta.csv'\n",
    "dir_ex_t = np.loadtxt(dir_ex_coeff_fname, ndmin=1)\n",
    "\n",
    "print('Direct excitation coefficient (dir_ex_t):', dir_ex_t)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Analysis parameters:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "gamma_sel = 0.44     # Used to compute burst size during burst selection\n",
    "donor_ref = False    # False -> gamma correction is: g*nd + na\n",
    "                     # True  -> gamma correction is: nd + na/g\n",
    "\n",
    "hist_weights = 'size'\n",
    "\n",
    "## Background fit parameters\n",
    "bg_kwargs_auto = dict(fun=bg.exp_fit,\n",
    "                 time_s = 30,\n",
    "                 tail_min_us = 'auto',\n",
    "                 F_bg=1.7,\n",
    "                 )\n",
    "\n",
    "## Burst search\n",
    "F=6\n",
    "dither = False\n",
    "size_th = 30    # Burst size threshold (selection on corrected burst sizes)\n",
    "\n",
    "## FRET fit parameters\n",
    "bandwidth = 0.03        # KDE bandwidth\n",
    "E_range = {'7d':  (0.7, 1.0), '12d': (0.4, 0.8), '17d': (0.2, 0.4), \n",
    "           '22d': (0.0, 0.1), '27d': (0.0, 0.1), 'DO': (0.0, 0.1)}\n",
    "E_axis_kde = np.arange(-0.2, 1.2, 0.0002)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Processing and plot options:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# Data load options\n",
    "reload_data = 1\n",
    "burst_search = 1\n",
    "delete_ph_times = 0"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "# Plot output options\n",
    "save_figure = 0\n",
    "nosuptitle = False\n",
    "\n",
    "plt.rc('savefig', dpi=75)  # Changes the figure size in the notebook\n",
    "savefig_kwargs =  dict(dpi=200, bbox_inches='tight')  # default save-figure options\n",
    "\n",
    "fret_plot_kw = dict(bins = np.r_[-0.2:1.201:bandwidth],\n",
    "                    fit_color='k', fit_alpha=0.9, fit_lw=2, #fit_fillcolor='#888888',\n",
    "                    sharey = False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Utility functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def print_fit_report(E_pr, gamma=1, leakage=0, dir_ex_t=0, math=True):\n",
    "    \"\"\"Print fit and standard deviation for both corrected and uncorrected E\n",
    "    Returns d.E_fit.\n",
    "    \"\"\"\n",
    "    E_corr = fretmath.correct_E_gamma_leak_dir(E_pr, gamma=gamma, leakage=leakage, dir_ex_t=dir_ex_t)\n",
    "    \n",
    "    E_pr_mean = E_pr.mean()*100\n",
    "    E_pr_delta = (E_pr.max() - E_pr.min())*100\n",
    "    \n",
    "    E_corr_mean = E_corr.mean()*100\n",
    "    E_corr_delta = (E_corr.max() - E_corr.min())*100\n",
    "    if math:\n",
    "        display(Math(r'\\text{Pre}\\;\\gamma\\quad\\langle{E}_{fit}\\rangle = %.1f\\%% \\qquad'\n",
    "                     '\\Delta E_{fit} = %.2f \\%%' % \\\n",
    "                     (E_pr_mean, E_pr_delta)))\n",
    "        display(Math(r'\\text{Post}\\;\\gamma\\quad\\langle{E}_{fit}\\rangle = %.1f\\%% \\qquad'\n",
    "                     '\\Delta E_{fit} = %.2f \\%%' % \\\n",
    "                     (E_corr_mean, E_corr_delta)))\n",
    "    else:\n",
    "        print('Pre-gamma  E (delta, mean):  %.2f  %.2f' % (E_pr_mean, E_pr_delta))\n",
    "        print('Post-gamma E (delta, mean):  %.2f  %.2f' % (E_corr_mean, E_corr_delta))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5-samples processing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "df = pd.DataFrame(index=['7d', '12d', '17d', '22d', '27d'], columns=range(8), dtype=float)\n",
    "df.index.name = 'Sample'\n",
    "df.columns.name = 'Channel'\n",
    "\n",
    "E_pr_fret = df.copy()\n",
    "E_pr_fret_sig = df.copy()\n",
    "nbursts = df.copy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7bp sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data_id = '7d'\n",
    "if reload_data:\n",
    "    d7 = loader.photon_hdf5(files_dict[data_id])\n",
    "    d7.calc_bg(**bg_kwargs_auto)\n",
    "if burst_search:\n",
    "    d7.burst_search(m=10, F=F, dither=dither)\n",
    "    if delete_ph_times: d7.delete('ph_times_m')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "gamma_sel, donor_ref"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dfs7 = Sel(d7, select_bursts.size, th1=30, gamma=gamma_sel, donor_ref=donor_ref)\n",
    "dx = dfs7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "nbursts.loc[data_id] = dx.num_bursts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fitter = bext.bursts_fitter(dx)\n",
    "fitter.histogram(bins=np.r_[-0.2 : 1.2 : bandwidth])\n",
    "fitter.model = mfit.factory_two_gaussians(add_bridge=True, p2_center=0.8)\n",
    "fitter.fit_histogram()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "E_pr_fret.loc[data_id] = fitter.params['p2_center']\n",
    "display(E_pr_fret.loc[[data_id]])\n",
    "print_fit_report(E_pr_fret.loc[data_id], gamma=0.42, leakage=leakage, dir_ex_t=dir_ex_t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dplot(dx, hist_fret, show_model=True, \n",
    "      show_fit_stats=True, fit_from='p2_center', show_fit_value=True);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(4, 2, figsize=(14, 12), sharex=True, sharey=True)\n",
    "fig.subplots_adjust(left=0.08, right=0.96, top=0.93, bottom=0.07,\n",
    "                    wspace=0.06, hspace=0.15)\n",
    "\n",
    "for ich, ax in enumerate(axes.ravel()):\n",
    "    mfit.plot_mfit(fitter, ich=ich, ax=ax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compare the effect of burst-size weights"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# bext.bursts_fitter(dx, weights=None)\n",
    "# dx.E_fitter.fit_histogram(mfit.factory_two_gaussians())\n",
    "# dplot(dx, hist_fret, weights=None, show_model=True, show_fit_stats=True, fit_from='p2_center');\n",
    "# ylim(0, 6)\n",
    "# fig_no_w = gcf()\n",
    "# plt.close(fig_no_w)\n",
    "\n",
    "# bext.bursts_fitter(dx, weights='size', gamma=0.43)\n",
    "# dx.E_fitter.fit_histogram(mfit.factory_two_gaussians())\n",
    "# dplot(dx, hist_fret, weights='size', gamma=0.43,\n",
    "#       show_model=True, show_fit_stats=True, fit_from='p2_center');\n",
    "# fig_w = gcf()\n",
    "# plt.close(fig_w)\n",
    "\n",
    "# def _plot(weights=False):\n",
    "#     if weights:\n",
    "#         display(fig_w)\n",
    "#     else:\n",
    "#         display(fig_no_w)\n",
    "\n",
    "# interact(_plot, weights=False);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 12bp sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data_id = '12d'\n",
    "if reload_data: \n",
    "    d12 = loader.photon_hdf5(files_dict[data_id])\n",
    "    d12.calc_bg_cache(**bg_kwargs_auto)\n",
    "if burst_search:\n",
    "    d12.burst_search(m=10, F=F, dither=dither)\n",
    "    if delete_ph_times: d12.delete('ph_times_m')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dfs12 = Sel(d12, select_bursts.size, th1=30, gamma=gamma_sel, donor_ref=donor_ref)\n",
    "dx = dfs12"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "nbursts.loc[data_id] = dx.num_bursts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fitter = bext.bursts_fitter(dx)\n",
    "fitter.histogram(bins=np.r_[-0.2 : 1.2 : bandwidth])\n",
    "fitter.model = mfit.factory_two_gaussians(add_bridge=True, p2_center=0.65)\n",
    "fitter.fit_histogram()\n",
    "E_pr_fret.loc[data_id] = fitter.params['p2_center']\n",
    "display(E_pr_fret.loc[[data_id]])\n",
    "print_fit_report(E_pr_fret.loc[data_id], gamma=0.42, leakage=leakage, dir_ex_t=dir_ex_t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dplot(dx, hist_fret, show_model=True, \n",
    "      show_fit_stats=True, fit_from='p2_center', show_fit_value=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 17bp sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data_id = '17d'\n",
    "if reload_data:\n",
    "    d17 = loader.photon_hdf5(files_dict[data_id])\n",
    "    d17.calc_bg_cache(**bg_kwargs_auto)\n",
    "if burst_search:\n",
    "    d17.burst_search(m=10, F=F, dither=dither)\n",
    "    if delete_ph_times: d17.delete('ph_times_m')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dfs17 = Sel(d17, select_bursts.size, th1=30, gamma=gamma_sel, donor_ref=donor_ref)\n",
    "dx = dfs17"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "nbursts.loc[data_id] = dx.num_bursts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fitter = bext.bursts_fitter(dx)\n",
    "fitter.histogram(bins=np.r_[-0.2 : 1.2 : bandwidth])\n",
    "fitter.model = mfit.factory_two_gaussians(add_bridge=False, p2_center=0.4)\n",
    "fitter.fit_histogram()\n",
    "E_pr_fret.loc[data_id] = fitter.params['p2_center']\n",
    "display(E_pr_fret.loc[[data_id]])\n",
    "print_fit_report(E_pr_fret.loc[data_id], gamma=0.42, leakage=leakage, dir_ex_t=dir_ex_t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dplot(dx, hist_fret, show_model=True, \n",
    "      show_fit_stats=True, fit_from='p2_center', show_fit_value=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 22bp sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data_id = '22d'\n",
    "if reload_data:\n",
    "    d22 = loader.photon_hdf5(files_dict[data_id])\n",
    "    d22.calc_bg_cache(**bg_kwargs_auto)\n",
    "if burst_search:\n",
    "    d22.burst_search(m=10, F=F, dither=dither)\n",
    "    if delete_ph_times: d22.delete('ph_times_m')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dfs22 = Sel(d22, select_bursts.size, th1=30, gamma=gamma_sel, donor_ref=donor_ref)\n",
    "dx = dfs22"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "nbursts.loc[data_id] = dx.num_bursts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fitter = bext.bursts_fitter(dx)\n",
    "fitter.histogram(bins=np.r_[-0.2 : 1.2 : bandwidth])\n",
    "fitter.model = mfit.factory_gaussian()\n",
    "fitter.fit_histogram()\n",
    "E_pr_fret.loc[data_id] = fitter.params['center']\n",
    "display(E_pr_fret.loc[[data_id]])\n",
    "print_fit_report(E_pr_fret.loc[data_id], gamma=0.42, leakage=leakage, dir_ex_t=dir_ex_t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dplot(dx, hist_fret, show_model=True, \n",
    "      show_fit_stats=True, fit_from='center', show_fit_value=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 27bp sample"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data_id = '27d'\n",
    "if reload_data:\n",
    "    d27 = loader.photon_hdf5(files_dict[data_id])\n",
    "    d27.calc_bg_cache(**bg_kwargs_auto)\n",
    "if burst_search:\n",
    "    d27.burst_search(m=10, F=F, dither=dither)#, ph_sel=Ph_sel(Dex='Dem'))\n",
    "    if delete_ph_times: d27.delete('ph_times_m')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dfs27 = Sel(d27, select_bursts.size, th1=30, gamma=gamma_sel, donor_ref=donor_ref)\n",
    "dx = dfs27"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "nbursts.loc[data_id] = dx.num_bursts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fitter = bext.bursts_fitter(dx)\n",
    "fitter.histogram(bins=np.r_[-0.2 : 1.2 : bandwidth])\n",
    "fitter.model = mfit.factory_asym_gaussian()\n",
    "fitter.fit_histogram()\n",
    "E_pr_fret.loc[data_id] = fitter.params['center']\n",
    "display(E_pr_fret.loc[[data_id]])\n",
    "print_fit_report(E_pr_fret.loc[data_id], gamma=0.42, leakage=leakage, dir_ex_t=dir_ex_t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dplot(dx, hist_fret, show_model=True, \n",
    "      show_fit_stats=True, fit_from='center', show_fit_value=True);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fitter = bext.bursts_fitter(dx)\n",
    "fitter.histogram(bins=np.r_[-0.2 : 1.2 : bandwidth])\n",
    "fitter.model = mfit.factory_gaussian()\n",
    "fitter.fit_histogram()\n",
    "E_pr_fret.loc[data_id] = fitter.params['center']\n",
    "display(E_pr_fret.loc[[data_id]])\n",
    "print_fit_report(E_pr_fret.loc[data_id], gamma=0.42, leakage=leakage, dir_ex_t=dir_ex_t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dplot(dx, hist_fret, show_model=True, \n",
    "      show_fit_stats=True, fit_from='center', show_fit_value=True);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5-samples: global analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Proximity ratios fitted from multispot data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "E_pr_fret = E_pr_fret.round(6)\n",
    "E_pr_fret"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "nbursts = nbursts.astype(int)\n",
    "nbursts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "E_pr_fret.to_csv('results/Multi-spot - dsDNA - PR - all_samples all_ch.csv')\n",
    "nbursts.to_csv('results/Multi-spot - dsDNA - nbursts - all_samples all_ch.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "norm = (E_pr_fret.T - E_pr_fret.mean(1))#/E_pr_fret.mean(1)\n",
    "norm_rel = (E_pr_fret.T - E_pr_fret.mean(1))/E_pr_fret.mean(1)\n",
    "norm.plot()\n",
    "norm_rel.plot()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "mch_plot_bg(dx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **NOTE:** The 27d and DO samples have a trend that correlates with the A-ch background. For these samples could be beneficial to use a D-only burst search.\n",
    "\n",
    "> **NOTE 2**: Like observed during the leakage fit, even a D-only burst search results in +2% offset (DO) in CH1, this cannot be correlation with the A-background and must be slightly different gamma in the spot."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Corrected $E$ from μs-ALEX data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "data_file = 'results/usALEX-5samples-E-corrected-all-ph.csv'\n",
    "data_alex = pd.read_csv(data_file).set_index('sample')#[['E_pr_fret_kde']]\n",
    "data_alex.round(6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "E_alex = data_alex.E_gauss_w\n",
    "E_alex"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Merging the channels"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Merge the data of the different channels:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "dfs7c = dfs7.collapse()\n",
    "dfs12c = dfs12.collapse()\n",
    "dfs17c = dfs17.collapse()\n",
    "dfs22c = dfs22.collapse()\n",
    "dfs27c = dfs27.collapse()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define samples lists"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Define list of results and labels:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "d_samples = [dfs7, dfs12, dfs17, dfs22, dfs27]#, dfso]\n",
    "d_samples_c = [dfs7c, dfs12c, dfs17c, dfs22c, dfs27c ]\n",
    "d_labels = ['7d', '12d', '17d', '22d', '27d']#, 'DO']\n",
    "CH = np.arange(8)\n",
    "CH_labels = ['CH%d' % i for i in CH]\n",
    "dist_s_bp = [7, 12, 17, 22, 27]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Print a summary of current processed data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def print_params(d_samples, d_labels, status=False):\n",
    "    print('Sample              Model             Ph_sel')\n",
    "    for dx, name in zip(d_samples, d_labels):\n",
    "        print(\"%3s %25s %35s\" % (name, dx.E_fitter.model.name, dx.ph_sel))\n",
    "    if status:\n",
    "        print()\n",
    "        for dx, name in zip(d_samples, d_labels):\n",
    "            print(dx.status())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "print_params(d_samples, d_labels, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot FRET vs distance"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {
    "collapsed": false
   },
   "source": [
    "fontsize = 12\n",
    "text_pos = {'7d': (10.5, 0.92), '12d': (15.5, 0.74), '17d': (16, 0.43), '22d': (21, 0.18), '27d': (29, 0.15)}\n",
    "text_kwargs = dict(ha='right', va='center', \n",
    "                   bbox=dict(facecolor='white', edgecolor='white'),\n",
    "                   zorder=2, fontsize=fontsize-1,\n",
    "                   )\n",
    "dfun = lambda E_fit: 100*(E_fit.max()-E_fit.min())\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(dist_s_bp, E_fret_mch, '+', lw=2, mew=1.2, ms=8, zorder=4)\n",
    "ax.plot(dist_s_bp, E_alex, '-', lw=3, mew=0, alpha=0.5, color='r', zorder=3)\n",
    "\n",
    "plt.title('Multi-spot smFRET dsDNA, Gamma = %.2f' % multispot_gamma)\n",
    "plt.xlabel('Distance in base-pairs', fontsize=fontsize); \n",
    "plt.ylabel('E', fontsize=fontsize)\n",
    "plt.ylim(0, 1); plt.xlim(0, 30)\n",
    "plt.grid(True)\n",
    "plt.legend(['CH1','CH2','CH3','CH4','CH5','CH6','CH7','CH8', u'μsALEX'], \n",
    "       fancybox=True, prop={'size':fontsize-1},\n",
    "       #loc='upper right', bbox_to_anchor=(1.05, 1))\n",
    "       loc='best')\n",
    "\n",
    "if deviance_plot:\n",
    "    for sample in E_fret_mch.index:\n",
    "        delta = dfun(E_fret_mch.loc[sample])\n",
    "        text_kwargs.update(\n",
    "                          s=\"%.1f%%\" % delta, \n",
    "                          #s=\"$\\Delta = %.1f%%$\" % delta, \n",
    "                          x=text_pos[sample][0],\n",
    "                          y=text_pos[sample][1],\n",
    "                          )\n",
    "        ax.text(**text_kwargs)\n",
    "\n",
    "#ax.set_axisbelow(True)\n",
    "#if save_figure: plt.savefig(\"FRET vs distance - %s SPW F%d.png\" % (fret_gen_fit_func.__name__, F))\n",
    "#plt.savefig(PLOT_DIR+\"FRET vs distance - Gpost_th30.png\", dpi=200, bbox_inches='tight')"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [default]",
   "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}