Example usage of neuralib.imaging in cellular imaging
Example of Scanbox acqusition system and suite2p registration
Example of calculating dF/F
Example of converting dF/F to spike (using either
OASISdeconvolution orCascademodel evaluation)
[1]:
import logging
import absl
import matplotlib.pyplot as plt
import numpy as np
from neuralib.imaging.scanbox.viewer import SBXViewer
from neuralib.imaging.spikes import oasis_dcnv, cascade_predict
from neuralib.imaging.suite2p import *
from neuralib.io.dataset import load_example_scanbox, load_example_suite2p
[2]:
%load_ext autoreload
%autoreload
print the Scanbox information as dictionary
[3]:
# replace to your own scanbox .mat output file
# sbx = SBXInfo.load(...)
sbx = load_example_scanbox()
sbx.print_asdict()
{'abort_bit': 0,
'area_line': 1,
'ballmotion': array([], dtype=uint8),
'bytesPerBuffer': 9504000,
'calibration': [{'delta': array([ 99.9927, -99.9927]),
'gain_resonant_mult': 1,
'uv': array([50, 57], dtype=uint8),
'x': 1.7542578947368421,
'y': -1.999854},
{'delta': array([ 99.9927, -99.9927]),
'gain_resonant_mult': 1,
'uv': array([54, 62], dtype=uint8),
'x': 1.6127854838709677,
'y': -1.8517166666666667},
{'delta': array([ 79.6727, -79.6727]),
'gain_resonant_mult': 1,
'uv': array([55, 62], dtype=uint8),
'x': 1.2850435483870968,
'y': -1.4485945454545455},
{'delta': array([ 69.5127, -69.5127]),
'gain_resonant_mult': 1,
'uv': array([55, 62], dtype=uint8),
'x': 1.1211725806451611,
'y': -1.2638672727272726},
{'delta': array([ 59.3725, -59.3725]),
'gain_resonant_mult': 1,
'uv': array([59, 66], dtype=uint8),
'x': 0.8995833333333334,
'y': -1.006313559322034},
{'delta': array([ 49.9864, -49.9864]),
'gain_resonant_mult': 1,
'uv': array([57, 63], dtype=uint8),
'x': 0.7934349206349207,
'y': -0.8769543859649124},
{'delta': array([ 39.8264, -39.8264]),
'gain_resonant_mult': 1,
'uv': array([55, 62], dtype=uint8),
'x': 0.6423612903225806,
'y': -0.7241163636363637},
{'delta': array([ 34.3694, -34.3694]),
'gain_resonant_mult': 1,
'uv': array([56, 63], dtype=uint8),
'x': 0.5455460317460318,
'y': -0.6137392857142857},
{'delta': array([ 29.6863, -29.6863]),
'gain_resonant_mult': 1,
'uv': array([58, 66], dtype=uint8),
'x': 0.4497924242424242,
'y': -0.5118327586206897},
{'delta': array([ 24.9833, -24.9833]),
'gain_resonant_mult': 1,
'uv': array([59, 65], dtype=uint8),
'x': 0.38435846153846154,
'y': -0.4234457627118644},
{'delta': array([ 19.5263, -19.5263]),
'gain_resonant_mult': 1,
'uv': array([54, 58], dtype=uint8),
'x': 0.33666034482758617,
'y': -0.3615981481481481},
{'delta': array([ 19.5263, -19.5263]),
'gain_resonant_mult': 1,
'uv': array([65, 71], dtype=uint8),
'x': 0.27501830985915493,
'y': -0.3004046153846154},
{'delta': array([ 14.8431, -14.8431]),
'gain_resonant_mult': 1,
'uv': array([59, 62], dtype=uint8),
'x': 0.2394048387096774,
'y': -0.2515779661016949}],
'channels': 2,
'config': {'agc': {'agc_prctile': array([1.e-05, 1.e-04]),
'enable': 0,
'threshold': 250},
'coord_abs': array([ -312.975625 , 499.28859375, -11021.09375 , 0. ]),
'coord_rel': array([ -312.975625 , 499.28859375, -11021.09375 , 0. ]),
'frame_times': array([4.800000e-02, 8.100000e-02, 1.140000e-01, ..., 3.606479e+03,
3.606515e+03, 3.606549e+03]),
'frames': 108128,
'host_name': '',
'knobby': {'pos': {'a': 0, 'x': 0, 'y': 0, 'z': -300.78},
'schedule': array([[ 0, 0, 10, 0, 30],
[ 0, 0, 10, 0, 60],
[ 0, 0, 10, 0, 90],
[ 0, 0, 10, 0, 120],
[ 0, 0, 10, 0, 150],
[ 0, 0, 10, 0, 180]], dtype=uint8)},
'laser_power': 204.00000000000003,
'laser_power_perc': ' 80%',
'lines': 528,
'magnification': 4,
'magnification_list': array(['1.0', '1.2', '1.4', '1.7', '2.0', '2.4', '2.8', '3.4', '4.0',
'4.8', '5.7', '6.7', '8.0'], dtype='<U3'),
'objective': {'name': 'Nikon 16x/0.8w/WD3.0'},
'objective_type': 1,
'pmt0_gain': 0.68,
'pmt1_gain': 0.55,
'wavelength': 920},
'messages': array([], dtype=object),
'nchan': None,
'objective': 'Nikon 16x/0.8w/WD3.0',
'opto2pow': array([], dtype=uint8),
'otparam': array([], dtype=uint8),
'otwave': array([], dtype=uint8),
'otwave_um': array([], dtype=uint8),
'otwavestyle': 1,
'postTriggerSamples': 9000,
'power_depth_link': 0,
'recordsPerBuffer': 264,
'resfreq': 7915,
'scanbox_version': '2',
'scanmode': 0,
'sz': array([528, 796], dtype=uint16),
'usernotes': array([], dtype='<U1'),
'volscan': 0}
View the .sbx file
[ ]:
directory = ... # directory contains the .sbx and .mat output from scanbox
sbx_viewer = SBXViewer(directory)
# play 100 to 200 frames
sbx_viewer.play(slice(100, 200), plane=0, channel=0)
# save as tiff sequences
sbx_viewer.to_tiff(slice(100, 200), plane=0, channel=0, output='test.tiff')
Load Suite2p results
[4]:
# s2p = Suite2PResult.load(..., cell_prob=0.0, channel=0) # replace ... to suite2p base directory (*/suite2p/plane*)
s2p = load_example_suite2p()
See the Suite2p registered ROIs
[5]:
# see ROIs
_, ax = plt.subplots(1, 2)
pix = get_soma_pixel(s2p, color_diff=True)
pix[pix == 0] = np.nan
ax[0].imshow(s2p.image_mean, cmap='grey')
ax[1].imshow(pix, cmap='nipy_spectral')
plt.show()
Calculate dF/F
[6]:
dff = get_neuron_signal(s2p, n=np.arange(10), signal_type='df_f', normalize=True, dff=True, correct_neuropil=True)[0]
# pick the first 100 sec
points = int(100 * s2p.fs)
dff = dff[:, :points]
n_neurons = dff.shape[0]
y = 0
for i in range(n_neurons):
plt.plot(dff[i] + y)
y += 1
plt.xlabel('n_points')
plt.ylabel('# neuron ID')
plt.show()
Calculate deconvolution spikes using OASIS method
[7]:
spks = oasis_dcnv(dff, tau=s2p.indicator_tau, fs=s2p.fs)
spks = spks / np.max(spks, axis=1, keepdims=True)
n_neurons = spks.shape[0]
y = 0
for i in range(n_neurons):
plt.plot(spks[i] + y)
y += 1
plt.xlabel('n_points')
plt.ylabel('# neuron ID')
plt.show()
Predict spikes from dF/F using pre-trained model from Cascade
[8]:
# verbose disable
absl.logging.set_verbosity('error')
logging.getLogger('tensorflow').setLevel(logging.ERROR)
spks = cascade_predict(dff, model_type='Global_EXC_30Hz_smoothing100ms', verbose=False)
spks = spks / np.max(spks, axis=1, keepdims=True)
n_neurons = spks.shape[0]
y = 0
for i in range(n_neurons):
plt.plot(spks[i] + y)
y += 1
plt.xlabel('n_points')
plt.ylabel('# neuron ID')
plt.show()
[INFO][25-03-13 11:23:03] - The selected model was trained on 18 datasets, with 5 ensembles for each noise level, at a sampling rate of 30 Hz, with a resampled ground truth that was smoothed with a Gaussian kernelof a standard deviation of 200 ms.
[INFO][25-03-13 11:23:03] - Loaded model was trained at frame rate 30 Hz
[INFO][25-03-13 11:23:03] - Given argument traces contains 10 neurons and 2998 frames.
[INFO][25-03-13 11:23:03] - Noise levels (mean, std; in standard units): 0.46, 0.14
Predictions for noise level 1
... ensemble 0
4/4 ━━━━━━━━━━━━━━━━━━━━ 0s 26ms/step
/opt/homebrew/Caskroom/miniconda/base/envs/py310/lib/python3.10/site-packages/keras/src/optimizers/base_optimizer.py:33: UserWarning: Argument `decay` is no longer supported and will be ignored.
warnings.warn(
... ensemble 1
4/4 ━━━━━━━━━━━━━━━━━━━━ 0s 28ms/step
... ensemble 2
4/4 ━━━━━━━━━━━━━━━━━━━━ 0s 27ms/step
... ensemble 3
4/4 ━━━━━━━━━━━━━━━━━━━━ 0s 25ms/step
... ensemble 4
4/4 ━━━━━━━━━━━━━━━━━━━━ 0s 27ms/step
Predictions for noise level 2
No neurons for this noise level: 2
Predictions for noise level 3
No neurons for this noise level: 3
Predictions for noise level 4
No neurons for this noise level: 4
Predictions for noise level 5
No neurons for this noise level: 5
Predictions for noise level 6
No neurons for this noise level: 6
Predictions for noise level 7
No neurons for this noise level: 7
Predictions for noise level 8
No neurons for this noise level: 8
