Example - usALEX histogram

This notebook is part of smFRET burst analysis software FRETBursts.

In this notebook shows how to plot different styles of us-ALEX histograms and $E$ and $S$ marginal distributions. For a complete tutorial on burst analysis see FRETBursts - us-ALEX smFRET burst analysis.

Load software

FRETBursts

In [1]:
from fretbursts import *
 - Optimized (cython) burst search loaded.
 - Optimized (cython) photon counting loaded.
--------------------------------------------------------------
 You are running FRETBursts (version 0.5.7+0.g7267863.dirty).

 If you use this software please cite the following paper:

   FRETBursts: An Open Source Toolkit for Analysis of Freely-Diffusing Single-Molecule FRET
   Ingargiola et al. (2016). http://dx.doi.org/10.1371/journal.pone.0160716 

--------------------------------------------------------------
In [2]:
sns = init_notebook()

Interactive widgets

In [3]:
from ipywidgets import widgets, interact, interactive, fixed
from IPython.display import display, display_png, display_svg, clear_output
from IPython.core.pylabtools import print_figure

Download the sample file

In [4]:
url = 'http://files.figshare.com/2182601/0023uLRpitc_NTP_20dT_0.5GndCl.hdf5'
download_file(url, save_dir='./data')
URL:  http://files.figshare.com/2182601/0023uLRpitc_NTP_20dT_0.5GndCl.hdf5
File: 0023uLRpitc_NTP_20dT_0.5GndCl.hdf5
 
File already on disk: /Users/anto/src/FRETBursts/notebooks/data/0023uLRpitc_NTP_20dT_0.5GndCl.hdf5 
Delete it to re-download.

Select a data file

In [5]:
file_name = "0023uLRpitc_NTP_20dT_0.5GndCl.hdf5"

# Here the folder is the subfolder "data" of current notebook folder
folder_name = './data/'
full_fname = folder_name + file_name
full_fname
Out[5]:
'./data/0023uLRpitc_NTP_20dT_0.5GndCl.hdf5'

Let's check that the file exists:

In [6]:
import os
if os.path.isfile(full_fname):
    print("Perfect, I found the file!")
else:
    print("Sorry, I can't find the file:\n", full_fname)
Perfect, I found the file!

Load the selected file

In [7]:
d = loader.photon_hdf5(full_fname)
#d.add(det_donor_accept=(0, 1), alex_period=4000, 
#      offset=700, D_ON=(2180, 3900), A_ON=(200, 1800))
bpl.plot_alternation_hist(d)
loader.alex_apply_period(d)
d.calc_bg(bg.exp_fit, time_s=1000, tail_min_us=(800, 4000, 1500, 1000, 3000))
# Total photons (after ALEX selection):   2,259,522
#  D  photons in D+A excitation periods:    721,537
#  A  photons in D+A excitation periods:  1,537,985
# D+A photons in  D  excitation period:   1,434,842
# D+A photons in  A  excitation period:     824,680

 - Calculating BG rates ... [DONE]

Burst search

In [8]:
d.burst_search(L=10, m=10, F=6)
ds = Sel(d, select_bursts.size, add_naa=True, th1=30)
 - Performing burst search (verbose=False) ...[DONE]
 - Calculating burst periods ...[DONE]
 - Counting D and A ph and calculating FRET ... 
   - Applying background correction.
   - Applying leakage correction.
   - Applying direct excitation correction.
   [DONE Counting D/A]
In [9]:
#dplot(ds, hist_fret)
#dplot(ds, hist_S);
In [10]:
#ds.E_fitter.fit_histogram(mfit.factory_three_gaussians())

ALEX plots

Basics

We can make a simple E-S scatter plot with scatter_alex:

In [11]:
dplot(ds, scatter_alex, figsize=(4,4), mew=1, ms=4, mec='black', color='purple');
/Users/anto/miniconda3/lib/python3.5/site-packages/matplotlib/cbook.py:2644: UserWarning: Saw kwargs ['mew', 'markeredgewidth'] which are all aliases for 'markeredgewidth'.  Kept value from 'markeredgewidth'
  seen=seen, canon=canonical, used=seen[-1]))
/Users/anto/miniconda3/lib/python3.5/site-packages/matplotlib/cbook.py:2644: UserWarning: Saw kwargs ['mec', 'markeredgecolor'] which are all aliases for 'markeredgecolor'.  Kept value from 'markeredgecolor'
  seen=seen, canon=canonical, used=seen[-1]))
/Users/anto/miniconda3/lib/python3.5/site-packages/matplotlib/cbook.py:2644: UserWarning: Saw kwargs ['ms', 'markersize'] which are all aliases for 'markersize'.  Kept value from 'markersize'
  seen=seen, canon=canonical, used=seen[-1]))

