FRETBursts - 8-spot smFRET burst analysis

This notebook is part of a tutorial series for the FRETBursts burst analysis software.

For a step-by-step introduction to FRETBursts usage please refer to us-ALEX smFRET burst analysis.

In this notebook we present a typical FRETBursts workflow for multi-spot smFRET burst analysis. Briefly, we show how to perform background estimation, burst search, burst selection, FRET histograms, and FRET efficiency fit using different methods.

Loading the software



In [ ]:

    
from fretbursts import *



In [ ]:

    
sns = init_notebook()



In [ ]:

    
import lmfit; lmfit.__version__



In [ ]:

    
import phconvert; phconvert.__version__

Downloading the sample data file

The complete example dataset can be downloaded from here.

Here we download an 8-spot smFRET measurement file using the download_file in FRETBursts:



In [ ]:

    
url = 'http://files.figshare.com/2182604/12d_New_30p_320mW_steer_3.hdf5'



In [ ]:

    
download_file(url, save_dir='./data')

Selecting a data file



In [ ]:

    
filename = "./data/12d_New_30p_320mW_steer_3.hdf5"



In [ ]:

    
import os
assert os.path.exists(filename)

Data load and Burst search

Load and process the data:



In [ ]:

    
d = loader.photon_hdf5(filename)

For convenience we can set the correction coefficients right away so that they will be used in the subsequent analysis. The correction coefficients are:

leakage or bleed-through: leakage
direct excitation: dir_ex (ALEX-only)
gamma-factor gamma

The direct excitation cannot be applied to non-ALEX (single-laser) smFRET measurements (like the current one).



In [ ]:

    
d.leakage = 0.038
d.gamma = 0.43

NOTE: at any later moment, after burst search, a simple reassignment of these coefficient will update the burst data with the new correction values.

Compute background and burst search:



In [ ]:

    
d.calc_bg(bg.exp_fit, time_s=30, tail_min_us='auto', F_bg=1.7)
d.burst_search(L=10, m=10, F=7)

Perform a background plot as a function of the channel:



In [ ]:

    
mch_plot_bg(d)

Let's take a look at the photon waiting times histograms and at the fitted background rates:



In [ ]:

    
dplot(d, hist_bg);

Using dplot exactly in the same way as for the single-spot data has now generated 8 subplots, one for each channel.

Let's plot a timetrace for the background to see is there are significant variations during the measurement:



In [ ]:

    
dplot(d, timetrace_bg);

We can look at the timetrace of the photon stream (binning):



In [ ]:

    
dplot(d, timetrace)
xlim(2, 3); ylim(-100, 100);

We can also open the same plot in an interactive window that allows scrolling (uncomment the following lines):



In [ ]:

    
#%matplotlib qt



In [ ]:

    
#dplot(d, timetrace, scroll=True);



In [ ]:

    
#ylim(-100, 100)



In [ ]:

    
#%matplotlib inline

Burst selection and FRET

Selecting bursts by burst size (`select_bursts.size`)



In [ ]:

    
gamma = d.gamma
gamma



In [ ]:

    
d.gamma = 1
ds = d.select_bursts(select_bursts.size, th1=30, gamma=1)
dplot(ds, hist_fret);



In [ ]:

    
ds = d.select_bursts(select_bursts.size, th1=25, gamma=gamma, donor_ref=False)
dplot(ds, hist_fret);



In [ ]:

    
ds = d.select_bursts(select_bursts.size, th1=25, gamma=gamma)
dplot(ds, hist_fret, weights='size', gamma=gamma);



In [ ]:

    
dplot(ds, scatter_fret_nd_na); ylim(0,200);

FRET Fitting

2-Gaussian mixture

Let's fit the $E$ histogram with a 2-Gaussians model:



In [ ]:

    
ds.gamma = 1.
bext.bursts_fitter(ds, weights=None)
ds.E_fitter.fit_histogram(mfit.factory_two_gaussians(), verbose=False)

The fitted parameters are stored in a pandas DataFrame:



In [ ]:

    
ds.E_fitter.params



In [ ]:

    
dplot(ds, hist_fret, weights=None, show_model=True,
      show_fit_stats=True, fit_from='p2_center');

Weighted Expectation Maximization

The expectation maximization (EM) method is particularly suited to resolve population mixtures. Note that the EM algorithm does not fit the histogram but the $E$ distribution with no binning.

FRETBursts include a weighted version of the EM algorithm that can take into account the burst size. The algorithm and benchmarks with the 2-Gaussian histogram fit are reported here.

You can find the EM algorithm in fretbursts/fit/gaussian_fit.py or typing:

bl.two_gaussian_fit_EM??



In [ ]:

    
# bl.two_gaussian_fit_EM??



In [ ]:

    
EM_results = ds.fit_E_two_gauss_EM(weights=None, gamma=1.)
EM_results

The fitted parameters for each channel are stored in the fit_E_res attribute:



In [ ]:

    
ds.fit_E_name, ds.fit_E_res

The model function is stored in:



In [ ]:

    
ds.fit_E_model

Let's plot the histogram and the model with parameters from the EM fit:



In [ ]:

    
AX = dplot(ds, hist_fret, weights=None)

x = np.r_[-0.2: 1.2 : 0.01]
for ich, (ax, E_fit) in enumerate(zip(AX.ravel(), EM_results)):
    ax.axvline(E_fit, ls='--', color='r')
    ax.plot(x, ds.fit_E_model(x, ds.fit_E_res[ich]))

print('E mean: %.2f%%   E delta: %.2f%%' %\
      (EM_results.mean()*100, (EM_results.max() - EM_results.min())*100))

Comparing 2-Gaussian and EM fit

To quickly compare the 2-Gaussians with the EM fit we convert the EM fit results in a DataFrame:



In [ ]:

    
import pandas as pd



In [ ]:

    
EM_results = pd.DataFrame(ds.fit_E_res, columns=['p1_center', 'p1_sigma', 'p2_center', 'p2_sigma', 'p1_amplitude'])
EM_results * 100



In [ ]:

    
ds.E_fitter.params * 100

And we compute the difference between the two sets of parameters:



In [ ]:

    
(ds.E_fitter.params - EM_results) * 100

NOTE: The EM method follows more the "asymmetry" of the peaks because the center is a weighted mean of the bursts. On the contrary the 2-Gaussians histogram fit tends to follows more the peak position an less the "asymmetric" tails.



In [ ]: