OpenSMFS / FRETBursts

Burst analysis software for single and multi-spot single-molecule FRET (smFRET) data.
https://opensmfs.github.io/FRETBursts/
GNU General Public License v2.0
12 stars 9 forks source link

Weird Background correction #27

Closed BenjaminAmbrose closed 5 years ago

BenjaminAmbrose commented 5 years ago

Hi, I'm having some strange issues and I was wondering if you could provide insight into what's going on

I've recently started seeing a strange high FRET population (>1) in all my experiments, but for some reason these are only appearing in FRETBursts, our acquisition software's live view doesn't register them, and when I analyse the data in PAM it doesn't find them either.

I isolated a few of them on a time trace (pic of one attached) to see what was going on and it looks like there's plenty of DD photons there, but when queried FRETBursts is telling me the burst has negative counts in that channel. The background FRETBursts is calculating doesn't seem nearly high enough to account for this either.

If you have any insight into what might be going on here it would be much appreciated as it's got us scratching our heads. I can send hdf5's if that helps.

ES plot bg burst trace

tritemio commented 5 years ago

I agree that by looking at the timetrace, that burst should have FRET < 1.

You can look at the output of:

bext.burst_data(d, include_bg=True)

this will give you a nicely formatted table with both corrected and raw burst counts as well as the background component. I would look there too see if the numbers make sense and where the negative nd comes from.

Note that bext.burst_data(d) calls d.burst_data_ich(0) (code). The actual background component for each burst is computed in d.expand() (code). The correction of raw counts happens in d.background_correction() (code).

You can compare the FRET histograms before and after background correction to see the difference. Also, I would align the FRET histogram bins so that 0 and 1 fall in the middle of a bin in order to reduce visual artifacts. For example:

hist_fret(d, bins=np.arange(-0.2, 1.2, 0.025) + 0.025/2)

Hope this helps.

BenjaminAmbrose commented 5 years ago

Thanks for the help! I'm still not sure how to plot without background correction, however I've looked at the tables with bext.burst_data like you said and found something odd

The regular bursts with 0-1 FRET look like this:

size_raw t_start t_stop bg_period width_ms bg_ad bg_dd bg_aa bg_da E S nd na nt nda naa spot
156 2.030744 2.036433 0 5.68936 2.141787 3.532776 1.853006 0.203481 0.814376 0.561218 15.467224 67.858213 148.472431 -0.203481 65.146994 0
353 3.374009 3.381278 0 7.26967 2.736702 4.514060 2.367708 0.260001 0.804404 0.440296 29.485940 121.263298 342.381529 0.739999 191.632292 0
94 8.560643 8.562619 0 1.97606 0.743897 1.227023 0.643596 0.070674 0.844631 0.547451 7.772977 42.256103 91.385483 -0.070674 41.356404 0

But the >1 FRET ones look like this:

size_raw t_start t_stop bg_period width_ms bg_ad bg_dd bg_aa bg_da E S nd na nt nda naa spot
130 55.149405 55.153257 0 3.85167 1.449980 2.391673 1.254477 0.137755 1.033611 0.569705 -2.391673 73.550020 124.903871 -0.137755 53.745523 0
169 59.451233 59.457389 0 6.15636 2.317592 3.822757 2.005107 0.220183 1.073714 0.322401 -3.822757 55.682408 160.854545 -0.220183 108.994893 0
115 124.287127 124.289594 2 2.46788 0.960638 1.610612 0.892982 0.078901 1.022869 0.631445 -1.610612 72.039362 111.535769 -0.078901 41.107018 0

I've put timetraces of the bursts below. If you look at the table it seems that it's just subtracting the background from 0, as if it's not registering any DD counts at all. Maybe this isn't a background issue after all, and FRETBursts is just randomly ignoring the DD photons in some bursts? I still can't make sense of it.

mysteryhighfretbursts

Thanks again, your help is much appreciated

BenjaminAmbrose commented 5 years ago

Okay I think I may have found the problem, a while back we updated our acquisition software and in doing so messed around a bit with the way it saves files. I've gone into the raw files and looked specifically at the regions where FRETBursts is messing up and it looks like our software is occasionally printing the timestamps non chronologically (actually in blocks of each detector, explaining why it's counting A photons but no D photons).

This would make sense as to why these photons are appearing on the timetrace but not counted in bursts, as they are still there but they aren't actually positioned between the start and stop photons that define the burst.

This is our fault and I'll get to fixing it, but it might be useful to you anyway to know FRETBursts does this when the timestamps aren't properly sequential.

Thanks anyway!

tritemio commented 5 years ago

I'm glad you found the issue, yes monotonicity of timestamps is assumed pretty much everywhere in FRETBursts.

On a side note, if you convert the files to Photon-HDF5, all these checks of monotonicity are automatically done before saving.

As for the FRET histogram with no background, you are right, I think there is no pre-made function in FRETBursts, but you can use plain matplotlib:

df = bext.burst_data(d, include_bg=True)
nd_raw = df.nd + df.bg_dd
na_raw = df.na + df.bg_ad
E_raw = na_raw / (nd_raw + na_raw)
plt.hist(E_raw, bins=bins=np.arange(-0.2, 1.2, 0.025) + 0.025/2))

(I didn't test it, there may be some typo)

BenjaminAmbrose commented 5 years ago

Ah that's useful, do the checks happen in the notebooks you've provided or are they part of phconvert itself? Because I have been using:

phconvert.hdf5.save_photon_hdf5(data, h5_fname=savename, overwrite=True)

and this has not thrown up any flags

tritemio commented 5 years ago

Checking it ... the test looks like this:

assert (np.diff(timestamps) >= 0).all()

and it is only done for PicoQuant files. There is one check in the notebook and one in the unit tests, but no realtime check using phc.save_photon_hdf5, I was recalling it wrong.

I think such test could be added to phc.save_photon_hdf5, but I would prefer having it optional for performance reasons. Something like extensive_tests=True could be a catch-all flag for additional (potentially slow) tests.