tritemio / fretbursts

Burst analysis software for smFRET. **Moved to OpenSMFS organization**

Hi there 👋

• 🔭 Below you find some of my Open Source packages and contributions.

Thanks for visiting!

fretbursts's Issues

Lifetime Fitting

I was wondering if there is already some function implemented to calculate lifetimes from the TCSPC histogram?

Not getting same FRET hist with and without PIE

We measured the same sample using one excitation laser and PIE, and after analysis we obtain different FRET histograms. Shouldn't we obtain the same FRET histogram, and then with PIE having the ability to exclude the zero efficiency peak caused by molecules with only donor fluorophore?

In both cases we have used an "All photons" search:
d.burst_search(F = 6, m = 10, ph_sel = Ph_sel('all'))
ds = d.select_bursts(select_bursts.size, th1 = 30)

With PIE, we have also tried to select only donor and acceptor photons emitted during donor excitation, to correspond to the photons for non-ALEX data:
d.burst_search(F = 6, m = 10, ph_sel = Ph_sel(Dex=’DAem’))
but we still get the same result.

Do you have some insight into what might be happening?

API: Make leakage, gamma and dir_ex properties

The usage would be simplified if these 3 correction parameters would be properties.

In this way, an assignemt:

d.leakage = 0.05

would cause recorrecting all the burst data.

Questions:

• How to handle when we do not want to correct (like burst search with nofret=True)?
• Are there any method argumets that cannot be tweaked this way?

ymin and ymax deprication warning for Matplotlib 3.2

I am getting a persistent deprecation warning for many of the dplot() functions.

This is the error:

/Users/maya/anaconda3/envs/py36/lib/python3.6/site-packages/matplotlib/axes/_base.py:3610: MatplotlibDeprecationWarning: The ymaxargument was deprecated in Matplotlib 3.0 and will be removed in 3.2. Usetop instead. alternative='top', obj_type='argument') /Users/maya/anaconda3/envs/py36/lib/python3.6/site-packages/matplotlib/axes/_base.py:3604: MatplotlibDeprecationWarning: The yminargument was deprecated in Matplotlib 3.0 and will be removed in 3.2. Usebottom instead. alternative='bottom', obj_type='argument')

Repeated many times.

I have gone into my copy of burst_plot.py and changed all instances of ymax to top and ymin to bottom. I saved the file and relaunched the notebook, but the warning persists.

Initial Update

Hi 👊

This is my first visit to this fine repo, but it seems you have been working hard to keep all dependencies updated so far.

Once you have closed this issue, I'll create separate pull requests for every update as soon as I find one.

That's it for now!

Happy merging! 🤖

Make L=None default in Data.burst_search

Add link to docs in installation notebook

Improve definition of photon selections

Problem

Currently different photons selections are identified by a string. For example 'DA' indicates all photons, 'D' donor channel photons during donor excitation, 'A' donor channel photons during donor excitation and 'AA' acceptor channel photons during acceptor excitation.

These names are ambiguous without reading the docs or the code and don't allow definition of arbitrary combinations.

Proposal

A simple and general way to handle the definition of photon selection would be defining an object with two keyword arguments:

Ph_sel(Dex='DAem', Aex='DAem')

Possible values for Dex or Aex are 'DAem' (default), 'Dem' and 'Aem'. This way the definition would be explicit and unambiguous but still compact and completely general.

To save some typing and increase legibility a useful shortcut would be Ph_sel('all') meaning all-photons (equivalent but more readable than Ph_sel(Dex='DAem', Aex='DAem')) .

Implementation

The easiest implementation would be using a dictionary or a named tuple. The named tuple has the advantage of defining the keywords (the dictionary would accept any key) and to be immutable. The latter property allow to use the Ph_sel object as dictionary key in cases where the strings 'DA', 'D', 'A' are currently used.

The only disadvantage of the named tuple is that we need to subclass it in order to support the special shortcut ('all') or default keyword values.

A quick search reveals possible solutions for the default values:

