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](FRETBursts - us-ALEX smFRET burst analysis.ipynb).
from fretbursts import *
- Optimized (cython) burst search loaded. - Optimized (cython) photon counting loaded. ------------------------------------------------------------- You are running FRETBursts (version 0.4rc4-20-g4ea8dd4). If you use this software in a publication, please cite it as: FRETBursts - An opensource single-molecule FRET burst analysis toolkit. A. Ingargiola 2014. https://github.com/tritemio/FRETBursts -------------------------------------------------------------
import os
%matplotlib inline
import fretbursts.style
from IPython.html.widgets import interact, interactive, fixed
from IPython.html import widgets
from IPython.display import display, display_png, display_svg, clear_output
from IPython.core.pylabtools import print_figure
import matplotlib as mpl
import seaborn as sns
sns.set_style('darkgrid')
fs = 13
rc={'font.size': fs, 'axes.labelsize': fs, 'legend.fontsize': fs,
'axes.titlesize': fs*1.1, 'xtick.labelsize': fs, 'ytick.labelsize': fs,
'savefig.dpi': 70,
#'figure.figsize': (6, 4),
#'axes.facecolor': '0.95', 'axes.edgecolor': '0.85', 'grid.color': '0.85',
#'axes.linewidth': 1,
}
sns.set(rc=rc)
blue = '#0055d4'
green = '#2ca02c'
red = "#e74c3c"
purple = "#9b59b6"
color_brewer = sns.color_palette("Set1", 9)
colors = np.array(color_brewer)[(1,0,2,3,4,8,6,7), :]
colors = list(colors)
colors[:3] = (blue, colors[1], green)
sns.set_palette(colors, 8)
#sns.palplot(sns.color_palette(colors, 8))
url = 'http://files.figshare.com/1839121/0023uLRpitc_NTP_20dT_0.5GndCl.hdf5'
download_file(url, save_dir='./data')
URL: http://files.figshare.com/1839121/0023uLRpitc_NTP_20dT_0.5GndCl.hdf5 File: 0023uLRpitc_NTP_20dT_0.5GndCl.hdf5 File already on disk: C:\Data\Antonio\software\src\fretbursts_notebooks\notebooks\data/0023uLRpitc_NTP_20dT_0.5GndCl.hdf5 Delete it to re-download.
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
'./data/0023uLRpitc_NTP_20dT_0.5GndCl.hdf5'
Let's check that the file exists:
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!
d = loader.hdf5(fname=full_fname)
d.add(det_donor_accept=(0, 1), alex_period=4000, D_ON=(2850, 580), A_ON=(900, 2580))
bpl.plot_alternation_hist(d)
loader.usalex_apply_period(d)
d.calc_bg(bg.exp_fit, time_s=1000, tail_min_us=(800, 4000, 1500, 1000, 3000))
#donor: 727351 #acceptor: 1579957 - Calculating BG rates ... [DONE]
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]
#dplot(ds, hist_fret)
#dplot(ds, hist_S);
#ds.E_fitter.fit_histogram(mfit.factory_three_gaussians())
We can make a simple E-S scatter plot with scatter_alex
:
dplot(ds, scatter_alex, figsize=(4,4), mew=1, ms=4, mec='black', color='purple');
We can also plot the ALEX histogram with a scatterplot overlay using hist2d_alex
:
dplot(ds, hist2d_alex);
def _alex_plot_style(g):
"""Set plot style and colorbar for an ALEX joint plot.
"""
g.set_axis_labels(xlabel="E", ylabel="S")
g.ax_marg_x.grid(True)
g.ax_marg_y.grid(True)
plt.setp(g.ax_marg_y.get_xticklabels(), visible=True)
plt.setp(g.ax_marg_x.get_yticklabels(), visible=True)
g.ax_marg_x.locator_params(axis='y', tight=True, nbins=3)
g.ax_marg_y.locator_params(axis='x', tight=True, nbins=3)
pos = g.ax_joint.get_position().get_points()
X, Y = pos[:, 0], pos[:, 1]
cax = plt.axes([1., Y[0], (X[1] - X[0])*0.045, Y[1]-Y[0]])
plt.colorbar(cax=cax)
def _hist_bursts(arr, dx, **kwargs):
"""Wrapper function for calling hist_burst_data() from seaborn plot_marginals().
"""
vertical = kwargs.get('vertical', False)
data_name = 'S' if vertical else 'E'
hist_burst_data(dx, data_name=data_name, **kwargs)
def _alex_hexbin_vmax(patches, vmax_fret=True, Smax=0.8):
"""Return max count E-S hexbin plot in `patches`.
When `vmax_fret` is True, returns the max count for S < Smax.
Otherwise returns the max count in all the histogram.
"""
counts = patches.get_array()
if vmax_fret:
offset = patches.get_offsets()
xoffset, yoffset = offset[:, 0], offset[:, 1]
mask = yoffset < Smax
vmax = counts[mask].max()
else:
vmax = counts.max()
return vmax
def _calc_vmin(vmax, vmax_threshold, vmin_default):
if vmax <= vmax_threshold:
vmin = vmin_default - 0.5*vmax
elif vmax_threshold < vmax < 2*vmax_threshold:
vmin = vmin_default - 0.5*vmax*((2*vmax_threshold - vmax)/(1.*vmax_threshold))
else:
vmin = vmin_default
return vmin
def alex_jointplot(d, i=0, gridsize=50, cmap='YlGnBu_crop', zero_transparent=True, vmax_fret=True,
vmax_threshold=10, vmin_default=0, vmin=None, cmap_compensate=False,
joint_kws = {},
marginal_kws = dict(show_kde=True, bandwidth=0.03, binwidth=0.03)):
"""Plot an ALEX join plot: E-S 2-D histograms and marginal plots.
"""
g = sns.JointGrid(x=d.E[i], y=d.S[i], ratio=3, space=0.2, xlim=(-0.2, 1.2), ylim=(-0.2, 1.2))
joint_kws.update(gridsize=gridsize, cmap=cmap, extent=(-0.2, 1.2, -0.2, 1.2))
jplot = g.plot_joint(plt.hexbin, mincnt=zero_transparent, **joint_kws)
poly = jplot.ax_joint.get_children()[2]
vmax = _alex_hexbin_vmax(poly, vmax_fret=vmax_fret)
if vmin is None:
if cmap_compensate:
#vmin = _calc_vmin(vmax, vmax_threshold, vmin_default)
vmin = -vmax/3
else:
vmin = vmin_default
poly.set_clim(vmin, vmax)
g.plot_marginals(_hist_bursts, dx=d, **marginal_kws)
g.annotate(lambda x, y: x.size, stat='# Bursts', template='{stat}: {val}')
_alex_plot_style(g)
As an example let's create a cutom colormap and plot a plain 2D ALEX histogram:
sns.palplot(sns.color_palette('nipy_spectral', 64))
c = sns.color_palette('nipy_spectral', 64)[2:43]
sns.palplot(c)
cmap = mpl.colors.LinearSegmentedColormap.from_list('alex', c)
cmap.set_under(alpha=0)
mpl.cm.register_cmap(name='alex_lv', cmap=cmap)
cmap
<matplotlib.colors.LinearSegmentedColormap at 0x1a58a780>
squarehist, _, _ = np.histogram2d(ds.E_, ds.S_, range=((-0.2, 1.2), (-0.2, 1.2)), bins=np.arange(-0.2, 1.2, 0.018))
im = plt.imshow(squarehist[:,::-1].T, extent=(-0.2, 1.2, -0.2, 1.2), vmin=1, cmap='alex_lv', zorder=10)
gca().set_axisbelow(True)
plt.colorbar()
<matplotlib.colorbar.Colorbar instance at 0x000000001A098F88>
sns.palplot(sns.color_palette('YlGnBu', 64))
c = sns.color_palette('YlGnBu', 64)[16:]
sns.palplot(c)
cmap = mpl.colors.LinearSegmentedColormap.from_list('YlGnBu_crop', c)
cmap.set_under(alpha=0)
mpl.cm.register_cmap(name='YlGnBu_crop', cmap=cmap)
cmap
<matplotlib.colors.LinearSegmentedColormap at 0x1b7823c8>
# 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())
<fretbursts.mfit.MultiFitter at 0x1ba605c0>
sns.set_style('whitegrid')
alex_jointplot(ds)
alex_jointplot(ds, cmap = 'GnBu',
cmap_compensate = True,
gridsize = 50,
zero_transparent = True,
vmax_fret = True,
joint_kws = dict(edgecolor='grey'))
sns.set_style('darkgrid')
alex_jointplot(ds)
alex_jointplot(ds,
#cmap = 'GnBu',
#cmap_compensate = True,
gridsize = 50,
zero_transparent = True,
vmax_fret = True,
joint_kws = dict(edgecolor='grey'))
sns.set_style('darkgrid')
sns.set_style('whitegrid')
@interact(cmap_compensate = False,
overlay = widgets.RadioButtonsWidget(values=['fit model', 'KDE']),
binwidth = widgets.FloatTextWidget(value=0.03, min=0.01, max=1),
bandwidth = widgets.FloatTextWidget(value=0.03, min=0.01, max=1),
gridsize = (10, 100),
min_size=(10, 500, 5),
cmap=widgets.DropdownWidget(value='YlGnBu_crop',
values=['YlGnBu_crop', 'YlOrRd', 'Blues', 'PuBuGn',
'PuBu', 'GnBu', 'YlGnBu', 'afmhot', 'alex_lv',
'copper', 'summer', 'winter', 'cubehelix']),
zero_transparent = True,
reverse_cmap = False,
vmax_fret = True,
)
def plot_(min_size=50, cmap_compensate=False, overlay='KDE', binwidth=0.03, bandwidth=0.03, gridsize=50, cmap='YlGnBu_crop',
reverse_cmap=False, vmax_fret=True, zero_transparent=True):
dx = Sel(d, select_bursts.size, add_naa=True, th1=min_size)
bext.bursts_fitter(dx, 'E', binwidth=binwidth, bandwidth=bandwidth,
model=mfit.factory_three_gaussians())
bext.bursts_fitter(dx, 'S', binwidth=binwidth, bandwidth=bandwidth,
model=mfit.factory_two_gaussians())
if reverse_cmap: cmap += '_r'
if binwidth < 0.01: binwidth = 0.01
if bandwidth < 0.01: bandwidth = 0.01
if overlay == 'fit model':
marginal_kws = dict(binwidth=binwidth, show_model=True)
else:
marginal_kws = dict(binwidth=binwidth, show_kde=True, bandwidth=bandwidth)
alex_jointplot(dx, cmap=cmap, gridsize=gridsize, zero_transparent=zero_transparent, vmax_fret=vmax_fret,
cmap_compensate=cmap_compensate, marginal_kws=marginal_kws, joint_kws=dict())#edgecolor='grey'))
fig = gcf()
plt.close()
display(fig)
from IPython.core.display import HTML
HTML(open("./styles/custom2.css", "r").read())