We can also plot the ALEX histogram with a scatterplot overlay using hist2d_alex:

In [12]:
dplot(ds, hist2d_alex);

Colormaps

Here we plot a small selection of Matplotlib's colormaps:

In [13]:
cmaps = ['viridis', 'plasma', 'inferno', 'magma',
         'afmhot', 'Blues', 'BuGn', 'BuPu', 'GnBu', 'YlGnBu',
         'coolwarm', 'RdYlBu', 'RdYlGn', 'Spectral']
In [14]:
def plot_color_gradients(cmap_list):
    gradient = np.linspace(0, 1, 256)
    gradient = np.vstack((gradient, gradient))
    fig, axes = plt.subplots(nrows=len(cmap_list))
    fig.subplots_adjust(top=0.95, bottom=0.01, left=0.2, right=0.99)
    for ax, name in zip(axes, cmap_list):
        ax.imshow(gradient, aspect='auto', cmap=plt.get_cmap(name))
        pos = list(ax.get_position().bounds)
        x_text = pos[0] - 0.01
        y_text = pos[1] + pos[3]/2.
        fig.text(x_text, y_text, name, va='center', ha='right')
    for ax in axes:
        ax.set_axis_off()
In [15]:
plot_color_gradients(cmaps)

ALEX Joint Plot

In [16]:
sns.set_style('darkgrid')
In [17]:
# Fit E and S to a model and compute KDE
bext.bursts_fitter(ds, 'E', binwidth=0.03, bandwidth=0.03, model=mfit.factory_three_gaussians())
bext.bursts_fitter(ds, 'S', binwidth=0.03, bandwidth=0.03, model=mfit.factory_two_gaussians())
Out[17]:
<fretbursts.mfit.MultiFitter at 0x11bd56c50>
In [18]:
# Limiting the range of E and S works better for 2D KDE plots
dss = ds.select_bursts(select_bursts.ES, E1=-0.2, E2=1.2, S1=-0.2, S2=1.2)
In [19]:
# Remove D-only and A-only through burst selection
dssf = dss.select_bursts(select_bursts.size, th1=20)
dssf = dssf.select_bursts(select_bursts.naa, th1=20)

Just a test plot:

In [20]:
E, S = dss.E_, dss.S_
fig, ax = plt.subplots(figsize=(6, 6))
sns.kdeplot(E, S, shade=True, n_levels=30, cmap='alex_dark', shade_lowest=False)
plt.xlim(-0.19, 1.19)
plt.ylim(-0.19, 1.19)
Out[20]:
(-0.19, 1.19)

Joint plots

Here we show several jointplots. The argument joint_kws is reported with typical keywords used to customize the 2-D (inner) plot. FRETBursts defines 'alex_light' and 'alex_dark' colormaps in addition to the matplotlib and seaborn ones.

You can inver matplotlib colormaps appending _r to the name. For the full list see matplotlib colormaps.

In [21]:
alex_jointplot(ds, kind='scatter', cmap='Blues_r',
               joint_kws=dict(s=40, alpha=0.1, linewidths=0))
Out[21]:
<seaborn.axisgrid.JointGrid at 0x11ad01940>
In [22]:
alex_jointplot(dss, kind='kde', cmap='viridis',
               joint_kws=dict(shade=True, shade_lowest=False, n_levels=30))
Out[22]:
<seaborn.axisgrid.JointGrid at 0x11bcd5dd8>
In [23]:
alex_jointplot(ds, kind='hex', cmap='alex_dark',
               joint_kws=dict(edgecolor='grey', linewidth=0.2, gridsize=50))
Out[23]:
<seaborn.axisgrid.JointGrid at 0x11af96d30>
In [24]:
alex_jointplot(ds, kind='hex', cmap='inferno',
               joint_kws=dict(edgecolor='none', linewidth=0.001, gridsize=50))
Out[24]:
<seaborn.axisgrid.JointGrid at 0x11bc55f98>
In [25]:
alex_jointplot(dss, kind='kde', cmap='RdYlBu_r')
Out[25]:
<seaborn.axisgrid.JointGrid at 0x11afe2f28>
In [26]:
alex_jointplot(dssf, kind='kde', cmap='inferno',
               joint_kws=dict(shade=True, shade_lowest=False, n_levels=20))
Out[26]:
<seaborn.axisgrid.JointGrid at 0x11ad4e710>
In [27]:
alex_jointplot(ds, joint_kws = dict(edgecolor='gray'))

alex_jointplot(ds, joint_kws = {})
Out[27]:
<seaborn.axisgrid.JointGrid at 0x11afe20b8>
In [28]:
alex_jointplot(ds)
#plt.savefig('ALEX_jointplot.png', dpi=160, bbox_inches='tight')
Out[28]:
<seaborn.axisgrid.JointGrid at 0x11af96cf8>