Example - FRET histogram fitting

This notebook is part of smFRET burst analysis software FRETBursts.

In this notebook shows how to fit a FRET histogram. For a complete tutorial on burst analysis see FRETBursts - us-ALEX smFRET burst analysis.



In [ ]:

    
from fretbursts import *



In [ ]:

    
sns = init_notebook(apionly=True)
import lmfit
print('lmfit version:', lmfit.__version__)



In [ ]:

    
# Tweak here matplotlib style
import matplotlib as mpl
mpl.rcParams['font.sans-serif'].insert(0, 'Arial')
mpl.rcParams['font.size'] = 12
%config InlineBackend.figure_format = 'retina'

Get and process data



In [ ]:

    
url = 'http://files.figshare.com/2182601/0023uLRpitc_NTP_20dT_0.5GndCl.hdf5'
download_file(url, save_dir='./data')
full_fname = "./data/0023uLRpitc_NTP_20dT_0.5GndCl.hdf5"

d = loader.photon_hdf5(full_fname)
loader.alex_apply_period(d)
d.calc_bg(bg.exp_fit, time_s=1000, tail_min_us=(800, 4000, 1500, 1000, 3000))
d.burst_search(L=10, m=10, F=6)
ds = d.select_bursts(select_bursts.size, add_naa=True, th1=30)

Fitting the FRET histogram

We start defining the model. Here we choose a 3-Gaussian model:



In [ ]:

    
model = mfit.factory_three_gaussians()
model.print_param_hints()

The previsou cell prints all the model parameters. Each parameters has an initial value and bounds (min, max). The column vary tells if a parameter is varied during the fit (if False the parameter is fixed). Parameters with an expression (Expr column) are not free but the are computed as a function of other parameters.

We can modify the paramenters constrains as follows:



In [ ]:

    
model.set_param_hint('p1_center', value=0.1, min=-0.1, max=0.3)
model.set_param_hint('p2_center', value=0.4, min=0.3, max=0.7)
model.set_param_hint('p2_sigma', value=0.04, min=0.02, max=0.18)
model.set_param_hint('p3_center', value=0.85, min=0.7, max=1.1)

Then, we fit and plot the model:



In [ ]:

    
E_fitter = bext.bursts_fitter(ds, 'E', binwidth=0.03)
E_fitter.fit_histogram(model=model, pdf=False, method='nelder')
E_fitter.fit_histogram(model=model, pdf=False, method='leastsq')
dplot(ds, hist_fret, show_model=True, pdf=False);



In [ ]:

    
# dplot(ds, hist_fret, show_model=True, pdf=False, figsize=(6, 4.5));
# plt.xlim(-0.1, 1.1)
# plt.savefig('fret_hist_fit.png', bbox_inches='tight', dpi=200, transparent=False)

The results are in E_fitter:



In [ ]:

    
res = E_fitter.fit_res[0]
res.params.pretty_print()

To get a dictionary of values:



In [ ]:

    
res.values

This is the startndard lmfit's fit report:



In [ ]:

    
print(res.fit_report(min_correl=0.5))

The previous cell reports error ranges computed from the covariance matrix. More accurare confidence intervals can be obtained with:



In [ ]:

    
ci = res.conf_interval()
lmfit.report_ci(ci)

Tidy fit results

It is convenient to put the fit results in a DataFrame for further analysis. A dataframe of fitted parameters is already in E_fitter:



In [ ]:

    
E_fitter.params

With pybroom we can get a "tidy" DataFrame with more complete fit results:



In [ ]:

    
import pybroom as br



In [ ]:

    
df = br.tidy(res)
df

Now, for example, we can easily select parameters by name:



In [ ]:

    
df.loc[df.name.str.contains('center')]



In [ ]:

    
df.loc[df.name.str.contains('sigma')]