Merge pull request #41 from chungjjang80/master

Add a notebook for BVA

notebooks/Example - Burst Variance Analysis.ipynb

@@ -0,0 +1,318 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Example - Burst Variance Analysis\n",
+    "\n",
+    "*This notebook is part of smFRET burst analysis software [FRETBursts](http://tritemio.github.io/FRETBursts/).*"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "> In this notebook shows how to implement Burst Variance Analysis using FRETBursts.\n",
+    "\n",
+    "> For a complete tutorial on burst analysis see \n",
+    "> [FRETBursts - us-ALEX smFRET burst analysis](FRETBursts - us-ALEX smFRET burst analysis.ipynb)."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Loading the software"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "We start loading the **`FRETBursts`** software:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false,
+    "scrolled": false
+   },
+   "outputs": [],
+   "source": [
+    "from fretbursts import *\n",
+    "sns = init_notebook()"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The line after \"FRETBursts revision:\" tells us which FRETBursts revision is currently in use. This information can be used for long-term reproducibility, a crucial problem in science. Moreover, knowing the revision allows to easily track software regressions.\n",
+    "\n",
+    "Finally, a notice on how to cite FRETBursts is printed. Please follow this requirement by **always citing** FRETBursts in publications or presentations."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "url = 'http://files.figshare.com/2182601/0023uLRpitc_NTP_20dT_0.5GndCl.hdf5'\n",
+    "download_file(url, save_dir='./data')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "file_name = \"0023uLRpitc_NTP_20dT_0.5GndCl.hdf5\"\n",
+    "\n",
+    "# Here the folder is the subfolder \"data\" of current notebook folder\n",
+    "folder_name = './data/'\n",
+    "full_fname = folder_name + file_name\n",
+    "full_fname"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "import os\n",
+    "if os.path.isfile(full_fname):\n",
+    "    print (\"Perfect, I found the file!\")\n",
+    "else:\n",
+    "    print (\"Sorry, I can't find the file:\\n\", full_fname)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Load the selected file"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "d = loader.photon_hdf5(full_fname)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# us-ALEX parameters\n",
+    "\n",
+    "At this point, in `d`, we only have the timestamps (`ph_times_t`) and the detector numbers (`det_t`):"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "d.add(det_donor_accept=(0, 1), \n",
+    "      alex_period=4000, \n",
+    "      D_ON=(2100, 3900), \n",
+    "      A_ON=(150, 1900),\n",
+    "      offset=700)\n",
+    "bpl.plot_alternation_hist (d)\n",
+    "loader.usalex_apply_period(d)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The previous alternation histogram looks correct, so we can apply the current parameters:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "d.calc_bg(bg.exp_fit, time_s=50.1, tail_min_us='auto', F_bg=1.7)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Burst Variance Analysis"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "d.burst_search(m=10, computefret=False, ph_sel=Ph_sel(Dex='DAem'))\n",
+    "d.calc_fret(count_ph=True, corrections=False)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "ds = d.select_bursts(select_bursts.naa, th1=30, computefret=False)\n",
+    "ds1 = ds.select_bursts(select_bursts.size, th1=30, computefret=False)\n",
+    "ds_FRET = ds1.select_bursts(select_bursts.S, S1=0.25, S2=0.85, computefret=False)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "dx=ds_FRET\n",
+    "alex_jointplot(dx)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "bursts =ds_FRET.mburst[0]\n",
+    "ph_d = ds_FRET.get_ph_times(ph_sel=Ph_sel(Dex='DAem'))\n",
+    "AemDex_mask = ds_FRET.get_ph_mask(ph_sel=Ph_sel(Dex='Aem')) \n",
+    "Dex_mask = ds_FRET.get_ph_mask(ph_sel=Ph_sel(Dex='DAem'))   \n",
+    "AemDex_mask_d = AemDex_mask[Dex_mask]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "bursts_d = bursts.recompute_index_reduce(ph_d)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "n=7         # Determine the number of photons inside each sub-burst\n",
+    "\n",
+    "Std_small=[]\n",
+    "Avg_small=[]\n",
+    "for burst in bursts_d:\n",
+    "    E_small=[]\n",
+    "    startlist = range(burst.istart, burst.istop + 2 - n, n)\n",
+    "    stoplist = [i + n for i in startlist]\n",
+    "    for start, stop in zip(startlist, stoplist): \n",
+    "        A_D = AemDex_mask_d[start:stop].sum()\n",
+    "        assert stop-start == n\n",
+    "        E = A_D / n\n",
+    "        E_small.append(E)\n",
+    "    Std_small.append(np.std(E_small))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false,
+    "scrolled": false
+   },
+   "outputs": [],
+   "source": [
+    "sns.set_style(style='darkgrid')\n",
+    "plt.figure(figsize=(6,6))\n",
+    "x=np.arange(0,1.01,0.01)\n",
+    "y=np.sqrt((x*(1-x))/n)\n",
+    "plt.plot(x,y, lw=3, color='red')\n",
+    "im = sns.kdeplot(ds_FRET.E_, np.array(Std_small), shade=True, cmap='viridis', shade_lowest=False)\n",
+    "plt.xlim(0,1)\n",
+    "plt.ylim(0,0.4)\n",
+    "plt.xlabel('PR', fontsize=14)\n",
+    "plt.ylabel(r'$\\sigma_E$', fontsize=24)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": []
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "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
+}