This notebook is part of a tutorial series for the FRETBursts burst analysis software.
In this notebook, we present a typical FRETBursts workflow for μs-ALEX smFRET burst analysis. Briefly, we show how to perform background estimation, burst search, burst selection, compute FRET histograms and ALEX histograms, do sub-population selection and finally, FRET efficiency fit.
Before running the notebook, you can click on menu Cell -> All Output -> Clear to clear all previous output. This will avoid mixing output from current execution and the previously saved one.
We start by loading FRETBursts
:
from fretbursts import *
- Optimized (cython) burst search loaded. - Optimized (cython) photon counting loaded. -------------------------------------------------------------- You are running FRETBursts (version 0.7+46.ge31fadb.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 --------------------------------------------------------------
Note that FRETBursts version string tells you the exact FRETBursts version (and revision) in use. Storing the version in the notebook helps with reproducibility and tracking software regressions.
Next, we initialize the default plot style for the notebook (under the hood it uses seaborn):
sns = init_notebook()
Note that the previous command has no output. Finally, we print the version of some dependencies:
import lmfit; lmfit.__version__
'1.0.3'
import phconvert; phconvert.__version__
'0.9'
The full list of smFRET measurements used in the FRETBursts tutorials can be found on Figshare.
This is the file we will download:
url = 'http://files.figshare.com/2182601/0023uLRpitc_NTP_20dT_0.5GndCl.hdf5'
url
variable above to download your own data file.
This is useful if you are executing FRETBursts online and you want to use
your own data file. See
First Steps.
Here, we download the data file and put it in a folder named data
,
inside the notebook folder:
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: /home/paul/Disk/Python/OpenSMFS/FRETBursts_notebooks/notebooks/data/0023uLRpitc_NTP_20dT_0.5GndCl.hdf5 Delete it to re-download.
NOTE: If you modified the
url
variable providing an invalid URL the previous command will fail. In this case, edit the cell containing theurl
and re-try the download.
Use one of the following 2 methods to select a data file.
Here, we can directly define the file name to be loaded:
filename = "./data/0023uLRpitc_NTP_20dT_0.5GndCl.hdf5"
filename
'./data/0023uLRpitc_NTP_20dT_0.5GndCl.hdf5'
Now filename
contains the path of the file you just selected.
Run again the previous cell to select a new file. In a following cell
we will check if the file actually exists.
Alternatively, you can select a data file with an "Open File" windows. Note that, since this only works in a local installation, the next commands are commented (so nothing will happen when running the cell).
If you want to try the File Dialog, you need to remove the #
signs:
# filename = OpenFileDialog()
# filename
Let's check that the file exists:
import os
if os.path.isfile(filename):
print("Perfect, file found!")
else:
print("Sorry, file:\n%s not found" % filename)
Perfect, file found!
We can finally load the data and store it in a variable called d
:
d = loader.photon_hdf5(filename)
If you don't get any message, the file is loaded successfully.
We can also set the 3 correction coefficients:
leakage
dir_ex
(ALEX-only)gamma
d.leakage = 0.11
d.dir_ex = 0.04
d.gamma = 1.
NOTE: at any later moment, after burst search, a simple reassignment of these coefficient will update the burst data with the new correction values.
At this point, timestamps and detectors numbers are contained in the ph_times_t
and det_t
attributes of d
. Let's print them:
d.ph_times_t, d.det_t
([array([ 146847, 188045, 294124, ..., 47999863658, 47999877783, 47999955353])], [array([0, 1, 1, ..., 1, 1, 0], dtype=uint32)])
We need to define some ALEX parameters:
d.add(det_donor_accept = (0, 1),
alex_period = 4000,
offset = 700,
D_ON = (2180, 3900),
A_ON = (200, 1800))
Here the parameters are:
det_donor_accept
: donor and acceptor channelsalex_period
: length of excitation period (in timestamps units)D_ON
and A_ON
: donor and acceptor excitation windowsoffset
: the offset between the start of alternation and start of timestamping
(see also Definition of alternation periods).To check that the above parameters are correct, we need to plot the histogram of timestamps (modulo the alternation period) and superimpose the two excitation period definitions to it:
bpl.plot_alternation_hist(d)
If the previous alternation histogram looks correct, the corresponding definitions of the excitation periods can be applied to the data using the following command:
loader.alex_apply_period(d)
# 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
If the previous histogram does not look right, the parameters in the d.add(...)
cell can be modified and checked by running the histogram plot cell until everything looks fine. Don't forget to apply the
parameters with loader.usalex_apply_period(d)
as a last step.
NOTE: After applying the ALEX parameters a new array of timestamps containing only photons inside the excitation periods is created (name
d.ph_times_m
). To save memory, by default, the old timestamps array (d.ph_times_t
) is deleted. Therefore, in the following, when we talk about all-photon selection we always refer to all photons inside both excitation periods.
The entire measurement data is now stored in the variable d
. Printing it
will give a compact representation containing the file-name and additional parameters:
d
data_0023uLRpitc_NTP_20dT_0.5GndCl G1.000 Lk11.000 dir4.0
To check the measurement duration (in seconds) run:
d.time_max
599.9994331624999
In this section basic concepts of plotting with FRETBursts using the timetrace plot as an example.
To plot a timetrace of the measurement we use:
dplot(d, timetrace);
Here, dplot
is a generic wrapper (the same for all plots)
that takes care of setting up the figure, title and axis
(in the multispot case dplot
creates multi-panel plot).
The second argument, timetrace
, is the actual plot function.
All the eventual additional arguments passed to dplot
are,
in turn, passed to the plot function (e.g. timetrace
).
If we look at the documentation for timetrace
function we notice that it accepts a long list of arguments.
In python, when an argument is not specified, it will take the default
value specified in the function definition (see previus link).
As an example, to change the bin size (i.e. duration) of the timetrace histogram,
we can look up in the timetrace
documentation
and find that the argument we need to modify is binwidth
(we can also see that the default value is 0.001
seconds).
We can then re-plot the timetrace using a bin of 0.5 ms:
dplot(d, timetrace, binwidth=0.5e-3);
The timetrace is computed between tmin
and tmax
(by default 0 and 200s),
but as you can see is displayed only between 0 an 1 second, just because these
are the default x-axis limits. The axis limits can be changes by using the
standard matplotlib command plt.xlim()
.
On the other hand, to change the range where the timetrace is computed,
we pass the additional arguments tmin
and tmax
as follows:
dplot(d, timetrace, binwidth=0.5e-3, tmin=50, tmax=150)
plt.xlim(51, 52);
When using FRETBursts in a notebook, all plots are static by default.
This is because we use the so called inline
backend of matplotlib.
If you want to manipulate figures interactively, you can switch
to the interactive notebook
backend with:
%matplotlib notebook
to go back to inline use:
%matplotlib inline
NOTE: Currently, the notebook
backend is incompatible with the QT backend.
If in a session you activate the notebook
backend, then switching to the QT backend requires
restarting the notebook. Conversely, you can switch between inline
and notebook
or between inline
and qt4
backends in the same session wihtou issues.
See also:¶
- bpl.timetrace function documentation
- bpl.ratetrace function documentation
- Intensity timetrace and Rate-timetrace, a later section in this tutorial.
As a first step of the analysis, we need to estimate the background. Here we will compute the background using the recommended approach of using the auto-threshold.
For more details see [Background estimation](Example - Background estimation.ipynb).
It is a good practice to monitor background rates as a function of time. Here, we compute background in adjacent 30s windows (called background periods) and plot the estimated rates as a function of time.
d.calc_bg(bg.exp_fit, time_s=30, tail_min_us='auto')
- Calculating BG rates ... Channel 0 [DONE]
dplot(d, timetrace_bg)
<AxesSubplot:title={'center':'data_0023uLRpitc_NTP_20dT_0.5GndCl'}, xlabel='Time (s)', ylabel='BG rate (kcps)'>
NOTE: All background data is stored in
d.bg
. For details on how to to export it see the [Background estimation](Example - Background estimation.ipynb) notebook.
The first step of burst analysis is the burst search.
We will use the sliding-window algorithm on all photons. Note that "all photons", as mentioned before, means all photons selected in the alternation histogram. An important variation compared to the classical sliding-windows is that the threshold-rate for burst start is computed as a function of the background and changes when the background changes during the measurement.
To perform a burst search evaluating the photon rate with
10 photons (m=10
), and selecting a minimum rate 6 times larger than
the background rate (F=6) calculated with all photons (default):
d.burst_search(L=10, m=10, F=6)
- 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]
The previous command performs the burst search, corrects the bursts sizes for background, spectral leakage and direct excitation, and computes $\gamma$-corrected FRET and Stoichiometry.
See the
burst_search
documentation for more details.
We can plot the resulting FRET histogram using the following command:
dplot(d, hist_fret);
All pre-defined plots follow this pattern:
call the generic dplot()
function, passing 2 parameters:
d
in this case)hist_fret
)In some case we can add other optional parameters to tweak the plot.
All plot functions start with hist_
for histograms,
scatter_
for scatter-plots or timetrace_
for plots as a function
of measurement time. You can use autocompletion to find all
plot function or you can look in bursts_plot.py
where
all plot functions are defined.
Instead of hist_fret
we can use hist_fret_kde
to add a KDE overlay. Also, we can plot a weighted histogram by passing an additional parameter weights
:
dplot(d, hist_fret, show_kde=True);
dplot(d, hist_fret, show_kde=True, weights='size');
- Overwriting the old E_fitter object with the new weights.
You can experiment with different weighting schema (for all
supported weights see get_weigths()
function in fret_fit.py
).
When we performed the burst search, we specified L=10
without
explaining what this parameter means. L is traditionally the minimum size
(number of photons) for a burst: smaller bursts will be rejected.
By setting L=m (10 in this case) we are deciding to not discard
any burst (because the smallest detected burst has at least m counts).
Selecting the bursts in a second step, by applying a minimum burst size criterion, results in a more accurate and unbiased selection.
For example, we can select bursts with more than 30 photons (after
background, gamma, leakage and direct excitation corrections)
and store the result in a new
Data()
variable ds
:
ds = d.select_bursts(select_bursts.size, th1=30)
By defaults the burst size includes donor and acceptor photons
during donor excitation. To add acceptor photons during
acceptor excitation (naa
), we add the parameter add_naa=True
:
ds = d.select_bursts(select_bursts.size, add_naa=True, th1=30)
Similar to plot functions, all selection functions
are defined in select_bursts.py
and you can access them by typing
select_bursts.
and using the TAB key for autocompletion.
See also:
- Burst selection in the documentation.
In particular the function
select_bursts.size
andData.select_bursts
.
To replot the FRET histogram after selection (note that now
we are passing ds
to the plot function):
dplot(ds, hist_fret);
Note how the histogram exhibits much more clearly defined peaks after burst selection.
Under the hood the previous hist_fret
plot creates a MultiFitter
object for $E$ values. This object, stored as ds.E_fitter
, operates
on multi-channel data and computes the histogram, KDE and can fit
the histogram with a model (lmfit.Model).
Now, just for illustration purposes, we fit the previous histogram with 3 Gaussians, using the already created ds.E_fitter
object:
ds.E_fitter.fit_histogram(mfit.factory_three_gaussians(), verbose=False)
dplot(ds, hist_fret, show_model=True);
The bin width can be changed with binwidth
argument. Alternatively,
an arbitrary array of bin edges can be passed in bins
(overriding binwidth
).
We can customize the appearance of this plot (type
hist_fret?
for the complete set of arguments).
For example to change from a bar plot to a line-plot
we use the hist_style
argument:
dplot(ds, hist_fret, show_model=True, hist_style='line')
<AxesSubplot:title={'center':'data_0023uLRpitc_NTP_20dT_0.5GndCl, T=759μs, #bu=5409'}, xlabel='E', ylabel='PDF'>
We can customize the line-plot, bar-plot, the model plot and the KDE plot by passing dictionaries with matplotlib style. The name of the arguments are:
hist_plot_style
: style for the histogram line-plothist_bar_style
: style for the histogram bar-plotmodel_plot_style
: style for the model plotkde_plot_style
: style for the KDE plotAs an example:
dplot(ds, hist_fret, show_model=True, hist_style='bar', show_kde=True,
kde_plot_style = dict(linewidth=5, color='orange', alpha=0.6),
hist_plot_style = dict(linewidth=3, markersize=8, color='b', alpha=0.6))
plt.legend();
Similarly, we can plot the burst size using all photons
(type hist_size?
to learn about all plot options):
dplot(ds, hist_size, add_naa=True);
Or plot the burst size histogram for the different components:
dplot(ds, hist_size_all);
NOTE: The previous plot may generate a benign warning due to the presence of zeroes when switching to log scale. Just ignore it.
A scatterplot of Size vs FRET is created by:
dplot(ds, scatter_fret_nd_na)
xlim(-1, 2)
(-1.0, 2.0)
We can further select only bursts smaller than 300 photons to get rid of multi-molecule events:
ds2 = ds.select_bursts(select_bursts.size, th2=300)
and superimpose the two histograms before and after selection to see the difference:
ax = dplot(ds2, hist_fret, hist_style='bar', show_kde=True,
hist_bar_style = dict(facecolor='r', alpha=0.5,
label='Hist. no large bursts'),
kde_plot_style = dict(lw=3, color='m',
label='KDE no large bursts'))
dplot(ds, hist_fret, ax=ax, hist_style='bar', show_kde=True,
hist_bar_style = dict(label='Hist. with large bursts'),
kde_plot_style = dict(lw=3, label='KDE with large bursts'))
plt.legend();
NOTE: It is not necessarily true that bursts with more that 300 photons represents multiple molecules. To asses the valididty of this assumption it can be useful to plot the peak count rates in each burst. See
hist_burst_phrate
for this kind of plot.
We can find the KDE peak position in a range (let say 0.2 ... 0.6):
ds.E_fitter.find_kde_max(np.r_[0:1:0.0002], xmin=0.2, xmax=0.6)
and plot it with show_kde_peak=True
, we also use show_fit_value=True
to show a box with the fitted value:
dplot(ds, hist_fret, hist_style='line',
show_fit_value=True,
show_kde=True, show_kde_peak=True);
Instead of using the KDE, we can use the peak position as fitted from a gaussian model.
ds.E_fitter.fit_histogram(mfit.factory_three_gaussians(), verbose=False)
To select which peak to show we use fit_from='p1_center'
:
dplot(ds, hist_fret, hist_style='line',
show_fit_value=True, fit_from='p2_center',
show_model=True);
The string 'p2_center'
is the name of the parameter of the
gaussian fit that we want to show in the text box. To see all
the parameters of the model we look in:
ds.E_fitter.params # <-- pandas DataFrame, one row per channel
p1_amplitude | p1_center | p1_fwhm | p1_height | p1_sigma | p2_amplitude | p2_center | p2_fwhm | p2_height | p2_sigma | p3_amplitude | p3_center | p3_fwhm | p3_height | p3_sigma | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 0.485732 | 0.177278 | 0.470927 | 0.968971 | 0.199984 | 0.321084 | 0.457707 | 0.351419 | 0.858344 | 0.149234 | 0.218839 | 0.838622 | 0.254492 | 0.807827 | 0.108073 |
We can create 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);
Finally we can also plot an ALEX histogram and marginals (joint plots) as follow (for more options see: [Example - usALEX histogram](Example - usALEX histogram.ipynb)):
alex_jointplot(ds);
<class 'matplotlib.figure.Figure'>
To get rid of the large donor-only population, we can simply select bursts with at least 10 photons in the acceptor channel (during acceptor excitation). At the same time, with a burst size selection using Dex photons we get rid of the A-only population:
ds = d.select_bursts(select_bursts.size, th1=20)
ds2 = ds.select_bursts(select_bursts.naa, th1=10)
alex_jointplot(ds2);
<class 'matplotlib.figure.Figure'>
As you can see, we remained with only the FRET sub-populations.
See next sections show how to select a region on the E/S histogram to isolate a subpopulation.
#
) from the lines
containing the %matplotlib
command.
To select bursts graphically, we need to open the ALEX histogram in a new (QT) window, drag the mouse to define a selection and have it printed here in the notebook.
Here we describe how to do it in 3 steps.
Step 1 Switch the plot modality to QT (i.e. opens graphs in a separate window):
# Switches to open plot in external window
#%matplotlib qt
Step 2 Plot the E-S histogram in an external windows. There you can use the mouse to select a region:
# ALEX histogram with GUI selection enabled
dplot(ds, hist2d_alex, gui_sel=True)
<AxesSubplot:title={'center':'data_0023uLRpitc_NTP_20dT_0.5GndCl, T=759μs, #bu=4255'}, xlabel='E', ylabel='S'>
The region selection is printed above.
Step 3 Restore the normal inline plotting (no more external windows).
# Switch back to show plots inline in the notebook
#%matplotlib inline
To apply a selection based on E-S values, we can paste the values obtained from the previous plot (or we can type them in manually).
The selection function used here is select_bursts.ES
.
The same E and S boundaries can define either a rectangular
or an elliptical selection by using
respectively rect=True
or rect=False
.
Here we use the elliptical selection:
roi = dict(E1=-0.07, E2=1.17, S1=0.18, S2=0.70, rect=False)
d_fret_mix = ds.select_bursts(select_bursts.ES, **roi)
By plotting the FRET histogram we can double check that the selection has been applied:
g = alex_jointplot(d_fret_mix)
bpl.plot_ES_selection(g.axes[0], **roi);
<class 'matplotlib.figure.Figure'>
Now we can further separate high- and low FRET sub-populations.
roi_high_fret = dict(E1=0.65, E2=1.09, S1=-0.13, S2=0.96)
d_high_fret = d_fret_mix.select_bursts(select_bursts.ES, **roi_high_fret)
roi_low_fret = dict(E1=-0.19, E2=0.64, S1=-0.05, S2=0.92)
d_low_fret = d_fret_mix.select_bursts(select_bursts.ES, **roi_low_fret)
alex_jointplot(d_high_fret);
alex_jointplot(d_low_fret);
<class 'matplotlib.figure.Figure'> <class 'matplotlib.figure.Figure'>
We can for example compute the ratio of high- and low-fret bursts:
d_low_fret.num_bursts / d_high_fret.num_bursts
array([2.70819672])
To plot a burst-width histogram, we use hist_width
instead of hist_fret
:
dplot(d_low_fret, hist_width)
<AxesSubplot:title={'center':'data_0023uLRpitc_NTP_20dT_0.5GndCl, T=759μs, #bu=1652'}, xlabel='Burst width (ms)', ylabel='PDF'>
To use a larger bin size, plots two sub-populations (in different color) and add a legend:
ax = dplot(d_high_fret, hist_width, bins=(0, 10, 0.2))
dplot(d_low_fret, hist_width, bins=(0, 10, 0.2), ax=ax)
plt.legend(['High-FRET population', 'Low-FRET population']);
Finally, we compute the mean burst width for each subpopulation:
millisec = d.clk_p * 1e3
mean_b_width_low_fret = d_low_fret.mburst[0].width.mean() * millisec
mean_b_width_high_fret = d_high_fret.mburst[0].width.mean() * millisec
print('Mean burst width: %.1f ms (high-fret), %.1f (low-fret)' %
(mean_b_width_high_fret, mean_b_width_low_fret))
Mean burst width: 2.7 ms (high-fret), 3.5 (low-fret)
We can fit a FRET distribution to any model. For example,
we will fit the FRET selection (2 FRET sub-populations) with 2 Gaussians.
The fitting model is a Model
object
from the lmfit
library.
The first step, previously performed implicitely by the hist_fret()
plot function, is to create a MultiFitter
object for $E$ and/or $S$, by calling bext.bursts_fitter().
With MultiFitter
objects, we can compute histograms,
KDEs and fit the histogram in one single step, as in the following example:
E_fitter = bext.bursts_fitter(d_fret_mix, 'E', binwidth=0.03, bandwidth=0.03,
model=mfit.factory_two_gaussians())
S_fitter = bext.bursts_fitter(d_fret_mix, 'S', binwidth=0.03, bandwidth=0.03,
model=mfit.factory_gaussian())
However if we want to modify the model (for example to add a
constrain) we need to perform the fit in a second step.
To skip the fitting, we simply avoid passing a model
:
E_fitter = bext.bursts_fitter(d_fret_mix, 'E', binwidth=0.03, bandwidth=0.03)
S_fitter = bext.bursts_fitter(d_fret_mix, 'S', binwidth=0.03, bandwidth=0.03)
Now we create a model and initialize the parameters
using mfit.factory_two_gaussians()
(see Model factory functions
for a list of pre-defined models in FRETBursts):
model = mfit.factory_two_gaussians(add_bridge=True)
We can see the list of parameters, initial values and constraints:
model.print_param_hints()
Name Value Min Max Vary Expr br_amplitude 0.0001 0 inf True br_center1 nan 0 inf True p1_center br_center2 nan 0 inf True p2_center br_sigma1 nan 0 inf True p1_sigma br_sigma2 nan 0 inf True p2_sigma p1_amplitude 1 0.01 inf True p1_center 0.1 -1 2 True p1_fwhm nan -inf inf True 2.3548200*p1_sigma p1_height nan -inf inf True 0.3989423*p1_amplitude/max(1e-15, p1_sigma) p1_sigma 0.03 0.01 0.2 True p2_amplitude 1 0.01 inf True p2_center 0.9 -1 2 True p2_fwhm nan -inf inf True 2.3548200*p2_sigma p2_height nan -inf inf True 0.3989423*p2_amplitude/max(1e-15, p2_sigma) p2_sigma 0.03 0.01 0.2 True
We can change initial values and/or constrains (bounds):
model.set_param_hint('p1_center', value=0.3, min=-0.1, max=0.6)
model.set_param_hint('p2_center', value=0.85, min=0.5, max=1.1)
model.print_param_hints()
Name Value Min Max Vary Expr br_amplitude 0.0001 0 inf True br_center1 nan 0 inf True p1_center br_center2 nan 0 inf True p2_center br_sigma1 nan 0 inf True p1_sigma br_sigma2 nan 0 inf True p2_sigma p1_amplitude 1 0.01 inf True p1_center 0.3 -0.1 0.6 True p1_fwhm nan -inf inf True 2.3548200*p1_sigma p1_height nan -inf inf True 0.3989423*p1_amplitude/max(1e-15, p1_sigma) p1_sigma 0.03 0.01 0.2 True p2_amplitude 1 0.01 inf True p2_center 0.85 0.5 1.1 True p2_fwhm nan -inf inf True 2.3548200*p2_sigma p2_height nan -inf inf True 0.3989423*p2_amplitude/max(1e-15, p2_sigma) p2_sigma 0.03 0.01 0.2 True
Finally, we fit the histogram with one of the supported minimization methods (the default is least-squares):
E_fitter.fit_histogram(model=model, verbose=False) # default method is 'leastsq'
#E_fitter.fit_histogram(model=model, method='nelder') # example using simplex method
To plot the model with the fitted parameters on top of the FRET
histogram, we add show_model=True
as seen before:
dplot(d_fret_mix, hist_fret, show_model=True);
The fitting results of lmfit
models are stored in lmfit.ModelResult objects. In FRETBursts, these objects are available in MultiFitter.fit_res
:
E_fitter.fit_res[0]
fitting method | leastsq | |
# function evals | 106 | |
# data points | 46 | |
# variables | 7 | |
chi-square | 0.38773981 | |
reduced chi-square | 0.00994205 | |
Akaike info crit. | -205.698859 | |
Bayesian info crit. | -192.898369 |
name | value | standard error | relative error | initial value | min | max | vary | expression |
---|---|---|---|---|---|---|---|---|
p1_amplitude | 0.61479421 | 0.02047496 | (3.33%) | 1 | 0.01000000 | inf | True | |
p1_center | 0.35144517 | 0.00434333 | (1.24%) | 0.3 | -0.10000000 | 0.60000000 | True | |
p1_sigma | 0.11196388 | 0.00427400 | (3.82%) | 0.03 | 0.01000000 | 0.20000000 | True | |
p2_amplitude | 0.19882567 | 0.02129342 | (10.71%) | 1 | 0.01000000 | inf | True | |
p2_center | 0.83950345 | 0.00912688 | (1.09%) | 0.85 | 0.50000000 | 1.10000000 | True | |
p2_sigma | 0.07096411 | 0.00815099 | (11.49%) | 0.03 | 0.01000000 | 0.20000000 | True | |
br_center1 | 0.35144517 | 0.00434333 | (1.24%) | 0.3 | 0.00000000 | inf | False | p1_center |
br_center2 | 0.83950345 | 0.00912688 | (1.09%) | 0.85 | 0.00000000 | inf | False | p2_center |
br_sigma1 | 0.11196388 | 0.00427400 | (3.82%) | 0.03 | 0.00000000 | inf | False | p1_sigma |
br_sigma2 | 0.07096411 | 0.00815099 | (11.49%) | 0.03 | 0.00000000 | inf | False | p2_sigma |
br_amplitude | 0.66980918 | 0.06383628 | (9.53%) | 0.0001 | 0.00000000 | inf | True | |
p1_fwhm | 0.26365478 | 0.01006451 | (3.82%) | 0.0706446 | -inf | inf | False | 2.3548200*p1_sigma |
p1_height | 2.19059412 | 0.04737637 | (2.16%) | 13.298076666666669 | -inf | inf | False | 0.3989423*p1_amplitude/max(1e-15, p1_sigma) |
p2_fwhm | 0.16710769 | 0.01919412 | (11.49%) | 0.0706446 | -inf | inf | False | 2.3548200*p2_sigma |
p2_height | 1.11774778 | 0.05769748 | (5.16%) | 13.298076666666669 | -inf | inf | False | 0.3989423*p2_amplitude/max(1e-15, p2_sigma) |
p2_amplitude | p2_sigma | 0.8941 |
p1_amplitude | p1_sigma | 0.8254 |
p2_amplitude | p2_center | -0.7606 |
p2_center | p2_sigma | -0.7223 |
p1_amplitude | br_amplitude | -0.6419 |
p1_sigma | br_amplitude | -0.6097 |
p1_center | br_amplitude | -0.6089 |
p1_amplitude | p1_center | 0.5778 |
p1_center | p1_sigma | 0.5587 |
p2_amplitude | br_amplitude | -0.5581 |
p2_center | br_amplitude | 0.5345 |
p2_sigma | br_amplitude | -0.4969 |
p1_amplitude | p2_amplitude | 0.3311 |
p1_center | p2_amplitude | 0.3183 |
p1_amplitude | p2_center | -0.3179 |
p1_sigma | p2_amplitude | 0.3059 |
p1_center | p2_center | -0.3055 |
p1_sigma | p2_center | -0.2940 |
p1_amplitude | p2_sigma | 0.2920 |
p1_center | p2_sigma | 0.2811 |
p1_sigma | p2_sigma | 0.2691 |
Here the [0]
indicates CH=0. This index is used in multispot measurements, in which
there are several channels.
The values of the fitted parameters are available in the best_values
dictionary:
results = E_fitter.fit_res[0]
results.params.pretty_print(columns=['value'])
Name Value br_amplitude 0.6698 br_center1 0.3514 br_center2 0.8395 br_sigma1 0.112 br_sigma2 0.07096 p1_amplitude 0.6148 p1_center 0.3514 p1_fwhm 0.2637 p1_height 2.191 p1_sigma 0.112 p2_amplitude 0.1988 p2_center 0.8395 p2_fwhm 0.1671 p2_height 1.118 p2_sigma 0.07096
We can also print a complete report of fitted parameters including reduced $\chi^2$, error ranges ($\pm 1 \sigma$) and correlations:
print(results.fit_report(min_correl=0.5))
[[Model]] ((Model(gaussian, prefix='p1_') + Model(gaussian, prefix='p2_')) + Model(bridge_function, prefix='br_')) [[Fit Statistics]] # fitting method = leastsq # function evals = 106 # data points = 46 # variables = 7 chi-square = 0.38773981 reduced chi-square = 0.00994205 Akaike info crit = -205.698859 Bayesian info crit = -192.898369 [[Variables]] p1_amplitude: 0.61479421 +/- 0.02047496 (3.33%) (init = 1) p1_center: 0.35144517 +/- 0.00434333 (1.24%) (init = 0.3) p1_sigma: 0.11196388 +/- 0.00427400 (3.82%) (init = 0.03) p2_amplitude: 0.19882567 +/- 0.02129342 (10.71%) (init = 1) p2_center: 0.83950345 +/- 0.00912688 (1.09%) (init = 0.85) p2_sigma: 0.07096411 +/- 0.00815099 (11.49%) (init = 0.03) br_center1: 0.35144517 +/- 0.00434333 (1.24%) == 'p1_center' br_center2: 0.83950345 +/- 0.00912688 (1.09%) == 'p2_center' br_sigma1: 0.11196388 +/- 0.00427400 (3.82%) == 'p1_sigma' br_sigma2: 0.07096411 +/- 0.00815099 (11.49%) == 'p2_sigma' br_amplitude: 0.66980918 +/- 0.06383628 (9.53%) (init = 0.0001) p1_fwhm: 0.26365478 +/- 0.01006451 (3.82%) == '2.3548200*p1_sigma' p1_height: 2.19059412 +/- 0.04737637 (2.16%) == '0.3989423*p1_amplitude/max(1e-15, p1_sigma)' p2_fwhm: 0.16710769 +/- 0.01919412 (11.49%) == '2.3548200*p2_sigma' p2_height: 1.11774778 +/- 0.05769748 (5.16%) == '0.3989423*p2_amplitude/max(1e-15, p2_sigma)' [[Correlations]] (unreported correlations are < 0.500) C(p2_amplitude, p2_sigma) = 0.894 C(p1_amplitude, p1_sigma) = 0.825 C(p2_amplitude, p2_center) = -0.761 C(p2_center, p2_sigma) = -0.722 C(p1_amplitude, br_amplitude) = -0.642 C(p1_sigma, br_amplitude) = -0.610 C(p1_center, br_amplitude) = -0.609 C(p1_amplitude, p1_center) = 0.578 C(p1_center, p1_sigma) = 0.559 C(p2_amplitude, br_amplitude) = -0.558 C(p2_center, br_amplitude) = 0.535
To customize the printed report, see lmfit.fit_report() docs.
We can also take a look at the initial parameters:
results.init_params.pretty_print()
Name Value Min Max Stderr Vary Expr Brute_Step br_amplitude 0.0001 0 inf None True None None br_center1 0.3 0 inf None False p1_center None br_center2 0.85 0 inf None False p2_center None br_sigma1 0.03 0 inf None False p1_sigma None br_sigma2 0.03 0 inf None False p2_sigma None p1_amplitude 1 0.01 inf None True None None p1_center 0.3 -0.1 0.6 None True None None p1_fwhm 0.07064 -inf inf None False 2.3548200*p1_sigma None p1_height 13.3 -inf inf None False 0.3989423*p1_amplitude/max(1e-15, p1_sigma) None p1_sigma 0.03 0.01 0.2 None True None None p2_amplitude 1 0.01 inf None True None None p2_center 0.85 0.5 1.1 None True None None p2_fwhm 0.07064 -inf inf None False 2.3548200*p2_sigma None p2_height 13.3 -inf inf None False 0.3989423*p2_amplitude/max(1e-15, p2_sigma) None p2_sigma 0.03 0.01 0.2 None True None None
We can compute more accurate confidence-intervals (note that it can take several seconds, see lmfit docs for the method details):
ci = results.conf_interval()
lmfit.report_ci(ci)
99.73% 95.45% 68.27% _BEST_ 68.27% 95.45% 99.73% p1_amplitude: -0.06526 -0.04276 -0.02132 0.61479 +0.02218 +0.04639 +0.07413 p1_center : -0.01320 -0.00860 -0.00427 0.35145 +0.00441 +0.00919 +0.01465 p1_sigma : -0.01385 -0.00916 -0.00461 0.11196 +0.00487 +0.01027 +0.01653 p2_amplitude: -0.06152 -0.03973 -0.01977 0.19883 +0.02082 +0.04376 +0.07133 p2_center : -0.03011 -0.01845 -0.00887 0.83950 +0.00869 +0.01785 +0.02857 p2_sigma : -0.02165 -0.01454 -0.00747 0.07096 +0.00833 +0.01798 +0.03016 br_amplitude: -0.26887 -0.15099 -0.06769 0.66981 +0.05953 +0.11580 +0.17213
Finally, FRETBursts's
MultiFitter object
(e.g. E_fitter
) stores the fitted parameters in E_fitter.params
as pandas DataFrame:
E_fitter.params
br_amplitude | br_center1 | br_center2 | br_sigma1 | br_sigma2 | p1_amplitude | p1_center | p1_fwhm | p1_height | p1_sigma | p2_amplitude | p2_center | p2_fwhm | p2_height | p2_sigma | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 0.669809 | 0.351445 | 0.839503 | 0.111964 | 0.070964 | 0.614794 | 0.351445 | 0.263655 | 2.190594 | 0.111964 | 0.198826 | 0.839503 | 0.167108 | 1.117748 | 0.070964 |
For more information on fitting see:
To export data computed by FRETBursts, you need to find first where the data is stored or, for data computed on fly, which function/method is used to computed it.
Most data computed by FRETBursts is stored in the Data
object
(see the Data docs for the list of attributes).
For example:
Data.nd
, Data.na
, Data.naa
Data.ph_times_m
All these attributes are list of arrays, one per excitation spot.
This means that, even for single-spot data, we need to use indexing ([0]
)
to obtain the array. For example:
d.nd[0]
array([29.99522069, 5.6423473 , 7.56005931, ..., 8.57505512, 7.49520504, 2.47597537])
Any numpy array can be saved to disk with the method .tofile()
. For example:
d.nd[0].tofile('n_dd.csv', sep=',')
The background data is stored in Data.bg
(and Data.bg_mean
) attribute (see Getting the background rates in this notebook). This attribute is dict containing lists containing arrays (scalars) and ca be saved to disk similarly.
Part of burst data is stored in Data
attributes (mburst
, nd
, na
, naa
)
and part is computed on fly (duration, raw counts, etc.).
To put all burst data in a single "table" (a pandas.DataFrame
), one row
per burst, we can use bext.burst_data
:
bursts = bext.burst_data(ds)
bursts
size_raw | t_start | t_stop | width_ms | E | S | nd | na | nt | nda | naa | spot | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 37 | 0.062756 | 0.064491 | 1.734962 | -0.075550 | 0.931906 | 29.995221 | -2.106952 | 29.926057 | -0.135784 | 2.037789 | 0 |
1 | 53 | 0.203880 | 0.204703 | 0.823350 | 0.974436 | 0.409229 | 0.523168 | 19.941636 | 50.008174 | -0.064438 | 29.543370 | 0 |
2 | 43 | 0.375110 | 0.377019 | 1.908562 | 0.868662 | 0.580282 | 2.894683 | 19.145258 | 37.981451 | -0.149371 | 15.941510 | 0 |
3 | 93 | 0.528828 | 0.530547 | 1.719013 | 0.725235 | 0.338042 | 8.004458 | 21.127545 | 86.178638 | -0.134536 | 57.046635 | 0 |
4 | 35 | 0.739682 | 0.741413 | 1.730937 | 0.153470 | 0.859937 | 20.997552 | 3.806706 | 28.844279 | -0.135469 | 4.040021 | 0 |
... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
4250 | 104 | 599.544345 | 599.548232 | 3.887150 | 0.919265 | 0.371930 | 2.771193 | 31.553294 | 92.287514 | 0.705163 | 57.963027 | 0 |
4251 | 139 | 599.547832 | 599.551430 | 3.597612 | 0.731258 | 0.408432 | 13.937208 | 37.923749 | 126.975709 | -0.272876 | 75.114752 | 0 |
4252 | 127 | 599.570555 | 599.573195 | 2.639850 | 0.584503 | 0.257582 | 12.486368 | 17.565309 | 116.668323 | -0.200231 | 86.616646 | 0 |
4253 | 118 | 599.726445 | 599.729008 | 2.562925 | 0.798192 | 0.342427 | 7.530475 | 29.784508 | 108.971939 | -0.194396 | 71.656957 | 0 |
4254 | 41 | 599.767195 | 599.769260 | 2.064862 | 0.680829 | 0.691640 | 7.816053 | 16.672542 | 35.406551 | -0.156618 | 10.917955 | 0 |
4255 rows × 12 columns
NOTE: For multi-spot data use the
ich
argument to get burst data for the different spots.
This table (as any pandas.DataFrame
) can be saved in a CSV text file with:
bursts.to_csv('bursts.csv')
For an initial inspection of a data file, it is common to compute an intensity timetrace plot (or timetrace for short). We also can compute a similar plot called ratetrace, which does not bin counts but shows the instantaneous count rate. In both cases, it is convenient to scroll along the plot interactively.
In FRETBursts, the two timetrace and ratetrace plots support interactive scrolling. We just need to switch from the inline backend to QT, as we did before.
#
.
Here are the corresponding commands:
#%matplotlib qt
#dplot(ds, ratetrace, scroll=True, bursts=True)
#dplot(ds, timetrace, tmax=600, scroll=True, bursts=True)
#%matplotlib inline
Note that, to save memory, the previous plots stops will show the first 200 s
of measurement.
For both timetrace
and ratetrace
, you can extent (or reduce) the plotted
region passing the tmin
and tmax
arguments (values in seconds).
Related documentation:¶