fname = 'data/pax-2017-05-23_04_22d.hdf5'
from pathlib import Path
fname = Path(fname)
assert fname.is_file(), 'File not found.'
mlabel = '_'.join(fname.stem.replace('pax-', '').replace('alex-', '').split('_')[:4])
mlabel
'2017-05-23_04_22d'
!date
Thu Oct 5 16:13:54 PDT 2017
import os
from pathlib import Path
import numpy as np
from IPython.display import display, HTML, Math
import pandas as pd
import matplotlib as mpl
mpl.rcParams['font.sans-serif'].insert(0, 'Arial')
%matplotlib inline
import matplotlib.pyplot as plt
from heatmap import heatmap48, spotsh, spotsv
import pybroom as br
/home/anto/miniconda3/envs/py36/lib/python3.6/site-packages/seaborn/apionly.py:6: UserWarning: As seaborn no longer sets a default style on import, the seaborn.apionly module is deprecated. It will be removed in a future version. warnings.warn(msg, UserWarning)
from fretbursts import *
- Optimized (cython) burst search loaded. - Optimized (cython) photon counting loaded. -------------------------------------------------------------- You are running FRETBursts (version 0.6.5). 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 --------------------------------------------------------------
sns = init_notebook(apionly=True)
plt.rcParams['font.size'] = 14
%config InlineBackend.figure_format = 'retina'
# Instead of importing, "load" some utility functions using `%run -i`.
# This gives these functions access to variables defined in the notebok (e.g. mlabel).
# Defines: savefig, save_name, info_html, cal_phrate, make_df_spots, make_df_bursts
# cal_phrate_alex, make_df_bursts_alex
%run -i utils.py
plot_timetraces = True
skip_ch = (12, 13)
save_figures = True
savefigdir = 'figures'
highres = True
d = loader.photon_hdf5(str(fname), ondisk=True)
info_html(d)
22d dsDNA , 20pM.
fig, ax = plt.subplots(figsize=(12, 8))
bpl.plot_alternation_hist_usalex(d, ax=ax, bins=np.arange(0, 4097, 16))
%%timeit -n1 -r1
loader.alex_apply_period(d)
# Total photons (after ALEX selection): 345,374,281 # D photons in D+A excitation periods: 170,470,305 # A photons in D+A excitation periods: 174,903,976 # D+A photons in D excitation period: 129,171,591 # D+A photons in A excitation period: 216,202,690 1min 14s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)
if plot_timetraces:
num_time_points = 5
kws = dict(figsize=(24, 8), xrotation=45)
# Timepoints equally distributed along the measurement
time_points = np.round(np.linspace(d.time_min+1, d.time_max-2, num=num_time_points))
for i in time_points:
dplot(d, timetrace, tmin=i, tmax=i+1, **kws);
plt.ylim(-100, 100)
savefig("%s_timetrace_t=%d-%d.png" % (mlabel, i, i+1))
Saved: figures/2017-05-23_04_22d_2017-05-23_04_22d_timetrace_t=1-2.png.png Saved hires: figures/2017-05-23_04_22d_2017-05-23_04_22d_timetrace_t=1-2_highres.png Saved: figures/2017-05-23_04_22d_2017-05-23_04_22d_timetrace_t=359-360.png.png Saved hires: figures/2017-05-23_04_22d_2017-05-23_04_22d_timetrace_t=359-360_highres.png Saved: figures/2017-05-23_04_22d_2017-05-23_04_22d_timetrace_t=717-718.png.png Saved hires: figures/2017-05-23_04_22d_2017-05-23_04_22d_timetrace_t=717-718_highres.png Saved: figures/2017-05-23_04_22d_2017-05-23_04_22d_timetrace_t=1074-1075.png.png Saved hires: figures/2017-05-23_04_22d_2017-05-23_04_22d_timetrace_t=1074-1075_highres.png Saved: figures/2017-05-23_04_22d_2017-05-23_04_22d_timetrace_t=1432-1433.png.png Saved hires: figures/2017-05-23_04_22d_2017-05-23_04_22d_timetrace_t=1432-1433_highres.png
d.calc_bg_cache(bg.exp_fit, time_s=5, tail_min_us='auto', F_bg=1.7)
* Loading BG rates from cache ... - Loading bakground data: [DONE]
ax = dplot(d, timetrace_bg, show_da=True, hspace=0, wspace=0.08,
plot_style=dict(marker=None), title='right top',
title_kws=dict(fontsize=16), xrotation=45);
plt.xlim(0)
plt.ylim(0, 8);
bg_AexDem = d.bg_from(Ph_sel(Aex='Dem'))
bg_AexAem = d.bg_from(Ph_sel(Aex='Aem'))
bg_DexDem = d.bg_from(Ph_sel(Dex='Dem'))
bg_DexAem = d.bg_from(Ph_sel(Dex='Aem'))
fig, ax = plt.subplots(1, 2, figsize=(12, 4), sharex=True, sharey=True)
fig.suptitle("BG A-excitation period")
t = np.arange(len(bg_AexDem[0])) * 5
ax[0].plot(t, np.array(bg_AexDem).T, color='g', alpha=0.2);
ax[1].plot(t, np.array(bg_AexAem).T, color='r', alpha=0.2);
plt.setp(ax, xlabel='Time (s)', ylabel='cps')
plt.subplots_adjust(wspace=0.07)
ax[0].text(0.05, 0.95, '$A_{EX}D_{EM}$', va='top', fontsize=16, transform=ax[0].transAxes)
ax[1].text(0.95, 0.95, '$A_{EX}A_{EM}$', va='top', ha='right', fontsize=16, transform=ax[1].transAxes);
fig, ax = plt.subplots(1, 2, figsize=(12, 4), sharex=True, sharey=True)
fig.suptitle("BG D-excitation period")
t = np.arange(len(bg_DexDem[0])) * 5
ax[0].plot(t, np.array(bg_DexDem).T, color='g', alpha=0.2);
ax[1].plot(t, np.array(bg_DexAem).T, color='r', alpha=0.2);
plt.setp(ax, xlabel='Time (s)', ylabel='cps')
plt.subplots_adjust(wspace=0.07)
ax[0].text(0.05, 0.95, '$D_{EX}D_{EM}$', va='top', fontsize=16, transform=ax[0].transAxes)
ax[1].text(0.95, 0.95, '$D_{EX}A_{EM}$', va='top', ha='right', fontsize=16, transform=ax[1].transAxes);
d.burst_search(min_rate_cps=50e3, pax=True)
- Performing burst search (verbose=False) ...[DONE] - Calculating burst periods ...[DONE] - Counting D and A ph and calculating FRET ... - Applying background correction. [DONE Counting D/A]
kws = dict(skip_ch=skip_ch, hspace=0, wspace=0, top=0.96,
title_bg=False, title_nbursts=False, title='in')
dplot(d, hist_size_all, **kws);
plt.xlim(0, 120)
plt.legend();
kws.update(title_nbursts=False)
ax = dplot(d, hist_size, which='na', **kws);
dplot(d, hist_size, which='naa', AX=ax, **kws);
plt.xlim(0, 100);
plt.legend(['$n_a$', '$n_{DA_{ex}A_{em}}$'], fontsize=18, loc='upper right')
<matplotlib.legend.Legend at 0x7f7eab37c668>
size_th = 60
size_sel_kws = dict(ph_sel=Ph_sel(Dex='DAem', Aex='DAem'),
na_comp=False, naa_comp=False, naa_aexonly=False)
dplot(d, hist_size, vline=size_th, **size_sel_kws, **kws);
plt.xlim(-10, 300)
plt.legend(fontsize=15, loc='upper right');
Math(d._burst_sizes_pax_formula(**size_sel_kws))
bs_mean = []
for i in range(48):
sizes_ch = d.burst_sizes_pax_ich(i, **size_sel_kws)
tail = sizes_ch > size_th
bs_mean.append(sizes_ch[tail].mean() - size_th)
utils.py:5: RuntimeWarning: Mean of empty slice. from pathlib import Path /home/anto/miniconda3/envs/py36/lib/python3.6/site-packages/numpy/core/_methods.py:80: RuntimeWarning: invalid value encountered in double_scalars ret = ret.dtype.type(ret / rcount)
heatmap48(bs_mean, cmap='viridis', skip_ch=skip_ch,
title="Mean Burst Sizes", vmin=10)
savefig('heatmap_burst_sizes')
Saved: figures/2017-05-23_04_22d_heatmap_burst_sizes.png Saved hires: figures/2017-05-23_04_22d_heatmap_burst_sizes_highres.png
width_th = 0.4
dplot(d, hist_width, vline=width_th, **kws);
plt.xlim(-0.5, 4)
(-0.5, 4)
bw_mean = []
for i in range(48):
widths_ch = d.mburst[i].width * d.clk_p * 1e3
tail = widths_ch > width_th
bw_mean.append(widths_ch[tail].mean() - width_th)
heatmap48(bw_mean, cmap='viridis', skip_ch=skip_ch,
title="Mean Burst Width (ms)", vmin=0.1)
savefig('heatmap_burst_widths')
Saved: figures/2017-05-23_04_22d_heatmap_burst_widths.png Saved hires: figures/2017-05-23_04_22d_heatmap_burst_widths_highres.png
# Use recompute=True if changing burst-search parameters
recompute = False
phrates = {}
streams = ('all', 'DexDem', 'AexDem', 'DexAem', 'AexAem', )
for stream in streams:
print(' - Computing peak photon rates for %6s stream.' % str(Ph_sel.from_str(stream)))
cal_phrate(d, stream=Ph_sel.from_str(stream), phrates=phrates, recompute=recompute)
- Computing peak photon rates for all stream. Valid fraction (mean of all ch): 100.0 % - Computing peak photon rates for DexDem stream. Valid fraction (mean of all ch): 29.8 % - Computing peak photon rates for AexDem stream. Valid fraction (mean of all ch): 30.1 % - Computing peak photon rates for DexAem stream. Valid fraction (mean of all ch): 1.2 % - Computing peak photon rates for AexAem stream. Valid fraction (mean of all ch): 10.1 %
kws = dict(figsize=(18, 6), skip_ch=skip_ch, wspace=0, hspace=0, grid=True,
title=None, title_nbursts=False, title_bg=False, xrotation=45)
ax = dplot(phrates['AexAemB'], hist_burst_phrate,
plot_style=dict(color=bpl.purple, ms=3), **kws);
dplot(phrates['AexDemB'], hist_burst_phrate,
plot_style=dict(color='C0', ms=3), AX=ax, **kws);
dplot(phrates['DexDemB'], hist_burst_phrate,
plot_style=dict(color=bpl.green, ms=3), AX=ax, **kws);
kws.update(vline=60, title='in', title_nbursts=False, title_bg=False,)
dplot(phrates['DexAemB'], hist_burst_phrate,
plot_style=dict(color=bpl.red, ms=3), AX=ax, **kws);
plt.xlim(-50, 1150);
plt.setp(ax[:, 0], ylabel='PDF')
plt.setp(ax[-1], xlabel='kcps')
savefig('peak_phrate_all')
Saved: figures/2017-05-23_04_22d_peak_phrate_all.png Saved hires: figures/2017-05-23_04_22d_peak_phrate_all_highres.png
phr_th = 60e3
for stream in streams:
phrates[stream]['mean'] = 0
phr = phrates[stream+'B']
for ich in range(48):
valid = ~pd.isnull(phr[ich])
phr_valid = phr[ich][valid]
m = phr_valid[phr_valid >= phr_th].mean() - phr_th
phrates[stream].loc[ich, 'mean'] = np.round(m*1e-3, 1)
streams
('all', 'DexDem', 'AexDem', 'DexAem', 'AexAem')
cmaps = ('bone', 'Greens_r', 'Blues_r', 'Reds_r', 'Purples_r')
label = {'DexDem': '$D_{ex}D_{em}$', 'DexAem': '$D_{ex}A_{em}$',
'AexDem': '$DA_{ex}D_{em}$', 'AexAem': '$DA_{ex}A_{em}$', 'all': 'all-photons'}
with plt.rc_context({'font.size': 15}):
for stream, cmap in zip(streams, cmaps):
heatmap48(phrates[stream]['mean'], cmap=cmap, skip_ch=skip_ch,
title=r"Peak Photon Rate in %s (kcps)" % label[stream])
savefig('peak_phrate_heatmap_%s' % stream)
Saved: figures/2017-05-23_04_22d_peak_phrate_heatmap_all.png Saved hires: figures/2017-05-23_04_22d_peak_phrate_heatmap_all_highres.png Saved: figures/2017-05-23_04_22d_peak_phrate_heatmap_DexDem.png Saved hires: figures/2017-05-23_04_22d_peak_phrate_heatmap_DexDem_highres.png Saved: figures/2017-05-23_04_22d_peak_phrate_heatmap_AexDem.png Saved hires: figures/2017-05-23_04_22d_peak_phrate_heatmap_AexDem_highres.png Saved: figures/2017-05-23_04_22d_peak_phrate_heatmap_DexAem.png Saved hires: figures/2017-05-23_04_22d_peak_phrate_heatmap_DexAem_highres.png Saved: figures/2017-05-23_04_22d_peak_phrate_heatmap_AexAem.png Saved hires: figures/2017-05-23_04_22d_peak_phrate_heatmap_AexAem_highres.png
size_th = 80
size_sel_kws = dict(ph_sel=Ph_sel(Dex='DAem', Aex='Dem'), gamma=0.5,
na_comp=True, naa_comp=False, naa_aexonly=False)
ds = d.select_bursts(select_bursts.size, th1=size_th, **size_sel_kws)
Su = [(nd + na)/(nd + na + naa) for nd, na, naa in zip(ds.nd, ds.na, ds.naa)]
ds.add(Su=Su)
Su2 = [(nd + nda + na)/(nd + nda + na + naa)
for nd, na, nda, naa in zip(ds.nd, ds.na, ds.nda, ds.naa)]
ds.add(Su2=Su2)
Math(d._burst_sizes_pax_formula(**size_sel_kws))
heatmap48(ds.num_bursts)
ax = dplot(ds, hist_fret, pdf=False, skip_ch=skip_ch);
plt.setp(ax[:, 0], ylabel='# Bursts');
# plt.ylim(0, 600);
savefig('48spot hist E all-bursts')
/home/anto/src/FRETBursts/fretbursts/mfit.py:329: RuntimeWarning: invalid value encountered in true_divide self.hist_pdf /= self.hist_counts.sum(1)[:, np.newaxis]
Saved: figures/2017-05-23_04_22d_48spot_hist_E_all-bursts.png Saved hires: figures/2017-05-23_04_22d_48spot_hist_E_all-bursts_highres.png
dplot(ds, hist_burst_data, data_name='Su', pdf=False, skip_ch=skip_ch);
plt.setp(ax[:, 0], ylabel='# Bursts');
savefig('48spot hist Su')
/home/anto/src/FRETBursts/fretbursts/mfit.py:329: RuntimeWarning: invalid value encountered in true_divide self.hist_pdf /= self.hist_counts.sum(1)[:, np.newaxis]
Saved: figures/2017-05-23_04_22d_48spot_hist_Su.png Saved hires: figures/2017-05-23_04_22d_48spot_hist_Su_highres.png
ax = dplot(ds, hexbin_alex, vmax=10, skip_ch=skip_ch, wspace=0, hspace=0, gridsize=60,
title='bottom left', title_bg=False);
ax0 = ax[0,0]
fig = ax0.figure
cax = fig.add_axes([0.97, 0.25, 0.01, 0.5])
plt.colorbar(cax=cax)
plt.setp(ax[:, 0], ylabel='S');
ax0.set_xticks([0, 0.5, 1])
ax0.set_yticks([0, 0.5, 1]);
ax0.set_xlim(-0.2, 1)
ax0.set_ylim(0, 1.2);
savefig('48spot alex hist S all-bursts')
Saved: figures/2017-05-23_04_22d_48spot_alex_hist_S_all-bursts.png Saved hires: figures/2017-05-23_04_22d_48spot_alex_hist_S_all-bursts_highres.png
ax = dplot(ds, hexbin_alex, S_name='Su2', vmax=15, skip_ch=skip_ch,
wspace=0, hspace=0, gridsize=60, title='bottom left', title_bg=False);
ax0 = ax[0, 0]
fig = ax0.figure
cax = fig.add_axes([0.97, 0.25, 0.01, 0.5])
plt.setp(ax[:, 0], ylabel='$S_{up}$');
plt.setp(ax[-1], xlabel='$E_{PR}$');
plt.colorbar(cax=cax)
ax0.set_xticks([0, 0.5, 1])
ax0.set_yticks([0, 0.5, 1]);
ax0.set_xlim(-0.2, 1)
ax0.set_ylim(0, 1.2);
ax = dplot(ds, hexbin_alex, S_name='Su', vmax=15, skip_ch=skip_ch,
wspace=0, hspace=0, gridsize=60, title='bottom left', title_bg=False);
ax0 = ax[0, 0]
fig = ax0.figure
cax = fig.add_axes([0.97, 0.25, 0.01, 0.5])
plt.setp(ax[:, 0], ylabel='$S_u$');
plt.setp(ax[-1], xlabel='$E_{PR}$');
plt.colorbar(cax=cax)
ax0.set_xticks([0, 0.5, 1])
ax0.set_yticks([0, 0.5, 1]);
ax0.set_xlim(-0.2, 1)
ax0.set_ylim(0, 1.2);
plt.suptitle('')
savefig('48spot alex hist Su all-bursts')
Saved: figures/2017-05-23_04_22d_48spot_alex_hist_Su_all-bursts.png Saved hires: figures/2017-05-23_04_22d_48spot_alex_hist_Su_all-bursts_highres.png
dplot(ds, hist_fret, pdf=False, dec=4, suptitle=False, title='in', figsize=(10, 3.5), wspace=0);
dplot(ds, hist_fret, pdf=False, dec=2, suptitle=False, title='in', wspace=0, hspace=0);
size_Dex = 80
size_sel_kws = dict(ph_sel=Ph_sel(Dex='DAem', Aex='Dem'), gamma=0.5,
na_comp=True, naa_comp=False, naa_aexonly=False)
Math(d._burst_sizes_pax_formula(**size_sel_kws))
size_Aex = 20
naa_sel_kws = dict(ph_sel=Ph_sel(Aex='Aem'), gamma=0.5,
na_comp=False, naa_comp=False, naa_aexonly=False)
Math(d._burst_sizes_pax_formula(**naa_sel_kws))
ds1 = d.select_bursts(select_bursts.size, th1=size_Dex, **size_sel_kws)
Su = [(nd + na)/(nd + na + naa) for nd, na, naa in zip(ds1.nd, ds1.na, ds1.naa)]
ds1.add(Su=Su)
dss = ds1.select_bursts(select_bursts.size, th1=size_Aex, **naa_sel_kws)
Su = [(nd + na)/(nd + na + naa) for nd, na, naa in zip(dss.nd, dss.na, dss.naa)]
dss.add(Su=Su)
g = alex_jointplot(ds1, i=30, vmax=22, S_name='Su')
g.ax_joint.set_ylabel('$S_u$')
g = alex_jointplot(dss, i=30, vmax=22, S_name='Su')
g.ax_joint.set_ylabel('$S_u$')
/home/anto/src/FRETBursts/fretbursts/mfit.py:329: RuntimeWarning: invalid value encountered in true_divide self.hist_pdf /= self.hist_counts.sum(1)[:, np.newaxis]
<matplotlib.text.Text at 0x7f7ea0bbca20>
ax = dplot(dss, hexbin_alex, vmax=15, S_name='Su', skip_ch=skip_ch,
wspace=0, hspace=0, gridsize=60, title='bottom left', title_bg=False);
ax0 = ax[0,0]
fig = ax0.figure
cax = fig.add_axes([0.97, 0.25, 0.01, 0.5])
plt.colorbar(cax=cax)
plt.setp(ax[:, 0], ylabel='$S_u$');
plt.setp(ax[-1], xlabel='$E_{PR}$');
ax0.set_xticks([0, 0.5, 1])
ax0.set_yticks([0, 0.5, 1]);
ax0.set_xlim(-0.2, 1)
ax0.set_ylim(0, 1.2);
savefig('48spot alex hist Su naa AND size selection')
Saved: figures/2017-05-23_04_22d_48spot_alex_hist_Su_naa_AND_size_selection.png Saved hires: figures/2017-05-23_04_22d_48spot_alex_hist_Su_naa_AND_size_selection_highres.png
ds_fret = dss
E_fitter = bext.bursts_fitter(ds_fret, 'E', skip_ch=skip_ch)
model = mfit.factory_gaussian()
model.set_param_hint('center', value=0.5, min=0, max=1)
model.set_param_hint('sigma', value=0.1, min=0.03, max=0.2)
model.print_param_hints()
Name Value Min Max Vary Expr amplitude 1 -inf inf True center 0.5 0 1 True fwhm nan -inf inf True 2.3548200*sigma height nan -inf inf True 0.3989423*amplitude/max(1.e-15, sigma) sigma 0.1 0.03 0.2 True
E_fitter.fit_histogram(model, pdf=False, method='mealder')
/home/anto/src/FRETBursts/fretbursts/mfit.py:329: RuntimeWarning: invalid value encountered in true_divide self.hist_pdf /= self.hist_counts.sum(1)[:, np.newaxis]
ax = dplot(ds_fret, hist_fret, pdf=False, show_model=True, skip_ch=skip_ch);
plt.setp(ax[:, 0], ylabel='# Bursts')
savefig('FRET pop - E hist fit')
Saved: figures/2017-05-23_04_22d_FRET_pop_-_E_hist_fit.png Saved hires: figures/2017-05-23_04_22d_FRET_pop_-_E_hist_fit_highres.png
S_fitter = bext.bursts_fitter(ds_fret, 'S', skip_ch=skip_ch)
model = mfit.factory_gaussian()
model.set_param_hint('center', value=0.5, min=0, max=1)
model.set_param_hint('sigma', value=0.1, min=0.03, max=0.2)
model.print_param_hints()
Name Value Min Max Vary Expr amplitude 1 -inf inf True center 0.5 0 1 True fwhm nan -inf inf True 2.3548200*sigma height nan -inf inf True 0.3989423*amplitude/max(1.e-15, sigma) sigma 0.1 0.03 0.2 True
S_fitter.fit_histogram(model, pdf=False, method='mealder')
/home/anto/src/FRETBursts/fretbursts/mfit.py:329: RuntimeWarning: invalid value encountered in true_divide self.hist_pdf /= self.hist_counts.sum(1)[:, np.newaxis]
ax = dplot(ds_fret, hist_S, pdf=False, show_model=True, skip_ch=skip_ch);
plt.setp(ax[:, 0], ylabel='# Bursts')
savefig('FRET pop - E hist fit');
Saved: figures/2017-05-23_04_22d_FRET_pop_-_E_hist_fit.png Saved hires: figures/2017-05-23_04_22d_FRET_pop_-_E_hist_fit_highres.png
Su_fitter = bext.bursts_fitter(ds_fret, 'Su', skip_ch=skip_ch)
model = mfit.factory_gaussian()
model.set_param_hint('center', value=0.5, min=0, max=1)
model.set_param_hint('sigma', value=0.1, min=0.03, max=0.2)
model.print_param_hints()
Name Value Min Max Vary Expr amplitude 1 -inf inf True center 0.5 0 1 True fwhm nan -inf inf True 2.3548200*sigma height nan -inf inf True 0.3989423*amplitude/max(1.e-15, sigma) sigma 0.1 0.03 0.2 True
Su_fitter.fit_histogram(model, pdf=False, method='mealder')
/home/anto/src/FRETBursts/fretbursts/mfit.py:329: RuntimeWarning: invalid value encountered in true_divide self.hist_pdf /= self.hist_counts.sum(1)[:, np.newaxis]
ax = dplot(ds_fret, hist_burst_data, data_name='Su', pdf=False, show_model=True, skip_ch=skip_ch);
plt.setp(ax[:, 0], ylabel='# Bursts')
savefig('FRET pop - Su hist fit');
Saved: figures/2017-05-23_04_22d_FRET_pop_-_Su_hist_fit.png Saved hires: figures/2017-05-23_04_22d_FRET_pop_-_Su_hist_fit_highres.png
def make_ES_fit_dataframe(E_df_fitt, S_df_fitt, Su_df_fitt, dx, skip_ch=skip_ch):
pop_peak = E_df_fitt.query('name == "center"').loc[:,('value', 'spot')].set_index('spot')
pop_peak.columns = ['E']
pop_peak['S'] = S_df_fitt.query('name == "center"').loc[:,('value', 'spot')].set_index('spot')
pop_peak['Su'] = Su_df_fitt.query('name == "center"').loc[:,('value', 'spot')].set_index('spot')
pop_peak['num_bursts'] = dx.num_bursts
pop_peak['E_sigma'] = E_df_fitt.query('name == "sigma"').loc[:,('value', 'spot')].set_index('spot')
pop_peak['S_sigma'] = S_df_fitt.query('name == "sigma"').loc[:,('value', 'spot')].set_index('spot')
pop_peak['Su_sigma'] = Su_df_fitt.query('name == "sigma"').loc[:,('value', 'spot')].set_index('spot')
pop_peak['E_err'] = pop_peak['E_sigma'] / np.sqrt(pop_peak['num_bursts'])
pop_peak['S_err'] = pop_peak['S_sigma'] / np.sqrt(pop_peak['num_bursts'])
pop_peak['Su_err'] = pop_peak['Su_sigma'] / np.sqrt(pop_peak['num_bursts'])
pop_peak.loc[skip_ch, (c for c in pop_peak.columns if c != 'num_bursts')] = np.nan
# Add categorical column for pixel grouping
# spotv is a global var from heatmatp.py
pixel_groups = {
'center': spotsv[4:8],
'top': spotsv[:4],
'bottom': spotsv[-4:]}
cat = {}
for i, row in pop_peak.iterrows():
for group, group_vals in pixel_groups.items():
if i in group_vals:
cat[i] = group
pop_peak['Pixel'] = pd.Series(cat, dtype="category")
return pop_peak
Efret_fitg = br.glance(E_fitter.fit_res, var_names='spot')
Sfret_fitg = br.glance(S_fitter.fit_res, var_names='spot')
Sufret_fitg = br.glance(Su_fitter.fit_res, var_names='spot')
Efret_fitg.head(3)
model | method | num_params | num_data_points | chisqr | redchi | AIC | BIC | num_func_eval | success | message | spot | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | Model(gaussian) | Nelder-Mead | 3 | 46 | 38.595662 | 0.897574 | -2.073069 | 3.412855 | 365 | True | Optimization terminated successfully. | 0 |
1 | Model(gaussian) | Nelder-Mead | 3 | 46 | 101.887264 | 2.369471 | 42.580375 | 48.066300 | 265 | True | Optimization terminated successfully. | 1 |
2 | Model(gaussian) | Nelder-Mead | 3 | 46 | 745.645867 | 17.340602 | 134.138032 | 139.623956 | 302 | True | Optimization terminated successfully. | 2 |
Sfret_fitg.head(3)
model | method | num_params | num_data_points | chisqr | redchi | AIC | BIC | num_func_eval | success | message | spot | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | Model(gaussian) | Nelder-Mead | 3 | 46 | 100.012545 | 2.325873 | 41.726095 | 47.212019 | 395 | True | Optimization terminated successfully. | 0 |
1 | Model(gaussian) | Nelder-Mead | 3 | 46 | 222.462245 | 5.173541 | 78.501336 | 83.987260 | 300 | True | Optimization terminated successfully. | 1 |
2 | Model(gaussian) | Nelder-Mead | 3 | 46 | 334.385396 | 7.776405 | 97.248029 | 102.733953 | 245 | True | Optimization terminated successfully. | 2 |
Sufret_fitg.head(3)
model | method | num_params | num_data_points | chisqr | redchi | AIC | BIC | num_func_eval | success | message | spot | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | Model(gaussian) | Nelder-Mead | 3 | 46 | 8.994937 | 0.209185 | -69.071060 | -63.585136 | 267 | True | Optimization terminated successfully. | 0 |
1 | Model(gaussian) | Nelder-Mead | 3 | 46 | 508.235528 | 11.819431 | 116.505965 | 121.991889 | 239 | True | Optimization terminated successfully. | 1 |
2 | Model(gaussian) | Nelder-Mead | 3 | 46 | 652.057617 | 15.164131 | 127.968610 | 133.454535 | 268 | True | Optimization terminated successfully. | 2 |
Efret_fitt = br.tidy(E_fitter.fit_res, var_names='spot')
Sfret_fitt = br.tidy(S_fitter.fit_res, var_names='spot')
Sufret_fitt = br.tidy(Su_fitter.fit_res, var_names='spot')
Efret_fitt.head()
name | value | min | max | vary | expr | stderr | init_value | spot | |
---|---|---|---|---|---|---|---|---|---|
0 | amplitude | 2.615994 | -inf | inf | True | None | NaN | 1.0 | 0 |
1 | center | 0.165646 | 0.000000 | 1.000000 | True | None | NaN | 0.5 | 0 |
2 | fwhm | 0.120900 | -inf | inf | False | 2.3548200*sigma | NaN | NaN | 0 |
3 | height | 20.327208 | -inf | inf | False | 0.3989423*amplitude/max(1.e-15, sigma) | NaN | NaN | 0 |
4 | sigma | 0.051342 | 0.030000 | 0.200000 | True | None | NaN | 0.1 | 0 |
FRET_peak = make_ES_fit_dataframe(Efret_fitt, Sfret_fitt, Sufret_fitt, ds_fret)
FRET_peak.head()
E | S | Su | num_bursts | E_sigma | S_sigma | Su_sigma | E_err | S_err | Su_err | Pixel | |
---|---|---|---|---|---|---|---|---|---|---|---|
spot | |||||||||||
0 | 0.165646 | 0.857917 | 0.759625 | 90 | 0.051342 | 0.066490 | 0.044569 | 0.005412 | 0.007009 | 0.004698 | bottom |
1 | 0.157866 | 0.843901 | 0.743290 | 313 | 0.061516 | 0.078325 | 0.061728 | 0.003477 | 0.004427 | 0.003489 | bottom |
2 | 0.168479 | 0.843876 | 0.739681 | 470 | 0.060791 | 0.083696 | 0.062733 | 0.002804 | 0.003861 | 0.002894 | bottom |
3 | 0.168380 | 0.831501 | 0.730692 | 637 | 0.069229 | 0.081629 | 0.070032 | 0.002743 | 0.003234 | 0.002775 | bottom |
4 | 0.177619 | 0.809611 | 0.706921 | 543 | 0.072249 | 0.088955 | 0.075574 | 0.003101 | 0.003817 | 0.003243 | center |
binwidth = 0.02
bins= np.arange(-0.1, 1.1, binwidth) + 0.5*binwidth
labels = sorted(set(FRET_peak.Pixel))
Efret_groups = [np.array(FRET_peak.E[FRET_peak.Pixel == grp]) for grp in labels]
plt.hist(Efret_groups, bins=bins, range=(bins.min(), bins.max()), histtype='bar', stacked=True);
plt.xlabel('E')
plt.ylabel('# Spots')
plt.xlim(-0.1, 1.1)
plt.title('Distributions of FRET peak positions')
plt.legend(labels, title='Pixel group')
print('FRET pop: E mean [min, max] = %5.3f [%5.3f, %5.3f], σ = %5.3f' %
(FRET_peak.E.mean(), FRET_peak.E.min(), FRET_peak.E.max(), FRET_peak.E.std()))
savefig('FRET_peaks_distributions_grouped.png')
FRET pop: E mean [min, max] = 0.154 [0.103, 0.178], σ = 0.014 Saved: figures/2017-05-23_04_22d_FRET_peaks_distributions_grouped.png.png Saved hires: figures/2017-05-23_04_22d_FRET_peaks_distributions_grouped_highres.png
binwidth = 0.02
bins= np.arange(-0.1, 1.1, binwidth) + 0.5*binwidth
labels = sorted(set(FRET_peak.Pixel))
Sfret_groups = [np.array(FRET_peak.S[FRET_peak.Pixel == grp]) for grp in labels]
plt.hist(Sfret_groups, bins=bins, range=(bins.min(), bins.max()), histtype='bar', stacked=True);
plt.xlabel('S')
plt.ylabel('# Spots')
plt.xlim(-0.1, 1.1)
plt.title('Distributions of S peak positions')
plt.legend(labels, title='Pixel group')
print('FRET pop: S mean [min, max] = %5.3f [%5.3f, %5.3f], σ = %5.3f' %
(FRET_peak.S.mean(), FRET_peak.S.min(), FRET_peak.S.max(), FRET_peak.S.std()))
savefig('S_peaks_distributions_grouped.png')
FRET pop: S mean [min, max] = 0.836 [0.810, 0.872], σ = 0.015 Saved: figures/2017-05-23_04_22d_S_peaks_distributions_grouped.png.png Saved hires: figures/2017-05-23_04_22d_S_peaks_distributions_grouped_highres.png
binwidth = 0.02
bins= np.arange(-0.1, 1.1, binwidth) + 0.5*binwidth
labels = sorted(set(FRET_peak.Pixel))
Sufret_groups = [np.array(FRET_peak.Su[FRET_peak.Pixel == grp]) for grp in labels]
plt.hist(Sufret_groups, bins=bins, range=(bins.min(), bins.max()), histtype='bar', stacked=True);
plt.xlabel('$S_u$')
plt.ylabel('# Spots')
plt.xlim(-0.1, 1.1)
plt.title('Distributions of $S_u$ peak positions')
plt.legend(labels, title='Pixel group')
print('FRET pop: Su mean [min, max] = %5.3f [%5.3f, %5.3f], σ = %5.3f' %
(FRET_peak.Su.mean(), FRET_peak.Su.min(), FRET_peak.Su.max(), FRET_peak.Su.std()))
savefig('Su_peaks_distributions_grouped.png')
FRET pop: Su mean [min, max] = 0.740 [0.707, 0.787], σ = 0.017 Saved: figures/2017-05-23_04_22d_Su_peaks_distributions_grouped.png.png Saved hires: figures/2017-05-23_04_22d_Su_peaks_distributions_grouped_highres.png
num_bursts = FRET_peak.num_bursts
x = FRET_peak.index
y = FRET_peak.E
yerr = np.array(Efret_fitt.query('name == "sigma"')['value']) / np.sqrt(num_bursts)
sort = num_bursts.argsort()
plt.figure(figsize=(10, 4))
x = FRET_peak.index / FRET_peak.index.max() * 100
y = FRET_peak.E
yerr = FRET_peak.E_err
plt.errorbar(x, y.loc[sort], yerr.loc[sort], lw=0, elinewidth=2, marker='o', label='Fitted E')
plt.xlabel('Fraction of spots (%)')
plt.ylabel('Fitted E')
plt.ylim(0, 1)
plt.legend(frameon=False, loc='upper left')
delta = 0.05
for ich, (spot, yi, xi) in enumerate(zip(np.arange(48)[sort], y.loc[sort], x)):
delta *= -1
plt.annotate('%d' % spot, (xi, yi + delta), ha='center', va='center')
ax2 = plt.twinx()
ax2.plot(x, num_bursts[sort], color='C1', ls='--', label='# bursts')
ax2.set_ylabel('# Bursts')
ax2.grid(True)
plt.legend(frameon=False, loc='upper right');
savefig('E fit vs num_bursts')
Saved: figures/2017-05-23_04_22d_E_fit_vs_num_bursts.png Saved hires: figures/2017-05-23_04_22d_E_fit_vs_num_bursts_highres.png
y = np.array(FRET_peak.E).reshape(4, 12)
yerr = np.array(FRET_peak.E_err).reshape(4, 12)
x = np.arange(12)
fig, ax = plt.subplots(1, 1, figsize=(10, 6), squeeze=True, sharex=True, sharey=True)
for irow in range(4):
ax.errorbar(x, y[irow], yerr[irow], lw=0, elinewidth=2, marker='o', label='Pixel Row %d' % irow)
ax.set_xticks(range(0, 12, 1))
ax.set_ylabel('Fitted E')
ax.set_xlabel('Pixel Column')
plt.legend()
plt.grid(True, axis='y')
plt.title('Fitted E vs Pixel Column')
savefig('Fitted E vs Pixel Column')
Saved: figures/2017-05-23_04_22d_Fitted_E_vs_Pixel_Column.png Saved hires: figures/2017-05-23_04_22d_Fitted_E_vs_Pixel_Column_highres.png
fig, ax = plt.subplots(figsize=(5,5))
ax.plot(FRET_peak.E, FRET_peak.S, lw=0, marker='+', mew=1, label='FRET peak')
#ax.plot(FRET_peak.E, FRET_peak.S, lw=0, marker='o', label='FRET peak')
ax.set_xlim(0, 1)
ax.set_ylim(0, 1)
(0, 1)
ddo = d.copy()
Deep copy executed.
ddo.burst_search(min_rate_cps=50e3, pax=True, ph_sel=Ph_sel(Dex='Dem', Aex='Dem'))
- Performing burst search (verbose=False) ... - Fixing burst data to refer to ph_times_m ... [DONE] [DONE] - Calculating burst periods ...[DONE] - Counting D and A ph and calculating FRET ... - Applying background correction. [DONE Counting D/A]
ds_dox = ddo.select_bursts(select_bursts.size, ph_sel=Ph_sel(Dex='Dem', Aex='Dem'),
th1 = 40)
ds_do2 = ds_dox.select_bursts(select_bursts.size,
ph_sel=Ph_sel(Aex='Aem'), naa_aexonly=False,
th2=4, th1=-1000)
Su = [(nd + na)/(nd + na + naa) for nd, na, naa in zip(ds_do2.nd, ds_do2.na, ds_do2.naa)]
ds_do2.add(Su=Su)
g = alex_jointplot(ds_do2, i=30, S_name='Su', vmax_fret=False)
g.ax_joint.set_ylabel('$S_u$')
/home/anto/src/FRETBursts/fretbursts/mfit.py:329: RuntimeWarning: invalid value encountered in true_divide self.hist_pdf /= self.hist_counts.sum(1)[:, np.newaxis]
<matplotlib.text.Text at 0x7f7eeba51240>
ax = dplot(ds_do2, hexbin_alex, vmax=30, S_name='Su', skip_ch=skip_ch,
wspace=0, hspace=0, gridsize=60, title='bottom left', title_bg=False);
ax0 = ax[0, 0]
fig = ax0.figure
cax = fig.add_axes([0.97, 0.25, 0.01, 0.5])
plt.setp(ax[:, 0], ylabel='$S_u$');
plt.colorbar(cax=cax)
ax0.set_xticks([0, 0.5, 1])
ax0.set_yticks([0, 0.5, 1]);
ax0.set_xlim(-0.2, 1)
ax0.set_ylim(0, 1.2);
savefig('48spot alex hist Su DonlyBS')
Saved: figures/2017-05-23_04_22d_48spot_alex_hist_Su_DonlyBS.png Saved hires: figures/2017-05-23_04_22d_48spot_alex_hist_Su_DonlyBS_highres.png
ds_do = ds_do2
E_fitter = bext.bursts_fitter(ds_do, 'E', skip_ch=skip_ch)
model = mfit.factory_gaussian()
model.set_param_hint('center', value=0.05, min=-0.2, max=0.4)
model.set_param_hint('sigma', value=0.1, min=0.03, max=0.2)
model.print_param_hints()
Name Value Min Max Vary Expr amplitude 1 -inf inf True center 0.05 -0.2 0.4 True fwhm nan -inf inf True 2.3548200*sigma height nan -inf inf True 0.3989423*amplitude/max(1.e-15, sigma) sigma 0.1 0.03 0.2 True
E_fitter.fit_histogram(model, pdf=False, method='mealder')
/home/anto/src/FRETBursts/fretbursts/mfit.py:329: RuntimeWarning: invalid value encountered in true_divide self.hist_pdf /= self.hist_counts.sum(1)[:, np.newaxis]
ax = dplot(ds_do, hist_fret, pdf=False, show_model=True, skip_ch=skip_ch);
plt.setp(ax[:, 0], ylabel='# Bursts')
savefig('DO pop - E hist fit')
Saved: figures/2017-05-23_04_22d_DO_pop_-_E_hist_fit.png Saved hires: figures/2017-05-23_04_22d_DO_pop_-_E_hist_fit_highres.png
S_fitter = bext.bursts_fitter(ds_do, 'S', skip_ch=skip_ch)
model = mfit.factory_gaussian()
model.set_param_hint('center', value=0.9, min=0, max=1.1)
model.set_param_hint('sigma', value=0.1, min=0.03, max=0.3)
model.print_param_hints()
Name Value Min Max Vary Expr amplitude 1 -inf inf True center 0.9 0 1.1 True fwhm nan -inf inf True 2.3548200*sigma height nan -inf inf True 0.3989423*amplitude/max(1.e-15, sigma) sigma 0.1 0.03 0.3 True
S_fitter.fit_histogram(model, pdf=False, method='mealder')
/home/anto/src/FRETBursts/fretbursts/mfit.py:329: RuntimeWarning: invalid value encountered in true_divide self.hist_pdf /= self.hist_counts.sum(1)[:, np.newaxis]
ax = dplot(ds_do, hist_S, pdf=False, show_model=True, skip_ch=skip_ch);
plt.setp(ax[:, 0], ylabel='# Bursts')
savefig('DO pop - S hist fit');
Saved: figures/2017-05-23_04_22d_DO_pop_-_S_hist_fit.png Saved hires: figures/2017-05-23_04_22d_DO_pop_-_S_hist_fit_highres.png
Su_fitter = bext.bursts_fitter(ds_do, 'Su', skip_ch=skip_ch)
model = mfit.factory_gaussian()
model.set_param_hint('center', value=0.9, min=0, max=1.2)
model.set_param_hint('sigma', value=0.1, min=0.03, max=0.2)
model.print_param_hints()
Name Value Min Max Vary Expr amplitude 1 -inf inf True center 0.9 0 1.2 True fwhm nan -inf inf True 2.3548200*sigma height nan -inf inf True 0.3989423*amplitude/max(1.e-15, sigma) sigma 0.1 0.03 0.2 True
Su_fitter.fit_histogram(model, pdf=False, method='mealder')
/home/anto/src/FRETBursts/fretbursts/mfit.py:329: RuntimeWarning: invalid value encountered in true_divide self.hist_pdf /= self.hist_counts.sum(1)[:, np.newaxis]
ax = dplot(ds_do, hist_burst_data, data_name='Su', pdf=False, show_model=True, skip_ch=skip_ch);
plt.setp(ax[:, 0], ylabel='# Bursts')
savefig('DO pop - Su hist fit');
Saved: figures/2017-05-23_04_22d_DO_pop_-_Su_hist_fit.png Saved hires: figures/2017-05-23_04_22d_DO_pop_-_Su_hist_fit_highres.png
Edo_fitg = br.glance(E_fitter.fit_res, var_names='spot')
Sdo_fitg = br.glance(S_fitter.fit_res, var_names='spot')
Sudo_fitg = br.glance(Su_fitter.fit_res, var_names='spot')
Edo_fitg.head(3)
model | method | num_params | num_data_points | chisqr | redchi | AIC | BIC | num_func_eval | success | message | spot | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | Model(gaussian) | Nelder-Mead | 3 | 46 | 2085.804981 | 48.507093 | 181.456362 | 186.942286 | 182 | True | Optimization terminated successfully. | 0 |
1 | Model(gaussian) | Nelder-Mead | 3 | 46 | 4525.757869 | 105.250183 | 217.089351 | 222.575275 | 201 | True | Optimization terminated successfully. | 1 |
2 | Model(gaussian) | Nelder-Mead | 3 | 46 | 5234.973709 | 121.743575 | 223.785882 | 229.271807 | 218 | True | Optimization terminated successfully. | 2 |
Sdo_fitg.head(3)
model | method | num_params | num_data_points | chisqr | redchi | AIC | BIC | num_func_eval | success | message | spot | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | Model(gaussian) | Nelder-Mead | 3 | 46 | 537.136139 | 12.491538 | 119.050068 | 124.535993 | 148 | True | Optimization terminated successfully. | 0 |
1 | Model(gaussian) | Nelder-Mead | 3 | 46 | 1872.626954 | 43.549464 | 176.496981 | 181.982906 | 193 | True | Optimization terminated successfully. | 1 |
2 | Model(gaussian) | Nelder-Mead | 3 | 46 | 1521.825278 | 35.391286 | 166.955120 | 172.441044 | 247 | True | Optimization terminated successfully. | 2 |
Sudo_fitg.head(3)
model | method | num_params | num_data_points | chisqr | redchi | AIC | BIC | num_func_eval | success | message | spot | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | Model(gaussian) | Nelder-Mead | 3 | 46 | 94.436148 | 2.196189 | 39.086996 | 44.572921 | 173 | True | Optimization terminated successfully. | 0 |
1 | Model(gaussian) | Nelder-Mead | 3 | 46 | 1559.830014 | 36.275117 | 168.089774 | 173.575698 | 203 | True | Optimization terminated successfully. | 1 |
2 | Model(gaussian) | Nelder-Mead | 3 | 46 | 539.470375 | 12.545823 | 119.249538 | 124.735462 | 210 | True | Optimization terminated successfully. | 2 |
Edo_fitt = br.tidy(E_fitter.fit_res, var_names='spot')
Sdo_fitt = br.tidy(S_fitter.fit_res, var_names='spot')
Sudo_fitt = br.tidy(Su_fitter.fit_res, var_names='spot')
Edo_fitt.head()
name | value | min | max | vary | expr | stderr | init_value | spot | |
---|---|---|---|---|---|---|---|---|---|
0 | amplitude | 7.708896 | -inf | inf | True | None | NaN | 1.00 | 0 |
1 | center | 0.028217 | -0.200000 | 0.400000 | True | None | NaN | 0.05 | 0 |
2 | fwhm | 0.092955 | -inf | inf | False | 2.3548200*sigma | NaN | NaN | 0 |
3 | height | 77.908654 | -inf | inf | False | 0.3989423*amplitude/max(1.e-15, sigma) | NaN | NaN | 0 |
4 | sigma | 0.039474 | 0.030000 | 0.200000 | True | None | NaN | 0.10 | 0 |
DO_peak = make_ES_fit_dataframe(Edo_fitt, Sdo_fitt, Sudo_fitt, ds_do)
DO_peak.head()
E | S | Su | num_bursts | E_sigma | S_sigma | Su_sigma | E_err | S_err | Su_err | Pixel | |
---|---|---|---|---|---|---|---|---|---|---|---|
spot | |||||||||||
0 | 0.028217 | 0.997646 | 0.956716 | 286 | 0.039474 | 0.054290 | 0.047883 | 0.002334 | 0.003210 | 0.002831 | bottom |
1 | 0.031751 | 1.001682 | 0.957605 | 575 | 0.036241 | 0.062280 | 0.048894 | 0.001511 | 0.002597 | 0.002039 | bottom |
2 | 0.040697 | 1.007145 | 0.959760 | 731 | 0.044569 | 0.062732 | 0.042136 | 0.001648 | 0.002320 | 0.001558 | bottom |
3 | 0.038256 | 1.006073 | 0.960267 | 841 | 0.042194 | 0.056118 | 0.044878 | 0.001455 | 0.001935 | 0.001548 | bottom |
4 | 0.039505 | 1.004416 | 0.958228 | 736 | 0.045398 | 0.062950 | 0.042887 | 0.001673 | 0.002320 | 0.001581 | center |
binwidth = 0.02
bins= np.arange(-0.1, 1.1, binwidth) + 0.5*binwidth
labels = sorted(set(DO_peak.Pixel))
Edo_groups = [np.array(DO_peak.E[DO_peak.Pixel == grp]) for grp in labels]
plt.hist(Edo_groups, bins=bins, histtype='bar', range=(bins.min(), bins.max()), stacked=True);
plt.hist(Efret_groups, bins=bins, histtype='bar', range=(bins.min(), bins.max()), stacked=True, color=('C0', 'C1', 'C2'));
plt.xlabel('E')
plt.ylabel('# Spots')
plt.xlim(-0.1, 1.1)
plt.title('Distributions of E peak centers for FRET and D-only populations')
plt.legend(labels)
print('DO pop : E mean [min, max] = %5.3f [%5.3f, %5.3f], σ = %5.3f' %
(DO_peak.E.mean(), DO_peak.E.min(), DO_peak.E.max(), DO_peak.E.std()))
print('FRET pop: E mean [min, max] = %5.3f [%5.3f, %5.3f], σ = %5.3f' %
(FRET_peak.E.mean(), FRET_peak.E.min(), FRET_peak.E.max(), FRET_peak.E.std()))
savefig('FRET_vs_DO_peaks_distributions_grouped.png')
DO pop : E mean [min, max] = 0.036 [0.022, 0.048], σ = 0.005 FRET pop: E mean [min, max] = 0.154 [0.103, 0.178], σ = 0.014 Saved: figures/2017-05-23_04_22d_FRET_vs_DO_peaks_distributions_grouped.png.png Saved hires: figures/2017-05-23_04_22d_FRET_vs_DO_peaks_distributions_grouped_highres.png
binwidth = 0.02
bins= np.arange(-0.1, 1.1, binwidth) + 0.5*binwidth
labels = sorted(set(DO_peak.Pixel))
Sudo_groups = [np.array(DO_peak.Su[DO_peak.Pixel == grp]) for grp in labels]
plt.hist(Sudo_groups, bins=bins, histtype='bar', range=(bins.min(), bins.max()), stacked=True);
plt.hist(Sufret_groups, bins=bins, histtype='bar', range=(bins.min(), bins.max()), stacked=True, color=('C0', 'C1', 'C2'));
plt.xlabel('$S_u$')
plt.ylabel('# Spots')
plt.xlim(-0.1, 1.1)
plt.title('Distributions of $S_u$ peak centers for FRET and D-only populations')
plt.legend(labels)
print('FRET pop: Su mean [min, max] = %5.3f [%5.3f, %5.3f], σ = %5.3f' %
(FRET_peak.Su.mean(), FRET_peak.Su.min(), FRET_peak.Su.max(), FRET_peak.Su.std()))
print('DO : Su mean [min, max] = %5.3f [%5.3f, %5.3f], σ = %5.3f' %
(DO_peak.Su.mean(), DO_peak.Su.min(), DO_peak.Su.max(), DO_peak.Su.std()))
savefig('Su_peaks_distributions_grouped.png')
FRET pop: Su mean [min, max] = 0.740 [0.707, 0.787], σ = 0.017 DO : Su mean [min, max] = 0.962 [0.952, 0.968], σ = 0.003 Saved: figures/2017-05-23_04_22d_Su_peaks_distributions_grouped.png.png Saved hires: figures/2017-05-23_04_22d_Su_peaks_distributions_grouped_highres.png
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(FRET_peak.E, FRET_peak.S, lw=0, marker='+', mew=1, label='FRET peak')
ax.plot(DO_peak.E, DO_peak.S, lw=0, marker='x', mew=1, label='D-only peak')
ax.set_xlim(-0.1, 1)
ax.set_ylim(0.4, 1.1)
plt.setp(ax, xlabel='E', ylabel='S')
plt.title('FRET vs DO fitted Epr-S peak position')
plt.legend()
savefig('FRET vs DO fitted Epr-S peak position')
Saved: figures/2017-05-23_04_22d_FRET_vs_DO_fitted_Epr-S_peak_position.png Saved hires: figures/2017-05-23_04_22d_FRET_vs_DO_fitted_Epr-S_peak_position_highres.png
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(FRET_peak.E, FRET_peak.Su, lw=0, marker='+', mew=1, label='FRET peak')
ax.plot(DO_peak.E, DO_peak.Su, lw=0, marker='x', mew=1, label='D-only peak')
ax.set_xlim(-0.1, 1)
ax.set_ylim(0.4, 1.1)
plt.setp(ax, xlabel='E', ylabel='$S_u$')
plt.title('Fitted $E_{PR}$-$S_u$ peak centers in different spots')
plt.legend()
savefig('FRET vs DO fitted Epr-Su peak position')
Saved: figures/2017-05-23_04_22d_FRET_vs_DO_fitted_Epr-Su_peak_position.png Saved hires: figures/2017-05-23_04_22d_FRET_vs_DO_fitted_Epr-Su_peak_position_highres.png
def figure48(figsize=(18, 6), annotate=True,
text_fmt='%d', text_x=0.9, text_y=0.9, text_data=spotsh, text_va='top', text_ha='right'):
fig, axes = plt.subplots(4, 12, figsize=figsize, sharex=True, sharey=True)
plt.subplots_adjust(hspace=0, wspace=0)
if annotate:
for irow, axrow in enumerate(axes):
for icol, ax in enumerate(axrow):
ax.text(text_x, text_y, text_fmt % text_data[irow, icol],
transform=ax.transAxes, ha=text_ha, va=text_va)
return fig, axes
err_kws = dict(capsize=6, capthick=1.5, marker='o', ms=6)
x = np.array(FRET_peak.E).reshape(4, 12)
xerr = np.array(FRET_peak.E_sigma).reshape(4, 12)
y = np.array(FRET_peak.Su).reshape(4, 12)
yerr = np.array(FRET_peak.Su_sigma).reshape(4, 12)
xmean = np.nanmean(x)
ymean = np.nanmean(y)
delta = np.max([np.nanmax(yerr), np.nanmax(xerr)])*1.5
fig, axes = figure48(text_x=0.80, text_y=0.85, text_ha='right')
for irow, axrow in enumerate(axes):
for icol, ax in enumerate(axrow):
ax.errorbar(x[irow, icol], y[irow, icol], xerr=xerr[irow,icol], yerr=yerr[irow,icol], **err_kws)
ax.plot(xmean, ymean, marker='o', zorder=10)
#sns.despine(left=True, bottom=True)
#plt.setp(axes[0], xlim=(np.round(xmean - delta, 1), np.round(xmean + delta, 1)));
#plt.setp(axes[:,0], ylim=(ymean - delta, ymean + delta));
plt.setp(axes[0], xlim=(0, 0.49));
plt.setp(axes[:,0], ylim=(0.4, 0.9));
plt.setp(axes[-1], xlabel='$E_{PR}$')
plt.setp(axes[:, 0], ylabel='$S_u$');
# for ax in axes[-1]:
# sns.despine(bottom=False, left=True, trim=True, ax=ax)
# for ax in axes[:, 0]:
# sns.despine(left=False, bottom=True, trim=True, ax=ax)
# sns.despine(left=False, bottom=False, trim=True, ax=axes[-1,0])
# for ax in axes[:-1].ravel():
# ax.tick_params(bottom=False)
# for ax in axes[:, 1:].ravel():
# ax.tick_params(left=False)
plt.suptitle('Deviation for fitted $E_{PR}$ and $S_u$ from mean values for the 48 spots (sigma peak)', fontsize=14, va='top');
savefig('48-spot Epr-Su peak position deviation from mean')
Saved: figures/2017-05-23_04_22d_48-spot_Epr-Su_peak_position_deviation_from_mean.png Saved hires: figures/2017-05-23_04_22d_48-spot_Epr-Su_peak_position_deviation_from_mean_highres.png
For convenience we reporte the formula for the relative-$\gamma$ factor $\chi_{ch}$ and the expression of $E$ as a function of $E_{PR}$:
$$ \chi_{ch} = \frac{\frac{1}{\langle E_{PR\,ch} \rangle_{ch}} - 1}{\frac{1}{E_{PR\,ch}} - 1} $$$$ E = \frac{1}{1 + \gamma \left( \frac{1}{E_{PR}} - 1 \right)} $$chi_ch = (1/FRET_peak.E.mean() - 1) / (1/FRET_peak.E - 1)
E_corr = 1 / (1 + chi_ch * (1/FRET_peak.E - 1))
E_corr.head()
spot 0 0.153842 1 0.153842 2 0.153842 3 0.153842 4 0.153842 Name: E, dtype: float64
assert (E_corr.max() - E_corr.min()) < 1e-9 # Test that E_corr are all "equal"
ds_fret.num_bursts
array([ 90, 313, 470, 637, 543, 663, 679, 778, 569, 588, 487, 297, 0, 0, 491, 614, 565, 627, 630, 409, 634, 507, 480, 395, 117, 367, 375, 498, 646, 718, 690, 836, 584, 560, 564, 342, 109, 270, 415, 352, 530, 570, 595, 628, 611, 534, 391, 287])
ds_fret.chi_ch = np.ones(48)
bursts0 = bext.burst_data(ds_fret)
bursts0.shape
(23055, 13)
%debug
ERROR:root:No traceback has been produced, nothing to debug.
ds_fret.chi_ch = chi_ch
heatmap48(chi_ch, annot=True, fmt='.2f')
chi_ch.hist(bins=np.arange(0.5, 1.5, 0.025))
<matplotlib.axes._subplots.AxesSubplot at 0x7f7eebcbbd68>
bursts = bext.burst_data(ds_fret)
bursts.shape
(23055, 13)
plt.rcParams['font.size'] = 12
kw = dict(bins=np.arange(-0.1,1.2, 0.02), histtype='step', lw=2, range=(-0.1, 1.1))
plt.hist(bursts0.E, label='raw', **kw);
plt.hist(bursts.E, label='calibrated', **kw);
plt.legend();
plt.title('48-spot E histogram, raw vs per-ch calibration');
plt.xlabel('$E_{PR}$')
savefig('48-spot calibration')
plt.rcParams['font.size'] = 14
Saved: figures/2017-05-23_04_22d_48-spot_calibration.png Saved hires: figures/2017-05-23_04_22d_48-spot_calibration_highres.png
bursts.head()
E | S | na | naa | nar | nd | nda | nt | size_raw | t_start | t_stop | width_ms | spot | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 0.057279 | 0.867516 | 3.332915 | 12.219002 | 3.332915 | 55.384825 | 45.084430 | 116.021172 | 122 | 16.305300 | 16.306744 | 1.444075 | 0 |
1 | 0.264611 | 0.731624 | 9.515794 | 22.707256 | 9.515794 | 19.827619 | 28.609576 | 80.660244 | 85 | 16.306873 | 16.307922 | 1.048187 | 0 |
2 | 0.222807 | 0.804819 | 9.488131 | 19.815531 | 9.488131 | 22.918953 | 37.699443 | 89.922059 | 94 | 38.477065 | 38.478054 | 0.988987 | 0 |
3 | 0.140544 | 0.849402 | 6.596321 | 14.917730 | 6.596321 | 40.974772 | 32.907075 | 95.395897 | 99 | 42.778437 | 42.779320 | 0.882937 | 0 |
4 | 0.154250 | 0.822340 | 8.479690 | 20.356280 | 8.479690 | 38.636711 | 46.520526 | 113.993207 | 119 | 77.963292 | 77.964613 | 1.320600 | 0 |
ich = 33
chi_ch[ich]
0.9467056663475103
ds_fret.chi_ch = np.ones(48)
E0 = ds_fret.E[ich]
ds_fret.chi_ch = chi_ch
E = ds_fret.E[ich]
kw = dict(bins=np.arange(-0.1,1.1, 0.0333), histtype='step', lw=2)
plt.hist(E0, label='raw', **kw);
plt.hist(E, label='calibrated', **kw);
plt.legend();
plt.title('Spot %d ($\chi_i$ = %.2f) E histogram, raw vs calibrated' % (ich, chi_ch[ich]));
ich
33
E0i = bursts0.E.loc[bursts0.spot == ich]
Ei = bursts.E.loc[bursts.spot == ich]
E0a = bursts0.E.loc[~pd.isnull(bursts0.E)]
Ea = bursts.E.loc[~pd.isnull(bursts.E)]
plt.rcParams['font.size'] = 12
fig, (ax0, ax1) = plt.subplots(1, 2, figsize=(10, 4))
kw = dict(bins=np.arange(-0.1,1.1, 0.0333), histtype='step', lw=2)
ax0.hist(E0a, label='raw', **kw);
ax0.hist(Ea, label='calibrated', **kw);
ax1.hist(E0i, label='raw', **kw);
ax1.hist(Ei, label='calibrated', **kw);
ax1.legend();
ax1.set_title('Spot %d, $\chi_i$ = %.2f' % (ich, chi_ch[ich]));
ax0.set_title('All Spots');
ax1.set_title('Spot %d, $\chi_i$ = %.2f' % (ich, chi_ch[ich]));
plt.rcParams['font.size'] = 14
(bursts.spot == ich).sum()
560
bursts.spot.shape
(23055,)
bursts0.E.loc[bursts.spot == ich].shape
(560,)
E.shape, bursts.E.loc[bursts.spot == ich].shape
((560,), (560,))
bursts.sample(10).sort_index()
E | S | na | naa | nar | nd | nda | nt | size_raw | t_start | t_stop | width_ms | spot | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|
5206 | 0.122435 | 0.822022 | 8.317135 | 23.024978 | 8.317135 | 50.913461 | 74.599207 | 156.854782 | 162 | 1076.893996 | 1076.895564 | 1.568288 | 9 |
5739 | 0.152062 | 0.813591 | 7.223231 | 18.106822 | 7.223231 | 44.040737 | 50.723958 | 120.094748 | 125 | 1156.250903 | 1156.252274 | 1.371250 | 10 |
6203 | 0.131823 | 0.998516 | 10.271580 | 10.387411 | 10.271580 | 65.741473 | 53.615952 | 140.016415 | 145 | 199.695614 | 199.696936 | 1.321937 | 14 |
13248 | 0.119972 | 0.945672 | 9.327317 | 13.793714 | 9.327317 | 63.420804 | 67.039989 | 153.581823 | 160 | 862.688717 | 862.690066 | 1.349425 | 28 |
15407 | 0.133347 | 0.731858 | 5.588891 | 20.945009 | 5.588891 | 40.280894 | 29.133797 | 95.948590 | 99 | 851.615705 | 851.616502 | 0.796375 | 31 |
16413 | 0.163287 | 1.004104 | 10.483836 | 10.221417 | 10.483836 | 65.806568 | 47.684395 | 134.196217 | 139 | 207.322591 | 207.323695 | 1.103887 | 33 |
17643 | 0.131583 | 0.727874 | 4.681963 | 17.984735 | 4.681963 | 29.422136 | 39.233076 | 91.321911 | 94 | 872.385366 | 872.386152 | 0.786225 | 35 |
18201 | 0.067742 | 0.791706 | 3.696733 | 18.053993 | 3.696733 | 45.181154 | 45.876927 | 112.808807 | 121 | 189.054209 | 189.055209 | 1.000375 | 38 |
18586 | 0.164623 | 0.822696 | 8.535910 | 19.710716 | 8.535910 | 41.664599 | 38.523864 | 108.435088 | 113 | 114.451531 | 114.452541 | 1.009512 | 39 |
22535 | 0.229643 | 0.777919 | 7.716824 | 17.309964 | 7.716824 | 28.521615 | 24.447046 | 77.995450 | 80 | 492.644958 | 492.645617 | 0.658575 | 46 |
bursts.t_start.max()
1434.1370831624999
df = d.fuse_bursts(ms=0)
- - - - - CHANNEL 1 - - - - --> END Fused 3822 bursts (26.8%, 10 iter) - - - - - CHANNEL 2 - - - - --> END Fused 4297 bursts (25.0%, 9 iter) - - - - - CHANNEL 3 - - - - --> END Fused 4644 bursts (24.6%, 8 iter) - - - - - CHANNEL 4 - - - - --> END Fused 5299 bursts (25.1%, 9 iter) - - - - - CHANNEL 5 - - - - --> END Fused 5129 bursts (24.3%, 9 iter) - - - - - CHANNEL 6 - - - - --> END Fused 5991 bursts (25.4%, 8 iter) - - - - - CHANNEL 7 - - - - --> END Fused 5921 bursts (25.5%, 9 iter) - - - - - CHANNEL 8 - - - - --> END Fused 7607 bursts (26.3%, 8 iter) - - - - - CHANNEL 9 - - - - --> END Fused 5477 bursts (24.9%, 9 iter) - - - - - CHANNEL 10 - - - - --> END Fused 4554 bursts (23.2%, 8 iter) - - - - - CHANNEL 11 - - - - --> END Fused 4592 bursts (24.2%, 8 iter) - - - - - CHANNEL 12 - - - - --> END Fused 4544 bursts (24.8%, 10 iter) - - - - - CHANNEL 13 - - - - --> END Fused 45 bursts (12.7%, 5 iter) - - - - - CHANNEL 14 - - - - --> END Fused 195 bursts (16.0%, 6 iter) - - - - - CHANNEL 15 - - - - --> END Fused 4753 bursts (25.0%, 8 iter) - - - - - CHANNEL 16 - - - - --> END Fused 6514 bursts (26.2%, 10 iter) - - - - - CHANNEL 17 - - - - --> END Fused 5408 bursts (25.1%, 11 iter) - - - - - CHANNEL 18 - - - - --> END Fused 5516 bursts (24.4%, 9 iter) - - - - - CHANNEL 19 - - - - --> END Fused 5616 bursts (25.3%, 8 iter) - - - - - CHANNEL 20 - - - - --> END Fused 5550 bursts (25.4%, 8 iter) - - - - - CHANNEL 21 - - - - --> END Fused 5068 bursts (24.2%, 8 iter) - - - - - CHANNEL 22 - - - - --> END Fused 4953 bursts (24.7%, 9 iter) - - - - - CHANNEL 23 - - - - --> END Fused 4403 bursts (23.6%, 9 iter) - - - - - CHANNEL 24 - - - - --> END Fused 4459 bursts (24.9%, 9 iter) - - - - - CHANNEL 25 - - - - --> END Fused 4612 bursts (27.3%, 8 iter) - - - - - CHANNEL 26 - - - - --> END Fused 5031 bursts (26.0%, 10 iter) - - - - - CHANNEL 27 - - - - --> END Fused 4771 bursts (25.6%, 9 iter) - - - - - CHANNEL 28 - - - - --> END Fused 4951 bursts (25.2%, 9 iter) - - - - - CHANNEL 29 - - - - --> END Fused 5945 bursts (25.8%, 11 iter) - - - - - CHANNEL 30 - - - - --> END Fused 7840 bursts (25.8%, 12 iter) - - - - - CHANNEL 31 - - - - --> END Fused 5634 bursts (24.7%, 8 iter) - - - - - CHANNEL 32 - - - - --> END Fused 5829 bursts (23.9%, 7 iter) - - - - - CHANNEL 33 - - - - --> END Fused 6688 bursts (26.0%, 9 iter) - - - - - CHANNEL 34 - - - - --> END Fused 5217 bursts (25.1%, 8 iter) - - - - - CHANNEL 35 - - - - --> END Fused 7791 bursts (21.9%, 8 iter) - - - - - CHANNEL 36 - - - - --> END Fused 4257 bursts (25.0%, 10 iter) - - - - - CHANNEL 37 - - - - --> END Fused 4411 bursts (27.1%, 10 iter) - - - - - CHANNEL 38 - - - - --> END Fused 5714 bursts (26.0%, 8 iter) - - - - - CHANNEL 39 - - - - --> END Fused 6590 bursts (25.3%, 9 iter) - - - - - CHANNEL 40 - - - - --> END Fused 4703 bursts (25.6%, 8 iter) - - - - - CHANNEL 41 - - - - --> END Fused 5036 bursts (24.8%, 8 iter) - - - - - CHANNEL 42 - - - - --> END Fused 5149 bursts (24.8%, 8 iter) - - - - - CHANNEL 43 - - - - --> END Fused 5161 bursts (24.8%, 9 iter) - - - - - CHANNEL 44 - - - - --> END Fused 5322 bursts (24.7%, 9 iter) - - - - - CHANNEL 45 - - - - --> END Fused 5101 bursts (24.5%, 12 iter) - - - - - CHANNEL 46 - - - - --> END Fused 4597 bursts (24.2%, 10 iter) - - - - - CHANNEL 47 - - - - --> END Fused 4454 bursts (24.5%, 8 iter) - - - - - CHANNEL 48 - - - - --> END Fused 4246 bursts (24.7%, 8 iter) - Counting D and A ph and calculating FRET ... - Applying background correction. [DONE Counting D/A and FRET]
dc = d.collapse(update_gamma=False, skip_ch=skip_ch) ## check that chi_ch is set to 1 here
d.num_bursts.sum()
973696
def sel_raw_naa(d, ich=0, th1=20, th2=np.inf, gamma=1., beta=1., donor_ref=True):
assert th1 <= th2, 'th1 (%.2f) must be <= of th2 (%.2f)' % (th1, th2)
kws = dict(ich=ich, gamma=gamma, beta=beta, donor_ref=donor_ref)
naa_term = d.naa[ich]
bursts_mask = (naa_term >= th1) * (naa_term <= th2)
return bursts_mask, ''
def sel_ndex_na(d, ich=0, th1=20, th2=np.inf, gamma=1):
assert th1 <= th2, 'th1 (%.2f) must be <= of th2 (%.2f)' % (th1, th2)
size = d.nd[ich] + d.nda[ich] + d.na[ich] / gamma
bursts_mask = (size >= th1) * (size <= th2)
return bursts_mask, ''
naa_sel_kws = dict(ph_sel=Ph_sel(Aex='Aem'), gamma=1,
na_comp=False, naa_comp=False, naa_aexonly=False)
Math(d._burst_sizes_pax_formula(**naa_sel_kws))
size_sel_kws = dict(ph_sel=Ph_sel(Dex='DAem', Aex='Dem'), gamma=0.5,
na_comp=True, naa_comp=False, naa_aexonly=False)
Math(d._burst_sizes_pax_formula(**size_sel_kws))
dcs_th1 = d.select_bursts(select_bursts.size, th1=10, **naa_sel_kws)
dcs_th2 = dcs_th1.select_bursts(select_bursts.size, th1=20, **size_sel_kws)
dcs_th = dcs_th2.collapse(update_gamma=False, skip_ch=skip_ch)
Th = (10, 20, 40, 80)
Th_naa = np.array((0, 10, 20, 40))
fig, AX = plt.subplots(len(Th), len(Th_naa), figsize=(3.5*len(Th_naa), 3.5*len(Th)), sharex=True, sharey=True)
for irow, th in enumerate(Th):
for icol, th_naa in enumerate(Th_naa):
ax = AX[irow,icol]
dcs_th1 = d.select_bursts(select_bursts.size, th1=th_naa, **naa_sel_kws)
dcs_th2 = dcs_th1.select_bursts(select_bursts.size, th1=th, **size_sel_kws)
dcs_th = dcs_th2.collapse(update_gamma=False, skip_ch=skip_ch)
Su = [(nd + na)/(nd + na + naa) for nd, na, naa in zip(dcs_th.nd, dcs_th.na, dcs_th.naa)]
dcs_th.add(Su=Su)
E = dcs_th.E[0]
E_std = E[E > 0].std()
dplot(dcs_th, hexbin_alex, S_name='Su', gridsize=60, ax=ax)
ax.text(0.05, 0.05, '$th_D$ = %d, $th_A$ = %d, σ = %.3f\n#B = %d' % (th, th_naa, E_std, dcs_th.num_bursts),
va='bottom', ha='left', transform=ax.transAxes, fontsize=12)
ax.set_title('')
plt.setp(AX[:-1], xlabel='')
plt.setp(AX[:, 1:], ylabel='')
plt.subplots_adjust(hspace=0, wspace=0)
plt.text(0.5, 0.89, 'ACBS', transform=fig.transFigure, fontsize=16);
Th = (10, 20, 40, 80)
Th_naa = np.array((0, 10, 20, 40))
fig, AX = plt.subplots(len(Th), len(Th_naa), figsize=(3.5*len(Th_naa), 3.5*len(Th)), sharex=True, sharey=True)
for irow, th in enumerate(Th):
for icol, th_naa in enumerate(Th_naa):
ax = AX[irow,icol]
dcs_th1 = d.select_bursts(select_bursts.size, th1=th_naa, **naa_sel_kws)
dcs_th2 = dcs_th1.select_bursts(select_bursts.size, th1=th, **size_sel_kws)
dcs_th = dcs_th2.collapse(update_gamma=False, skip_ch=skip_ch)
dplot(dcs_th, hist_fret, ax=ax)
E = dcs_th.E[0]
E_std = E[E > 0].std()
ax.text(0.05, 0.95, '$th_D$ = %d, $th_A$ = %d, σ = %.3f\n#B = %d' % (th, th_naa, E_std, dcs_th.num_bursts),
va='top', ha='left', transform=ax.transAxes, fontsize=12)
ax.set_title('')
plt.setp(AX[:-1], xlabel='')
plt.setp(AX[:, 1:], ylabel='')
plt.subplots_adjust(hspace=0, wspace=0)
plt.text(0.5, 0.89, 'ACBS', transform=fig.transFigure, fontsize=16)
<matplotlib.text.Text at 0x7f7e88e007f0>
d2 = bext.burst_search_and_gate(d, min_rate_cps=(30e3, 15e3))
Deep copy executed. Deep copy executed. Deep copy executed. - Performing burst search (verbose=False) ... - Fixing burst data to refer to ph_times_m ... [DONE] [DONE] - Calculating burst periods ...[DONE] - Performing burst search (verbose=False) ... - Fixing burst data to refer to ph_times_m ... [DONE] [DONE] - Calculating burst periods ...[DONE] - Calculating burst periods ...[DONE] - Counting D and A ph and calculating FRET ... - Applying background correction. [DONE Counting D/A]
Th = (20, 40, 60, 80, 160)
fig, AX = plt.subplots(1, len(Th), figsize=(3.5*len(Th), 3.5), sharex=True, sharey=True)
for irow, th in enumerate(Th):
ax = AX[irow]
dcs_th2 = d2.select_bursts(select_bursts.size, th1=th, **size_sel_kws)
dcs_th = dcs_th2.collapse(update_gamma=False, skip_ch=skip_ch)
Su = [(nd + na)/(nd + na + naa) for nd, na, naa in zip(dcs_th.nd, dcs_th.na, dcs_th.naa)]
dcs_th.add(Su=Su)
dplot(dcs_th, hexbin_alex, S_name='Su', gridsize=60, ax=ax)
E = dcs_th.E[0]
E_std = E[E > 0].std()
ax.text(0.05, 0.95, '$th_D$ = %d, σ = %.3f\n#B = %d' % (th, E_std, dcs_th.num_bursts), va='top', ha='left',
transform=ax.transAxes, fontsize=12)
ax.set_title('')
plt.setp(AX[1:], ylabel='')
AX[0].set_ylabel('$S_u$')
plt.subplots_adjust(hspace=0, wspace=0)
plt.text(0.5, 0.89, 'DCBS', transform=fig.transFigure, fontsize=16);
Th = (20, 40, 60, 80, 160)
fig, AX = plt.subplots(1, len(Th), figsize=(3.5*len(Th), 3.5), sharex=True, sharey=True)
for irow, th in enumerate(Th):
ax = AX[irow]
dcs_th2 = d2.select_bursts(select_bursts.size, th1=th, **size_sel_kws)
dcs_th = dcs_th2.collapse(update_gamma=False, skip_ch=skip_ch)
dplot(dcs_th, hist_fret, ax=ax)
E = dcs_th.E[0]
E_std = E[E > 0].std()
ax.text(0.05, 0.95, '$th_D$ = %d, σ = %.3f\n#B = %d' % (th, E_std, dcs_th.num_bursts), va='top', ha='left',
transform=ax.transAxes, fontsize=12)
ax.set_title('')
plt.setp(AX[1:], ylabel='')
plt.subplots_adjust(hspace=0, wspace=0)
plt.text(0.5, 0.89, 'DCBS', transform=fig.transFigure, fontsize=16);
dc = d.collapse(update_gamma=False, skip_ch=skip_ch)
size_Dex = 80
size_sel_kws = dict(ph_sel=Ph_sel(Dex='DAem', Aex='Dem'), gamma=0.5,
na_comp=True, naa_comp=False, naa_aexonly=False)
Math(d._burst_sizes_pax_formula(**size_sel_kws))
size_Aex = 20
naa_sel_kws = dict(ph_sel=Ph_sel(Aex='Aem'), gamma=0.5,
na_comp=False, naa_comp=False, naa_aexonly=False)
Math(d._burst_sizes_pax_formula(**naa_sel_kws))
dcs = dc.select_bursts(select_bursts.size, th1=size_Dex, **size_sel_kws)
Su = [(nd + na)/(nd + na + naa) for nd, na, naa in zip(dcs.nd, dcs.na, dcs.naa)]
dcs.add(Su=Su)
Math(d._burst_sizes_pax_formula(**size_sel_kws))
dcsf = dcs.select_bursts(select_bursts.size, th1=size_Aex, **naa_sel_kws)
Su = [(nd + na)/(nd + na + naa) for nd, na, naa in zip(dcsf.nd, dcsf.na, dcsf.naa)]
dcsf.add(Su=Su)
Math(d._burst_sizes_pax_formula(**naa_sel_kws))
g = alex_jointplot(dcs, S_name='Su', vmax=240, gridsize=80)
g.ax_joint.set_xlim(-0.19, 1);
g.ax_joint.set_ylim(0., 1.2);
g.ax_joint.set_ylabel('$S_u$');
g = alex_jointplot(dcsf, S_name='Su', vmax=240, gridsize=80)
g.ax_joint.set_xlim(-0.19, 1);
g.ax_joint.set_ylim(0., 1.2);
g.ax_joint.set_ylabel('$S_u$');
size_Dex = 40
size_sel_kws = dict(ph_sel=Ph_sel(Dex='DAem'), gamma=0.5,
na_comp=False, naa_comp=False, naa_aexonly=False)
Math(d._burst_sizes_pax_formula(**size_sel_kws))
size_Aex = 20
naa_sel_kws = dict(ph_sel=Ph_sel(Aex='Aem'), gamma=0.5,
na_comp=False, naa_comp=False, naa_aexonly=False)
Math(d._burst_sizes_pax_formula(**naa_sel_kws))
dcs = dc.select_bursts(select_bursts.size, th1=size_Dex, **size_sel_kws)
Su = [(nd + na)/(nd + na + naa) for nd, na, naa in zip(dcs.nd, dcs.na, dcs.naa)]
dcs.add(Su=Su)
d._aex_dex_ratio
1.0
dcsf = dcs.select_bursts(select_bursts.size, th1=size_Aex, **naa_sel_kws)
Su = [(nd + na)/(nd + na + naa) for nd, na, naa in zip(dcsf.nd, dcsf.na, dcsf.naa)]
E1 = [ na/(nd + na) for nd, na in zip(dcsf.nd, dcsf.na)]
Epax = [(2 * na)/(nd + nda + 2 * na) for nd, na, nda in zip(dcsf.nd, dcsf.na, dcsf.nda)]
dcsf.add(Su=Su, E1=E1)
E = dcsf.E[0]
E.std(), E1[0].std()
(0.074804367127414645, 0.07520919981327126)
g = alex_jointplot(dcsf, S_name='Su', vmax=240, gridsize=80)
g.ax_joint.set_xlim(-0.19, 1);
g.ax_joint.set_ylim(0., 1.2);
g.ax_joint.set_ylabel('$S_u$');
g = alex_jointplot(dcsf, E_name='E1', S_name='Su', vmax=240, gridsize=80)
g.ax_joint.set_xlim(-0.19, 1);
g.ax_joint.set_ylim(0., 1.2);
g.ax_joint.set_xlabel('$E_{pax}$');
g.ax_joint.set_ylabel('$S_u$');
g = alex_jointplot(dcs, S_name='Su', vmax=240, gridsize=80)
g.ax_joint.set_xlim(-0.19, 1);
g.ax_joint.set_ylim(0., 1.2);
g.ax_joint.set_ylabel('$S_u$');
g = alex_jointplot(dcsf, S_name='Su', vmax=240, gridsize=80)
g.ax_joint.set_xlim(-0.19, 1);
g.ax_joint.set_ylim(0., 1.2);
g.ax_joint.set_ylabel('$S_u$');
dcs15 = dcs.select_bursts(select_bursts.time, time_s1=0, time_s2=15)
Su = [(nd + na)/(nd + na + naa) for nd, na, naa in zip(dcs15.nd, dcs15.na, dcs15.naa)]
dcs15.add(Su=Su)
g = alex_jointplot(dcs15, S_name='Su')
g.ax_joint.set_ylabel('$S_u$');
dcs5 = dcs.select_bursts(select_bursts.time, time_s1=0, time_s2=5)
Su = [(nd + na)/(nd + na + naa) for nd, na, naa in zip(dcs5.nd, dcs5.na, dcs5.naa)]
dcs5.add(Su=Su)
g = alex_jointplot(dcs5, S_name='Su')
g.ax_joint.set_ylabel('$S_u$');