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

Polarization Data (again) #25

Closed ChristianGebhardt closed 5 years ago

ChristianGebhardt commented 5 years ago

Hello FRETBursts-developer, I've seen that there where quite some issues on that topic, but I am still struggling with getting your example notebook "FRETBursts - ns-ALEX example.ipynb" to run with polarization data. I have the file that I am using here: https://www.dropbox.com/s/qf9oxu3ybhno9rz/01_MalE_Alexa555_647_1.h5?dl=0 First problem is that the command d = loader.photon_hdf5(filename) doesn't detect that these are polarization data. This might be, because the /photon_data/measurement_specs/measurement_type is set to "smFRET-nsALEX". The flag "-2pol" is only added, if I use "generic" (loader.py, line 122/123). Shall I change my measurement_type to "generic"? If I add the flag "-2pol" to the measurement type by putting the if statement if setup.num_polarization_ch.read() > 1 out of the outer if statement, I can import the data properly with the polarization information stored in d.det_p_s_pol (array([0, 2], dtype=uint8), array([1, 3], dtype=uint8)).

However, when I now run the command bpl.plot_alternation_hist(d) (or also the next command loader.alex_apply_period(d)), I am facing the problem that the code cannot deal with the 2 polarization channels: In burst_plot.py, line 248/249, there is the selection nanotimes_d = d.nanotimes_t[ich][d.det_t[ich] == d_ch], but d_ch is now an array. The same is with loader.alex_apply_period(d) in loader.py, line 633/634, and also later.

Do I have to change something with my file-format? Or is this code not working with the polarization information? I fear that if I change the code to selecting both polarizations that I have to do this quite often.

Note: If I select just to channels by d.add(det_donor_accept=(0,2)), everything works fine. But of course, I lose half of my photons. Is there a simple way to merge the polarizations as long as I do not use the polarization information for the evaluation?

One further thing: the function _det_p_s_pol_multich (burstlib.py, line 1025) also raised an Attribute error with 'det_s_p_pol'. In other files it is called 'det_p_s_pol'. Maybe p and s are mixed up a bit here? Or should there be two different versions?

Thanks for your support, Christian

tritemio commented 5 years ago

Hi @ChristianGebhardt,

yes you have to change the measurement type to "generic" in order to use polarization. In Photon-HDF5, smFRET-nsALEX is a specific measurement type with 2 colors and no polarization. Initially we though about adding measurement types for each type of data, but later we figured that the number of combinations would have been too high. So, in Photon-HDF5 v 0.5, we introduced a measurement type called generic, where all the setup information is stored in the various fields in /setup. See the docs here:

https://photon-hdf5.readthedocs.io/en/0.5.dev/generic.html

FRETBursts polarization support is preliminary and only available on master. Make sure you are using the latest FRETBurst from master which has some fixes to load polarization data. The old issue is #18, and the PR with the fix was #19.

Regarding the command d.add(det_donor_accept=(0,2)) it should not be needed.

Finally, regarding the det_s_p_pol AttributeError, good catch, this is a typo and needs to be fixed, thanks for the report!

Please fix the data file using generic first, so we can see if there are other bugs.

ChristianGebhardt commented 5 years ago

Hi @tritemio ,

thanks for your support. I am getting closer to run everything.

I use the master branch now and I have changed the type to generic and loading the data works fine. The sample file can be found here, if needed:

https://www.dropbox.com/s/lrf90f0ux74zs8h/01_MalE_Alexa555_647_generic.h5?dl=0

I still have the problem: In burst_plot.py, line 248/249, there is the selection nanotimes_d = d.nanotimes_t[ich][d.det_t[ich] == d_ch]. But with d_ch as an array [0,2], the selection doesn't work. I could solve that, by using the same strategy like with the selection masks in loader.py and replacing the code by

values = np.atleast_1d(d_ch)
nanotimes_d = d.nanotimes_t[ich][d.det_t[ich] == values[0]]
for v in values[1:]:
   nanotimes_d = np.concatenate([nanotimes_d,d.nanotimes_t[ich][d.det_t[ich] == v]])

It is probably nicer to outsource it to a function, but it works.

With the command loader.alex_apply_period(d), I have other problems:

  1. The selection mask, which is created by _selection_mask(arr, values) (loader.py, line 220), was False for all photons. I could solve that by changing line 226 from mask *= arr == v to mask += arr == v. So, it's true if the photon is in one of the two polarization channels, but not in both.
  2. I get an assert error in _get_det_masks(det_t, det_ch1, det_ch2, valid, mask_ref=None, ich=0) (loader.py, line 689), when the channels do not cover the full time range. With my settings d.add(A_ON=(200, 1500), D_ON=(1750, 3200)), all photons out of that range throw the error. My workaround was to set it to d.add(A_ON=(-1, 1501), D_ON=(1500, 5000)). But if I really want to limit the range, this is not ideal.

One last thing, I found another mixing up of det_s_p_pol with det_p_s_pol in loader.py, line 304

After I changed those things, I was able to do all the following steps. Only the assert error wasn't solved completely.

I hope, that those inputs can help.

Regards, Christian

tritemio commented 5 years ago

Hey @ChristianGebhardt, good digging 👍. I made a PR (#26) with your proposed fixes.

Point by point:

I plan to merge the PR in a few days, let me know if you can test it.

In the future, feel free to open a PR with the proposed fixes ;-).

ChristianGebhardt commented 5 years ago

Thanks for the quick response!

I checked the new branch fretburst-pol-fix and everything works fine now. I couldn't reproduce the asset error anymore. Maybe, I changed to much in the code during debugging.

I think that we can close this issue.

tritemio commented 5 years ago

Great, thanks for reporting back!