-
Notifications
You must be signed in to change notification settings - Fork 0
/
2_cropping.py
143 lines (130 loc) · 7.07 KB
/
2_cropping.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
######################################
animal = '999'
alert_when_done = True
######################################
import os
import shlex
import subprocess
import time
import numpy as np
import pickle5 as pkl
from natsort import natsorted
from scipy.signal import butter, sosfiltfilt
from scipy.signal import hilbert
from utils import alert
butter_order = 3
theta_band = [4, 12]
spike_band = [300, 3000]
sampling_rate = 20000
frame_rate = 50 # frame rate of the video
start_time = time.time()
# folder preparation
animal_folder = r'E:/anxiety_ephys/' + animal + '/'
sub_folder = 'circus'
target_folder = animal_folder + sub_folder + '/'
experiment_names = natsorted(os.listdir(animal_folder))
if 'circus' in experiment_names:
experiment_names.remove('circus')
circus_entrypoint = target_folder + 'dat_files/' + experiment_names[
0] + '_0.dat' # file to run spyking circus on (first of the data files)
if os.listdir(target_folder + 'dat_files_mod'): #can only crop once
raise Exception('Croppig already done, to redo run config.py again first.')
mouse_is_late = {'2021-02-19_mBWfus010_EZM_ephys': 70,
# seconds delay from trigger signal to mouse in the maze #same as in 1_config.py
'2021-02-19_mBWfus009_EZM_ephys': 42,
'2021-02-26_mBWfus012_EZM_ephys': 35}
accomodation = 20 # time to let the late mice settle #same as in 1_config.py
logbook = np.zeros(len(experiment_names)) # will contain the length of the recordings for all sessions
for index, experiment_name in enumerate(experiment_names):
eventfile = animal_folder + experiment_name + '/ephys_processed/' + experiment_name + '_events.pkl'
##get events/movement
with open(eventfile, 'rb') as f:
events = pkl.load(f)
movement = events['movement']
xy = np.array([movement['calib_traj_x'], movement['calib_traj_y']],
dtype=np.float32) # first row: x coord, second row: y coord
movement_trigger_file = animal_folder + '/' + experiment_name + '/log.txt'
with open(movement_trigger_file, 'r') as file:
data = file.read().replace('\n', ' ')
point = data.index('.')
movement_trigger = int(float(data[point - 2:point + 3]) * frame_rate) # get the 50Hz movement trigger
# introduce offset for the movement #is done in 1_config.py for ephys
if experiment_name in mouse_is_late:
off = (mouse_is_late[experiment_name] + accomodation) * frame_rate
else:
off = 0
movement_trigger += off
xy = xy[:, movement_trigger:] # crop from movement trigger
mPFC_concatenated = np.load(target_folder + 'mPFC_raw/' + experiment_name + '.npy')
vHIP_concatenated = np.load(target_folder + 'vHIP_raw/' + experiment_name + '.npy')
##cut out the time segments specified in the cutter file
if os.path.exists(target_folder + 'cutter/' + experiment_name + '.npy'):
mPFC_spike_range = np.load(target_folder + 'mPFC_spike_range/' + experiment_name)
cutter = np.load(
target_folder + 'cutter/' + experiment_name + '.npy') # first row start time of cuts, second row stop time, dtype uint32
boolean_sampling_rate = np.ones(mPFC_concatenated.shape[1], dtype=np.bool)
boolean_frame_rate = np.ones(xy.shape[1], dtype=np.bool)
for cut in range(cutter.shape[1]):
boolean_sampling_rate[cutter[0, cut] * sampling_rate // frame_rate: cutter[
1, cut] * sampling_rate // frame_rate] = 0 # mask the values to cut out 20000Hz
boolean_frame_rate[cutter[0, cut]: cutter[1, cut]] = 0 # mask the values to cut out 20000Hz
##discard the masked values:
mPFC_concatenated = mPFC_concatenated[:, boolean_sampling_rate]
vHIP_concatenated = vHIP_concatenated[:, boolean_sampling_rate]
mPFC_spike_range = mPFC_spike_range[:, boolean_sampling_rate]
xy = xy[:, boolean_frame_rate]
np.save(target_folder + 'mPFC_spike_range/' + experiment_name,
mPFC_spike_range.astype(np.int16)) # used to calculate mean waveform in 4_sanity_check.py
if experiment_name[-9:-6] == 'EZM':
transitions = {}
for mode in events['transitions']:
event_indices = events['transitions'][mode] - off
event_boolean = np.zeros(mPFC_concatenated.shape[1], dtype=bool)
event_boolean[event_indices] = 1
event_boolean = event_boolean[boolean_frame_rate]
transitions[mode] = list(np.nonzero(event_boolean))
else:
if experiment_name[-9:-6] == 'EZM':
transitions = {}
for mode in events['transitions']:
transitions[mode] = events['transitions'][mode] - off
if experiment_name[-9:-6] == 'EZM':
with open(target_folder + 'transition_files/' + experiment_name + '.pkl', 'wb') as f:
pkl.dump(transitions, f)
np.save(target_folder + 'movement_files/' + experiment_name, xy.astype(np.float32))
##theta bandpass filtering:
sos_theta = butter(N=butter_order, Wn=theta_band, btype='bandpass', analog=False, output='sos', fs=sampling_rate)
theta_filtered = sosfiltfilt(sos_theta, vHIP_concatenated, axis=1)
hilbert_phase = np.angle(hilbert(theta_filtered, axis=1), deg=True) # get theta phase
data_for_spikesorting = np.transpose(mPFC_concatenated)
# save files:
data_for_spikesorting.tofile(
target_folder + 'dat_files/' + experiment_names[index] + '_' + str(index) + '.dat') # used by spykingcircus
np.save(target_folder + 'vHIP_phase/' + experiment_name, hilbert_phase.astype(np.int16)) # save theta phase
logbook[index] = data_for_spikesorting.shape[0] # logbook[i] contains the length of session i, 20000Hz
np.save(target_folder + 'utils/logbook', logbook.astype(np.uint32))
before_clustering_time = time.time()
##start clustering process:
cluster_command = 'spyking-circus ' + circus_entrypoint + ' -c 10' #-c param: number of CPU cores to be used for processing
os.system(cluster_command)
##convert output of spykingcircus for the phy viewer:
converter_command = 'spyking-circus ' + circus_entrypoint + ' -m converting -c 10'
args = shlex.split(converter_command)
converter = subprocess.run(args, stdout=subprocess.PIPE,
input='a\n', encoding='ascii')
print('Converter return_code: {}'.format(converter.returncode))
end_time = time.time()
print('Cropping for animal {} done! \nTime needed: {} minutes. Time for SpykingCircus: {}'.format(animal, (
end_time - start_time) / 60, (end_time - before_clustering_time) / 60))
if alert_when_done:
alert()
##create a batchfile to start the phy viewer by clicking on the file
viewer_command = '@echo off \ncall conda activate circus \ncircus-gui-python ' + circus_entrypoint
with open(target_folder + 'utils/start_viewer.bat', 'w') as f:
f.write(viewer_command)
##start the phy viewer:
viewer_command = 'circus-gui-python ' + circus_entrypoint
args = shlex.split(viewer_command)
viewer = subprocess.run(args, stdout=subprocess.PIPE,
encoding='ascii')
print('Viewer return_code: {}'.format(viewer.returncode))