>>> from collections import namedtuple
>>> Node = namedtuple('Node', 'val left right')
>>> Node.__new__.__defaults__ = (None, None, None)

or

class Node(namedtuple('Node', ['value', 'left', 'right'])):
    def __new__(cls, value=None, left=None, right=None):
        return super(Node, cls).__new__(cls, value, left, right)

(from Named tuple and optional keyword arguments)

The latter approach would easily allow defining shortcut argument 'all'.

Export "Inter-photon delay times" graph to .csv file

Hi,
I want to export the data of the "Inter-photon delay times" to a .csv file. How would I do it?

Outdated FRET/S histograms with multiple burst searches

Summary

A bug in Data may result in FRET or S histograms plots not to change after performing a second burst search on the same object (but no burst selection).

Many thanks to Danielis Rutkauskas for reporting this bug and providing an example notebook to reproduce it.

Steps to reproduce

In FRETBursts (up to version 0.5.9 included) the following steps reproduce the bug:

d.burst_search(F=6)
dplot(d, hist_fret)

[FRET histogram plot]

d.burst_search(F=12)  # change some parameter
dplot(d, hist_fret)

[FRET histogram plot does not change]

Notes:

• Although the second histogram plot is wrong, the data in d.E and d.S (and all the other burst attributes d.mburst, d.nd, d.na, etc...) is correct.

• This bug does not affect the common workflow where burst selection is performed and the plot is performed on the object containing the bursts selection.

• This bug does not affect fitting of FRET/S histograms when the fitter object is created by the user using bext.burst_fitter() as indicated in the docs.

Workaround

The easiest workaround is to reload the data each time you want to compare FRET or S histograms from different burst searches with no burst selection.

A second workaround that does not require reloading the data is replacing the FRET hist command:

dplot(d, hist_fret)

with:

if hasattr(d, 'E_fitter'): delattr(d, 'E_fitter')
dplot(d, hist_fret)

similarly the S hist command should be:

if hasattr(d, 'S_fitter'): delattr(d, 'S_fitter')
dplot(d, hist_S)

This bug will be fixed in FRETBursts version 0.6 and above.

Internal details

The bug is due to the incorrect handling of E_fitter and S_fitter attributes in Data. The Data object inherits from dict and all its attributes are also dictionary entries. The only exception are the E_fitter and S_fitter attributes that are not part of the dictionary and therefore are not copied when performing burst selection or Data.copy(). The fitter objects compute the FRET and S histograms (as well as computing the fitting and KDE).

When calling hist_fret or hist_S, new E_fitter or S_fitter attributes are created if not present. When the fitter attribute is already present (as in subsequent plots on the same objects) and other histograms parameters are not changes then the old object is used.

When performing a new burst search, the old E_fitter and S_fitter attributes should be deleted but they are not. This results in subsequent calls of the same plot to use the old histogram data.

Background histogram fit with lmfit

Use lmfit to fit the photon interval distribution tail.

• Refactor bg.exp_hist_fit to use lmfit.
• Add a new function for MLE histogram fit like implemented here.

burst_search and burst_search_and_gate

Hi,
I would like to understand what the searches: burst_search and burst_search_and_gate do.

In the following burst_search command:

d.burst_search(F = 6, m = 10, ph_sel = Ph_sel('all'))

We are looking for bursts where 10 (m) consecutive photons are in a time window ≤T (calculated from 6 (F) times the background rate). Do this 10 photons take into account Donor and Acceptor, (so we are not identifying the burst separately for donor and acceptor, but looking at them together)?
Does the burst_search function go through d.ph_times_m, and looks at 10 consecutive timestamps (donor and acceptor)?

This is the result we obtain after using this search:

We get a different result when using burst_search_and_gate:

d_dcbs = bext.burst_search_and_gate(d, F=6,m=10, ph_sel1=Ph_sel(Dex='Dem'), ph_sel2=Ph_sel(Dex='Aem'))

In this case are we identifying the bursts first separately for donor and acceptor, and then only taking into account the bursts present in both, and their duration is the intersection of the two overlapping bursts?

In our case, smFRET, should we use burst_search?
Is burst_search_and_gate used for PIE measurements but comparing the bursts after donor excitation to bursts after acceptor excitation?

Installation of Cython extensions

Add commands to compile Cython extensions in the Installation notebook.

Also add hook to setup.py to build the extensions and copy the binary in the installation dir.

Put installation instructions only in one place

Currently the installation instruction are reproduced both in the README.md and in the documentation.

There should be only one copy to ease maintenance.

Use versioneer to track version

Currently to track the current running FRETBursts version git must be installed.

Using versioneer the version string in each tarball or source tree is automagically generate based on tags and git revision id without the need of having git installed.

This also tracks the full revision (down to the revision id) for FRETBursts installed system-wide.

Add burst_plot module to sphinx docs

Reading Photon-HDF5 files of 'generic' type into FRETbursts

Hi,

I have data from a PIE-MFD setup (i.e. nsALEX with polarization) that I am trying to read from a Photon-HDF5 file into FRETbursts (file is attached). The file specifies the number of polarizations (/setup/num_polarization_ch = 2) and the polarization assignment (/photon_data/measurement_specs/detectors_specs/polarization_ch1-2), yet I get the following error upon trying to read the file:

Invalid_PhotonHDF5: The field /setup/num_polarization_chindicates more than one polarization. However, somedetectors_specs/polarization_ch* fields are missing.

The file was generated through phforge 0.1.1 using phconvert 0.8.2.

Since all the required fields are present when inspecting the file in HDFView, I'm confused why the error occurs.
Thanks for the help!

Best,
Anders

Link to File H20.h5

Save 8-spot data to HDF5

Problem

The 8-spot smFRET binary is poorly designed and space-wise very inefficient.

Would be better to save those files in HDF5 format. This would save space, increase the reading speed and make the data accessible to other programs.

Implementation

Currently 8-spot data files are pickled on first load for efficiency. Instead of pickling, they can be saved with pytables to HDF5 using zlib compression.

Similarly to what is done for the 48-spot format, the class PyTablesList defined in dataload/pytables_array_list.py can be used.

AssertionError when calling calc_bg from .h5 file

After converting my ptu file to hdf5, I am getting an error trying to calculate background rates.

From what I can tell, the hdf5 file is correct as it is loading into HDFView. Here is a link to download the file if that helps: https://www.dropbox.com/s/2n65978ei58kd50/test_DM1_smFRET.h5?dl=0

In [1]: import fretbursts as fb
 - Optimized (cython) burst search loaded.
 - Optimized (cython) photon counting loaded.
--------------------------------------------------------------
 You are running FRETBursts (version 0.5.9).

 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]: d = fb.loader.photon_hdf5('/Users/xxx/05_SM_data/test_DM1_smFRET.h5')

In [3]: d.calc_bg(fb.bg.exp_fit, time_s = 30, tail_min_us='auto')
 - Calculating BG rates ... ---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
<ipython-input-3-66994717892f> in <module>()
----> 1 d.calc_bg(fb.bg.exp_fit, time_s = 30, tail_min_us='auto')

/Users/nuriacodina/anaconda/lib/python3.5/site-packages/fretbursts/burstlib.py in calc_bg(self, fun, time_s, tail_min_us, F_bg, error_metrics)
   1467         if tail_min_us == 'auto':
   1468             bg_auto_th = True
-> 1469             Th_us = self._get_auto_bg_th_arrays(F_bg=F_bg)
   1470         else:
   1471             bg_auto_th = False

/Users/nuriacodina/anaconda/lib/python3.5/site-packages/fretbursts/burstlib.py in _get_auto_bg_th_arrays(self, F_bg, tail_min_us0)
   1364             for ich, ph in enumerate(self.iter_ph_times(ph_sel=ph_sel)):
   1365                 if ph.size > 0:
-> 1366                     bg_rate, _ = bg.exp_fit(ph, tail_min_us=tail_min_us0)
   1367                     th_us[ich] = 1e6 * F_bg / bg_rate
   1368             Th_us[ph_sel] = th_us

/Users/nuriacodina/anaconda/lib/python3.5/site-packages/fretbursts/background.py in exp_fit(ph, tail_min_us, clk_p, error_metrics)
    122     return _exp_fit_generic(ph, fit_fun=exp_fitting.expon_fit,
    123                             tail_min_us=tail_min_us, clk_p=clk_p,
--> 124                             error_metrics=error_metrics)
    125 
    126 def exp_cdf_fit(ph, tail_min_us=None, clk_p=12.5e-9, error_metrics=None):

/Users/nuriacodina/anaconda/lib/python3.5/site-packages/fretbursts/background.py in _exp_fit_generic(ph, fit_fun, tail_min_us, tail_min_p, clk_p, error_metrics)
     85         tail_min = tail_min_us*1e-6/clk_p
     86 
---> 87     res = fit_fun(dph, s_min=tail_min, calc_residuals=error_metrics is not None)
     88     Lambda, residuals, x_residuals, s_size = res
     89 

/Users/nuriacodina/anaconda/lib/python3.5/site-packages/fretbursts/fit/exp_fitting.py in expon_fit(s, s_min, offset, calc_residuals)
     75     """
     76     if s_min > 0: s = s[s >= s_min] - s_min
---> 77     assert s.size > 10
     78 
     79     # Maximum likelihood estimator of the waiting-time

AssertionError: 

In [4]

Document the HDF5-Ph-Data file format draft

Properly reference in the Sphinx docs:

https://github.com/tritemio/FRETBursts/wiki/HDF5-smFRET-format---Draft-0.2-(WIP)

Burst data

I'd like to extract the burst data after a burstsearch and store it in a file. As far as I can see, the burst start and stop information is stored in the burstlib.data object. How can I add the detection channel and excitation period information?

Add loader function for nsALEX

Reading nsALEX file requires only writing a loader function, similar to the usALEX one.

A function to decode the binary SPC Format (Beker & Hickl) is already present in dataload/spcreader.py.

Generate docs for the fitting function

Currently the documentation includes docstring only for the exponential fitting functions when discussing the background fit.

Add documentation for functions in gaussian_fitting.py used for FRET distribution fit.

Probably the docs should go here. Adding references for both exponential and gaussian fittings.

Improve visualisation with Seaborn

Seaborn builds upon Matplotlib and provides several high-quality preset plots.

Integrating Seaborn would improve the data representation and would simplify creation of complex plots.

For example, seaborn has builtin function to plot 2-D histograms with marginal histograms. Or distributions with histograms, rug-plots and kernel density estimation.

The violin plot is also interesting to plot the distribution of fitted FRET efficiencies in each channel:

http://www.stanford.edu/~mwaskom/software/seaborn/examples/violinplots.html

Hdf5 file for Alex

notebook and data.zip

Hi,

I am trying to process some data (data and notebook attached) using ALEX-related functions. However, I am not able to get past the bpl.plot_alternation_hist(). It all works, however, on the example data file 0023uLRpitc_NTP_20dT_0.5GndCl.hdf5. So the question is: what is so specific to the data file that it is recognized as of ALEX?
Also: why doesn't the data loaded from my file contain attribute d.ph_times_t but instead d.ph_times_m which should only be present after applying loader.alex_apply_period()?

Unrelated issue: after updating to 0.6.3, the combination
filename = OpenFileDialog()
d = loader.photon_hdf5(filename)
does not work. What is the reason?

Thank you,

Danielis

Move `data_dir` definition

Currently a data_dir variable is defined in path_def.py.

Setting this variable to the data folder eases defining file names to load for the analysis.
However modifying the variable in path_def.py will be detected as a source code modification by git.

We need to provide a way for the user to set data_dir without modifying the FRETBursts sources.

A possible solution would be defining data_dir in the helper script load_fretbursts. The user will copy the script in his notebook folder (or will find it already there if downloading FRETBursts_notebooks). So modifying it will not touch the FRETBursts sources.

Another solution would be using some type of configuration file in the user folder, similar to what ipython does.

Error on burstsearch when the leakage value is not specified

Without the specification on Lk, Burstsearch is not working and I have got the message below...
(It was running on Python 3.5 version.)

C:\Users\SangYoon\Anaconda3\lib\site-packages\fretbursts\burstlib.py in get_leakage_array(self)
   2275         leakage = self.leakage
   2276         Lk = np.r_[[leakage]*self.nch] if np.size(leakage) == 1 else leakage
-> 2277         Lk *= self.chi_ch
   2278         return Lk
   2279 

TypeError: Cannot cast ufunc multiply output from dtype('float64') to dtype('int32') with casting rule 'same_kind'

Print git commit date on loading

Loading FRETBursts in a notebook with load_fretbursts.py displays a one-liner description of the current HEAD commit showing SHA1 ID and commit message.

Would be useful to add also the commit date.

Direct A excitation correction

Hi,
a question about the definition of the correction factor for direct excitation of acceptor with donor excitation.
Excerpts from the tutorial notebook:
d = loader.photon_hdf5(filename)
d.dir_ex = 0.04

What is the meaning of dir_ex - is it d according to Lee et al. and also dAA according to "Derivation of FRET and S correction formulas", or is it dT according to "Derivation of FRET and S correction formulas" and also dir_ex_t according to "FRET Correction Formulas" in the online manual?

Consider using lmfit for fitting

Currently all the fitting are performed using the standard scipy functions.

In particular, 2-Gaussiam mixture fit of the FRET histogram is currently implemented using scipy.optimize.minimize. The only methods supporting boundaries (constrains) are SLSQP and L-BFGS-B.

Lmfit is a nice pythonic wrapper for non-linear least squares and other methods. Moreover lmfit displays nice fit summaries with several goodness-of-fit statistics.

Using lmfit instead of "home made" functions would make the fitting infrastructure more flexible and easy to use while reducing the maintenance burden.

PySide and PyQt4

Anaconda and WinPython distributions use PyQt4 by default. Anaconda switched to PyQt4 recently but PySide is still available for the foreseeable future.

FRETBursts currently uses PySide, so the user is required to change the default backend from PyQt4 to PySide (editing matplotlibrc) to avoid compatibility issues with non-inline plots and external GUI tools.

Switching to PyQt4 would remove this requirement on the user.

The switch to PyQt4 should be relatively easy. However, PyQt4 offers 2 APIs v1 and v2. v2 is much simpler and the default for python 3. The v2 API can also be used with python 2 but the compatibility with matplotlib and spyder has to be checked first.

Also, seems that IPython has slightly better support for PyQt4.

Bursts.recompute_index_reduce bug

In the rare case when the stop of a burst coincides with the last photon in the timestamps array AND the last two photons are the same timestamps (e.g. donor and acceptor channel registered a photon during the same 10ns time bin) then the method Bursts.recompute_index_reduce raises an error.

The problem is at this line:

• https://github.com/tritemio/FRETBursts/blob/0.5.1/fretbursts/burstsearch/burstsearchlib.py#L563

and happens because it becomes bigger than times_reduced.size.

Kudos to @chungjjang80 for reporting this.

Add test for histogram fitting with complex model

Add a test of fitting distributions with 2 Gaussian peaks + bridge model.

Background fitting

Executing

d.calc_bg(bg.exp_fit, time_s=3600, tail_min_us=4000)
dplot(d, hist_bg, show_fit=True)

produces the attached picture

where the approximation is obviously off the experimental curves. What is the reason for this and how could it be corrected?

Rework Notebook intialization

Currently FRETBursts is usually imported in a notebook through an helper script (load_fretbursts.py) that also performs some notebook-specific initialization and imports.

Any notebook-specific initialization should be done in a ad-hoc function like init_notebook().

Advantages:

• Allows importing FRETBursts from system installation and still add the notebook customizations
• Allows switching between different notebook backends for the figure: matplotlib inline, nbagg, mplotd3, etc...

Add fretmath module to the docs

AND gate burst search speed-up

AND-gate burst search can be speeded up avoiding counting photons, corrections and FRET for the preliminary burst search.

OpenFileDialog in 0.6.3

Hi,

In FRETBursts release 0.6.3 the combination:
filename = OpenFileDialog()
d = loader.photon_hdf5(filename)
does not work, whereas it functioned with older versions. Can this be fixed?

Danielis

Docs: updates Jupyter notebook info

The homepage needs to link more recent Jupyter Notebook docs and the reference beginner section in the reference docs need to be trimmed adding a link to:

http://jupyter-notebook-beginner-guide.readthedocs.org/

hist2d_alex bug with selected bursts

Scenario: hist2d_alex is used on a Data variable d, then a new variable ds is created by selecting bursts in d.

Problem: when using hist2d_alex on ds, the histogram for d is shown instead if the same binwidth of the first plot is used.

Why: the ES_binwidth attribute is not reset on burst selection. ES_binwidth is then used to check whether the ES histogram needs to be recomputed.

Figures

ALEX Jointplot

Error when selecting bursts with nofret=True

After bursts selection with "nofret=True", the program doesn't show a histogram correctly in some cases.
For example,
Case1)

ds=Sel(d, select_bursts.size, th1=50, nofret=True)
ds_FRET=Sel(ds, select_bursts.S, S1=0.3, S2=0.7, nofret=True)
hist_fret(ds)
hist_fret(ds_FRET)

Case2)

ds=Sel(d, select_bursts.size, th1=50, nofret=True)
hist_fret(ds)
ds_FRET=Sel(ds, select_bursts.S, S1=0.3, S2=0.7, nofret=True)
hist_fret(ds_FRET)

The program shows the right histograms in case1 but not in case2.
With "nofret=False", the program always shows the histogram correctly.

Use floating point division

Add:

from __future__ import division

to all the modules and check where the division is expected to return an integer, and use // instead.

Sort fit parameters

Best-fit values of fitted parameters are stored in the DataFrame mfit.MultiFitter.params.

Columns (i.e. parameter's names) should be sorted to avoid spurious changes in the notebook upon re-execution.

correction: Acceptor Direct Excitation

Hi,
I want to ask for advise on how to correct for the acceptor direct excitation. We've observed than when we have a solution of only acceptor dye, exciting with the donor laser, we still detect photons in the donor and acceptor channel:

In the "smFRET" mode (with only one excitation laser), is FRETBursts correcting for direct excitation?

In "nsALEX" (PIE measurement), we've seen that is applied: na -= naa*dir_ex
Is it correct to measure dir_ex, as:
(Acceptor-detector counts during Donor excitation) / (Acceptor-detector counts during Acceptor excitation)
in a sample of Acceptor only?

We've found this paper: Applying Corrections in Single-Molecule FRET
(http://biorxiv.org/content/early/2017/02/01/083287)
Could you explain me how to measure Dir from Definition 2?

Update installation docs

Add a note in the docs about running the notebook "Installation.ipynb" that saves a configuration file containing FRETBursts source path.

This allows load_fretbursts.py to find and import the software without relying on manual path entry from the user.

Refactor hist2d_alex

The plot function hist2d_alex is an old complex function with a visualization (2-D histogram + scatterplot) now superseded by the more powerful alex_jointplot.

hist2d_alex should be refactored (and maybe renamed to hist_ES or distribution_ES?) to produce a single E-S "distribution" plot (hexplot, kde, scatterplot) without marginals. Next, this can be used to simplify the implementation of alex_jointplot and to provide a more uniform API for both alex_jointplot and hist_ES.

Add loader.py to the docs

Add the loader functions to the sphinx documentation.