-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathpipeline_holo.py
1298 lines (1094 loc) · 58 KB
/
pipeline_holo.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
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
from __future__ import division
from __future__ import print_function
from macpath import basename
#from matplotlib.tests.test_backend_pgf import test_bbox_inches
#from __main__ import traceback
# Pipeline to obtain rois from 2p data based on Caiman
__author__ = 'Nuria'
from numpy.distutils.system_info import dfftw_info
try:
zip, str, map, range
except:
from builtins import zip
from builtins import str
from builtins import map
from builtins import range
from past.utils import old_div
from skimage import io
import tifffile
import glob
import matplotlib.pyplot as plt
import numpy as np
import os
import sys
import time
import scipy
import h5py
import pandas as pd
import cv2
try:
cv2.setNumThreads(0)
except:
pass
try:
if __IPYTHON__:
print((1))
# this is used for debugging purposes only. allows to reload classes
# when changed
get_ipython().magic('load_ext autoreload')
get_ipython().magic('autoreload 2')
except NameError:
print('Not IPYTHON')
pass
import caiman as cm
from caiman.utils.utils import download_demo
from caiman.utils.visualization import plot_contours, nb_view_patches, nb_plot_contour
from caiman.source_extraction.cnmf import cnmf as cnmf
from caiman.motion_correction import MotionCorrect
from caiman.source_extraction.cnmf.utilities import detrend_df_f
from caiman.components_evaluation import estimate_components_quality_auto
from caiman.components_evaluation import evaluate_components_CNN
from caiman.motion_correction import motion_correct_iteration
import bokeh.plotting as bpl
from skimage.feature import peak_local_max
from scipy.stats.mstats import zscore
from scipy import ndimage
import copy
from matplotlib import interactive
import sys, traceback
import imp
interactive(True)
def all_run(folder, animal, day, number_planes=4, number_planes_total=6):
"""
Function to run all the different functions of the pipeline that gives back the analyzed data
Folder (str): folder where the input/output is/will be stored
animal/day (str) to be analyzed
number_planes (int): number of planes that carry information
number_planes_total (int): number of planes given back by the recording system, it may differ from number_planes
to provide time for the objective to return to origen"""
folder_path = folder + 'raw/' + animal + '/' + day + '/'
folder_final = folder + 'processed/' + animal + '/' + day + '/'
err_file = open(folder_path + "errlog.txt", 'a+') # ERROR HANDLING
if not os.path.exists(folder_final):
os.makedirs(folder_final)
finfo = folder_path + 'wmat.mat' #file name of the mat
matinfo = scipy.io.loadmat(finfo)
ffull = [folder_path + matinfo['fname'][0]] # filename to be processed
fbase = [folder_path + matinfo['fbase'][0]]
try:
num_files, len_bmi = separate_planes(folder, animal, day, ffull, 'bmi', number_planes, number_planes_total)
num_files_b, len_base = separate_planes(folder, animal, day, fbase, 'baseline', number_planes, number_planes_total)
except Exception as e:
tb = sys.exc_info()[2]
err_file.write("\n{}\n".format(folder_path))
err_file.write("{}\n".format(str(e.args)))
traceback.print_tb(tb, file=err_file)
err_file.close()
sys.exit('Error in separate planes')
nam = folder_path + 'readme.txt'
readme = open(nam, 'w+')
readme.write("num_files_b = " + str(num_files_b) + '; \n')
readme.write("num_files = " + str(num_files)+ '; \n')
readme.write("len_base = " + str(len_base)+ '; \n')
readme.write("len_bmi = " + str(len_bmi)+ '; \n')
readme.close()
try:
analyze_raw_planes(folder, animal, day, num_files, num_files_b, number_planes, False)
except Exception as e:
tb = sys.exc_info()[2]
err_file.write("\n{}\n".format(folder_path))
err_file.write("{}\n".format(str(e.args)))
traceback.print_tb(tb, file=err_file)
err_file.close()
sys.exit('Error in analyze raw')
try:
shutil.rmtree(folder + 'raw/' + animal + '/' + day + '/separated/')
except OSError as e:
print ("Error: %s - %s." % (e.filename, e.strerror))
# This section is done nowadays in another computer to minimize time for analysis
# try:
# put_together(folder, animal, day, len_base, len_bmi, number_planes, number_planes_total)
# except Exception as e:
# tb = sys.exc_info()[2]
# err_file.write("\n{}\n".format(folder_path))
# err_file.write("{}\n".format(str(e.args)))
# traceback.print_tb(tb, file=err_file)
# err_file.close()
# sys.exit('Error in put together')
err_file.close()
def separate_planes(folder, animal, day, ffull, var='bmi', number_planes=4, number_planes_total=6, order='F', lim_bf=9000):
"""
Function to separate the different planes in the bigtiff file given by the recording system.
TO BECOME OBSOLETE
Folder (str): folder where the input/output is/will be stored
animal/day (str) to be analyzed
ffull (str): address of the file where the bigtiff is stored
var(str): variable to specify if performing BMI or baseline
number_planes(int): number of planes that carry information
number_planes_total(int): number of planes given back by the recording system, it may differ from number_planes
to provide time for the objective to return to origen
order(str): order to stablish the memmap C/F
lim_bf (int): limit of frames per split to avoid saving big tiff files"""
# function to separate a layered TIFF (tif generated with different layers info)
folder_path = folder + 'raw/' + animal + '/' + day + '/separated/'
fanal = folder + 'raw/' + animal + '/' + day + '/analysis/'
if not os.path.exists(folder_path):
os.makedirs(folder_path)
if not os.path.exists(fanal):
os.makedirs(fanal)
err_file = open("errlog.txt", 'a+') # ERROR HANDLING
print('loading image...')
ims = tifffile.TiffFile(ffull[0])
len_im = int(len(ims.pages)/number_planes_total)
dims = [len_im] + list(ims.pages[0].shape)
print('Image loaded')
num_files = int(np.ceil(len_im/lim_bf))
for plane in np.arange(number_planes):
len_im = int(len(ims.pages)/number_planes_total)
print ('length of tiff is: ' + str(len_im) + ' volumes')
for nf in np.arange(num_files):
# create the mmap file
fnamemm = folder_path + 'temp_plane_' + str(plane) + '_nf_' + str(nf) + '.mmap'
if len_im>lim_bf:
auxlen = lim_bf
len_im -= lim_bf
else:
auxlen = len_im
big_file = np.memmap(fnamemm, mode='w+', dtype=np.int16, shape=(np.prod(dims[1:]), auxlen), order=order)
# fill the mmap file
for ind in np.arange(auxlen):
new_img = ims.pages[int((ind + lim_bf*nf)*number_planes_total + plane)].asarray()
big_file[:, ind] = np.reshape(new_img, np.prod(dims[1:]), order=order)
#to plot the image before closing big_file (as a checkup that everything went smoothly)
if not os.path.exists(fanal + str(plane) + '/'):
os.makedirs(fanal + str(plane) + '/')
imgtosave = np.transpose(np.reshape(np.nanmean(big_file,1), [dims[1],dims[2]])) # careful here if the big_file is too big to average
plt.imshow(imgtosave)
plt.savefig(fanal + str(plane) + '/' + 'nf' + str(nf) + '_rawmean.png', bbox_inches="tight")
plt.close()
big_file.flush()
del big_file
# clean memory
ims.close()
# save the mmaps as tiff-files for caiman
for plane in np.arange(number_planes):
len_im = int(len(ims.pages)/number_planes_total)
print ('saving a tiff of: ' + str(len_im) + ' volumes')
for nf in np.arange(num_files):
fnamemm = folder_path + 'temp_plane_' + str(plane) + '_nf_' + str(nf) + '.mmap'
fnametiff = folder_path + var + '_plane_' + str(plane) + '_nf_' + str(nf) + '.tiff'
if len_im>lim_bf:
auxlen = lim_bf
len_im -= lim_bf
else:
auxlen = len_im
big_file = np.memmap(fnamemm, mode='r', dtype=np.int16, shape=(np.prod(dims[1:]), auxlen), order=order)
img_tosave = np.transpose(np.reshape(big_file, [dims[1], dims[2], int(auxlen)]))
io.imsave(fnametiff, img_tosave, plugin='tifffile') # saves each plane different tiff
del big_file
del img_tosave
try:
os.remove(fnamemm)
except OSError as e: ## if failed, report it back to the user ##
print ("Error: %s - %s." % (e.filename, e.strerror))
len_im = int(len(ims.pages)/number_planes_total)
del ims
err_file.close()
return num_files, len_im
def separate_planes_multiple_baseline(folder, animal, day, fbase1, fbase2, var='baseline', number_planes=4, number_planes_total=6, order='F', lim_bf=9000):
"""
Function to separate the different planes in the bigtiff file given by the recording system WHEN there is more than one baseline file.
TO BECOME OBSOLETE
Folder (str): folder where the input/output is/will be stored
animal/day (str) to be analyzed
ffull (str): address of the file where the bigtiff is stored
ffull2 (str): address of the file where the consecutive bigtiff is stored
var(str): variable to specify if performing BMI or baseline
number_planes(int): number of planes that carry information
number_planes_total(int): number of planes given back by the recording system, it may differ from number_planes
to provide time for the objective to return to origen
order(str): order to stablish the memmap C/F
lim_bf (int): limit of frames per split to avoid saving big tiff files"""
# function to separate a layered TIFF (tif generated with different layers info)
folder_path = folder + 'raw/' + animal + '/' + day + '/separated/'
fanal = folder + 'raw/' + animal + '/' + day + '/analysis/'
if not os.path.exists(folder_path):
os.makedirs(folder_path)
if not os.path.exists(fanal):
os.makedirs(fanal)
err_file = open("errlog.txt", 'a+') # ERROR HANDLING
# load the big file we will need to separate
print('loading image...')
imb1 = tifffile.TiffFile(fbase1[0])
imb2 = tifffile.TiffFile(fbase2[0])
print('Images loaded')
len_imb1 = int(len(imb1.pages)/number_planes_total)
len_imb2 = int(len(imb2.pages)/number_planes_total)
dims1 = [len_imb1] + list(imb1.pages[0].shape)
dims2 = [len_imb2] + list(imb2.pages[0].shape)
dims = dims1 + np.asarray([dims2[0], 0, 0])
len_im = copy.deepcopy(int(dims[0]))
num_files = int(np.ceil(len_im/lim_bf))
for plane in np.arange(number_planes):
first_images_left = int(dims1[0])
len_im = int(dims[0])
print ('length of tiff is: ' + str(len_im) + ' volumes')
for nf in np.arange(num_files):
# create the mmap file
fnamemm = folder_path + 'temp_plane_' + str(plane) + '_nf_' + str(nf) + '.mmap'
if len_im>lim_bf:
auxlen = lim_bf
len_im -= lim_bf
else:
auxlen = len_im
big_file = np.memmap(fnamemm, mode='w+', dtype=np.int16, shape=(np.prod(dims[1:]), auxlen), order=order)
# fill the mmap file
for ind in np.arange(auxlen):
if first_images_left > 0:
new_img = imb1.pages[int(plane+(ind + lim_bf*nf)*number_planes_total)].asarray()
first_images_left -= 1
else:
new_img = imb2.pages[int(plane+(ind - int(dims1[0]) + lim_bf*nf)*number_planes_total)].asarray()
big_file[:, ind] = np.reshape(new_img, np.prod(dims[1:]), order=order)
#to plot the image before closing big_file (as a checkup that everything went smoothly)
if not os.path.exists(fanal + str(plane) + '/'):
os.makedirs(fanal + str(plane) + '/')
imgtosave = np.transpose(np.reshape(np.nanmean(big_file,1), [dims[1],dims[2]]))
plt.imshow(imgtosave)
plt.savefig(fanal + str(plane) + '/' + 'nf' + str(nf) + '_rawmean.png', bbox_inches="tight")
plt.close()
big_file.flush()
del big_file
# save the mmaps as tiff-files for caiman
for plane in np.arange(number_planes):
len_im = int(dims[0])
print ('saving a tiff of: ' + str(len_im) + ' volumes')
for nf in np.arange(num_files):
fnamemm = folder_path + 'temp_plane_' + str(plane) + '_nf_' + str(nf) + '.mmap'
fnametiff = folder_path + var + '_plane_' + str(plane) + '_nf_' + str(nf) + '.tiff'
if len_im>lim_bf:
auxlen = lim_bf
len_im -= lim_bf
else:
auxlen = len_im
big_file = np.memmap(fnamemm, mode='r', dtype=np.int16, shape=(np.prod(dims[1:]), auxlen), order=order)
img_tosave = np.transpose(np.reshape(big_file, [dims[1], dims[2], int(auxlen)]))
io.imsave(fnametiff, img_tosave, plugin='tifffile') # saves each plane different tiff
del big_file
del img_tosave
try:
os.remove(fnamemm)
except OSError as e: ## if failed, report it back to the user ##
print ("Error: %s - %s." % (e.filename, e.strerror))
# clean memory
imb1.close()
imb2.close()
del imb1
del imb2
return num_files, int(dims[0])
def analyze_raw_planes(folder, animal, day, num_files, num_files_b, number_planes=4, dend=False, display_images=True):
"""
Function to analyze every plane and get the result in a hdf5 file. It uses caiman_main
Folder(str): folder where the input/output is/will be stored
animal/day(str) to be analyzed
num_files(int): number of files for the bmi file
num_files_b(int): number of files for the baseline file
number_planes(int): number of planes that carry information
number_planes_total(int): number of planes given back by the recording system, it may differ from number_planes
to provide time for the objective to return to origen
dend(bool): Boleean to change parameters to look for neurons or dendrites
display_images(bool): to display and save different plots"""
folder_path = folder + 'raw/' + animal + '/' + day + '/separated/'
finfo = folder + 'raw/' + animal + '/' + day + '/wmat.mat' #file name of the mat
matinfo = scipy.io.loadmat(finfo)
initialZ = int(matinfo['initialZ'][0][0])
fr = matinfo['fr'][0][0]
if dend:
sec_var = 'Dend'
else:
sec_var = ''
print('*************Starting with analysis*************')
neuron_mats = []
for plane in np.arange(number_planes):
dff_all = []
neuron_act_all = []
fnames = []
for nf in np.arange(int(num_files_b)):
fnames.append(folder_path + 'baseline' + '_plane_' + str(plane) + '_nf_' + str(nf) + '.tiff')
print('performing plane: ' + str(plane))
for nf in np.arange(int(num_files)):
fnames.append(folder_path + 'bmi' + '_plane_' + str(plane) + '_nf_' + str(nf) + '.tiff')
fpath = folder + 'raw/' + animal + '/' + day + '/analysis/' + str(plane) + '/'
if not os.path.exists(fpath):
os.makedirs(fpath)
try:
f = h5py.File(folder + 'raw/' + animal + '/' + day + '/' + 'bmi_' + sec_var + '_' + str(plane) + '.hdf5', 'w-')
except IOError:
print(" OOPS!: The file already existed ease try with another file, new results will NOT be saved")
continue
zval = calculate_zvalues(folder, plane)
print(fnames)
dff, com, cnm2, totdes, SNR = caiman_main(fpath, fr, fnames, zval, dend, display_images)
print ('Caiman done: saving ... plane: ' + str(plane) + ' file: ' + str(nf))
Asparse = scipy.sparse.csr_matrix(cnm2.estimates.A)
f.create_dataset('dff', data = dff) #activity
f.create_dataset('com', data = com) #distance
g = f.create_group('Nsparse') #neuron shape
g.create_dataset('data', data = Asparse.data)
g.create_dataset('indptr', data = Asparse.indptr)
g.create_dataset('indices', data = Asparse.indices)
g.attrs['shape'] = Asparse.shape
f.create_dataset('neuron_act', data = cnm2.estimates.S) #spikes
f.create_dataset('C', data = cnm2.estimates.C) #temporal activity
f.create_dataset('base_im', data = cnm2.estimates.b) #baseline image
f.create_dataset('tot_des', data = totdes) #total displacement during motion correction
f.create_dataset('SNR', data = SNR) #SNR of neurons
f.close()
print('... done')
def put_together(folder, animal, day, number_planes=4, number_planes_total=6, sec_var='', toplot=True, trial_time=30, tocut=False, len_experiment=30000):
"""
Function to put together the different hdf5 files obtain for each plane and convey all the information in one and only hdf5
it requires somo files in the original folder
Folder(str): folder where the input/output is/will be stored
animal/day(str) to be analyzed
len_base(int): length of the baseline file (in frames)
len_bmi(int): length of the bmi file (in frames)
number_planes(int): number of planes that carry information
number_planes_total(int): number of planes given back by the recording system, it may differ from number_planes
to provide time for the objective to return to origen
sec_var(str): secondary variable to save file. For extra information
toplot(bool): to allow plotting/saving of some results"""
# Folder to load/save
folder_path = folder + 'raw/' + animal + '/' + day + '/'
folder_dest = folder + 'processed/' + animal + '/'
fanal = folder_path + 'analysis/'
if not os.path.exists(folder_dest):
os.makedirs(folder_dest)
if not os.path.exists(fanal):
os.makedirs(fanal)
# Load information
print ('loading info')
vars = imp.load_source('readme', folder_path + 'readme.txt')
finfo = folder_path + 'wmat.mat' #file name of the mat
matinfo = scipy.io.loadmat(finfo)
ffull = [folder_path + matinfo['fname'][0]]
metadata = tifffile.TiffFile(ffull[0]).scanimage_metadata
fr = matinfo['fr'][0][0]
folder_red = folder + 'raw/' + animal + '/' + day + '/'
fmat = folder_red + 'red.mat'
redinfo = scipy.io.loadmat(fmat)
red = redinfo['red'][0]
com_list = []
neuron_plane = np.zeros(number_planes)
tot_des_plane = np.zeros((number_planes, 2))
for plane in np.arange(number_planes):
try:
f = h5py.File(folder_path + 'bmi_' + sec_var + '_' + str(plane) + '.hdf5', 'r')
except OSError:
break
auxb = np.nansum(np.asarray(f['base_im']),1)
bdim = int(np.sqrt(auxb.shape[0]))
base_im = np.transpose(np.reshape(auxb, [bdim,bdim]))
fred = folder_path + 'red' + str(plane) + '.tif'
red_im = tifffile.imread(fred)
auxdff = np.asarray(f['dff'])
auxC = np.asarray(f['C'])
neuron_plane[plane] = auxC.shape[0]
tot_des_plane[plane,:] = np.asarray(f['tot_des'])
if np.nansum(auxdff) == 0:
auxdff = auxC * np.nan
if plane == 0:
all_dff = auxdff
all_C = np.asarray(f['C'])
all_com = np.asarray(f['com'])
com_list.append(np.asarray(f['com']))
g = f['Nsparse']
all_neuron_shape = scipy.sparse.csr_matrix((g['data'][:], g['indices'][:], g['indptr'][:]), g.attrs['shape'])
all_base_im = np.ones((base_im.shape[0], base_im.shape[1], number_planes)) *np.nan
all_neuron_act = np.asarray(f['neuron_act'])
all_red_im = np.ones((red_im.shape[0], red_im.shape[1], number_planes)) *np.nan
all_SNR = np.asarray(f['SNR'])
else:
all_dff = np.concatenate((all_dff, auxdff), 0)
all_C = np.concatenate((all_C, np.asarray(f['C'])), 0)
all_com = np.concatenate((all_com, np.asarray(f['com'])), 0)
com_list.append(np.asarray(f['com']))
g = f['Nsparse']
gaux = scipy.sparse.csr_matrix((g['data'][:], g['indices'][:], g['indptr'][:]), g.attrs['shape'])
all_neuron_shape = scipy.sparse.hstack([all_neuron_shape, gaux])
all_neuron_act = np.concatenate((all_neuron_act, np.asarray(f['neuron_act'])), 0)
all_SNR = np.concatenate((all_SNR, np.asarray(f['SNR'])),0)
all_red_im[:, :, plane] = red_im
all_base_im[:, :, plane] = base_im
f.close()
print ('success!!')
auxZ = np.zeros((all_com.shape))
auxZ[:,2] = np.repeat(matinfo['initialZ'][0][0],all_com.shape[0])
all_com += auxZ
# Reorganize sparse matrix of spatial components
dims = all_neuron_shape.shape
dims = [int(np.sqrt(dims[0])), int(np.sqrt(dims[0])), all_neuron_shape.shape[1]]
Asparse = scipy.sparse.csr_matrix(all_neuron_shape)
Afull = np.reshape(all_neuron_shape.toarray(),dims)
# separates "real" neurons from dendrites
print ('finding neurons')
pred, _ = evaluate_components_CNN(all_neuron_shape, dims[:2], [4,4])
nerden = np.zeros(Afull.shape[2]).astype('bool')
nerden[np.where(pred[:,1]>0.75)] = True
# obtain the real position of components A
new_com = obtain_real_com(fanal, Afull, all_com, nerden, not tocut)
# sanity check of the neuron's quality
if not tocut:
plot_Cs(fanal, all_C, nerden)
print('success!!')
# identify ens_neur (it already plots sanity check in raw/analysis
online_data = pd.read_csv(folder_path + matinfo['fcsv'][0])
try:
mask = matinfo['allmask']
except KeyError:
mask = np.nan
print('finding ensemble neurons')
ens_neur = detect_ensemble_neurons(fanal, all_C, online_data, len(online_data.keys())-2,
new_com, metadata, neuron_plane, number_planes_total, vars.len_base)
nerden[ens_neur.astype('int')] = True
# obtain trials hits and miss
trial_end, trial_start, array_t1, array_miss, hits, miss = check_trials(matinfo, vars, fr, trial_time)
if np.sum(np.isnan(trial_end)) != 0:
print ("STOPPING nan's found in the trial_end")
return
print('finding red neurons')
# obtain the neurons label as red (controlling for dendrites)
redlabel = red_channel(red, neuron_plane, nerden, Afull, new_com, all_red_im, all_base_im, fanal, number_planes, toplot=toplot)
redlabel[ens_neur.astype('int')] = True
# obtain the frequency
try:
frequency = obtainfreq(matinfo['frequency'][0], vars.len_bmi)
except KeyError:
frequency = np.nan
cursor = matinfo['cursor'][0]
if tocut:
all_C, all_dff, all_neuron_act, trial_end, trial_start, hits, miss, array_t1, array_miss, cursor, frequency = \
cut_experiment(all_C, all_dff, all_neuron_act, trial_end, trial_start, hits, miss, cursor, frequency, vars.len_base, len_experiment)
# sanity checks
if toplot:
plt.figure()
plt.plot(np.nanmean(all_C,0)/10000)
plt.title('Cs')
plt.savefig(fanal + animal + '_' + day + '_Cs.png', bbox_inches="tight")
plt.figure()
plt.plot(matinfo['cursor'][0])
plt.title('cursor')
plt.savefig(fanal + animal + '_' + day + '_cursor.png', bbox_inches="tight")
plt.close('all')
#fill the file with all the correct data!
try:
fall = h5py.File(folder_dest + 'full_' + animal + '_' + day + '_' + sec_var + '_data.hdf5', 'w-')
print('saviiiiiing')
fall.create_dataset('dff', data = all_dff) # (array) (Ft - Fo)/Fo . Increment of fluorescence
fall.create_dataset('C', data = all_C) # (array) Relative fluorescence of each component
fall.create_dataset('SNR', data = all_SNR) # (array) Signal to noise ratio of each component
fall.create_dataset('com_cm', data = all_com) # (array) Position of the components as given by caiman
fall.attrs['blen'] = vars.len_base # (int) lenght of the baseline
gall = fall.create_group('Nsparse') # (sparse matrix) spatial filter of each component
gall.create_dataset('data', data = Asparse.data) # (part of the sparse matrix)
gall.create_dataset('indptr', data = Asparse.indptr) # (part of the sparse matrix)
gall.create_dataset('indices', data = Asparse.indices) # (part of the sparse matrix)
gall.attrs['shape'] = Asparse.shape # (part of the sparse matrix)
fall.create_dataset('neuron_act', data = all_neuron_act) # (array) Spike activity (S in caiman)
fall.create_dataset('base_im', data = all_base_im) # (array) matrix with all the average image of the baseline for each plane
fall.create_dataset('red_im', data = all_red_im) # (array) matrix with all the imagesfrom the red chanel for each plane
fall.create_dataset('online_data', data = online_data) # (array) Online recordings of the BMI
fall.create_dataset('ens_neur', data = ens_neur) # (array) Index of the ensemble neurons among the rest of components
fall.create_dataset('trial_end', data = trial_end) # (array) When a trial ended. Can be a hit or a miss
fall.create_dataset('trial_start', data = trial_start) # (array) When a trial started
fall.attrs['fr'] = matinfo['fr'][0][0] # (int) Framerate
fall.create_dataset('redlabel', data = redlabel) # (array-bool) True labels neurons as red
fall.create_dataset('nerden', data = nerden) # (array-bool) True labels components as neurons
fall.create_dataset('hits', data = hits) # (array) When the animal hit the target
fall.create_dataset('miss', data = miss) # (array) When the animal miss the target
fall.create_dataset('array_t1', data = array_t1) # (array) index of the trials that ended in hit
fall.create_dataset('array_miss', data = array_miss) # (array) Index of the trials that ended in miss
fall.create_dataset('cursor', data = cursor) # (array) Online cursor of the BMI
fall.create_dataset('freq', data = frequency) # (array) Frenquency resulting of the online cursor.
fall.close()
print('all done!!')
except IOError:
print(" OOPS!: The file already existed please try with another file, no results will be saved!!!")
def check_trials(matinfo, vars, fr, trial_time=30):
trial_end = (matinfo['trialEnd'][0] + vars.len_base).astype('int')
trial_start = (matinfo['trialStart'][0] + vars.len_base).astype('int')
if len(matinfo['hits']) > 0 :
hits = (matinfo['hits'][0] + vars.len_base).astype('float')
else:
hits = []
if len(matinfo['miss']) > 0 :
miss = (matinfo['miss'][0] + vars.len_base).astype('float')
else:
miss = []
# to remove false end of trials
if trial_start[0] > trial_end[0]:
trial_end = trial_end[1:]
if trial_start.shape[0] > trial_end.shape[0]:
trial_start = trial_start[:-1]
flag_correct = False
if (trial_end.shape[0] > trial_start.shape[0]): flag_correct = True
if not flag_correct:
if (len(np.where((trial_end - trial_start)<0)[0])>0): flag_correct = True
if flag_correct:
ind = 0
while ind < trial_start.shape[0]: # CAREFUL it can get stack foreveeeeeer
tokeep = np.ones(trial_end.shape[0]).astype('bool')
if trial_end.shape[0] < trial_start.shape[0]:
trial_start = trial_start[:trial_end.shape[0]]
elif (trial_end[ind] - trial_start[ind]) < 0 :
hitloc = np.where(trial_end[ind]==hits)[0]
misloc = np.where(trial_end[ind]==miss)[0]
if len(hitloc) > 0:
hits[hitloc[0]] = np.nan
if len(misloc) > 0:
miss[misloc[0]] = np.nan
tokeep[ind]=False
trial_end = trial_end[tokeep]
else:
ind += 1
# to remove ending trials that may have occur without trial start at the end of experiment
if trial_start.shape[0] != trial_end.shape[0]:
todel = trial_end>trial_start[-1]
todel[np.where(todel)[0][0]] = False
for tend in trial_end[todel]:
hitloc = np.where(tend==hits)[0]
misloc = np.where(tend==miss)[0]
if len(hitloc) > 0:
hits[hitloc[0]] = np.nan
if len(misloc) > 0:
miss[misloc[0]] = np.nan
trial_end = trial_end[:trial_start.shape[0]]
if np.sum((trial_end - trial_start) > trial_time*fr + 10) != 0: # make sure that no trial is more than it should be +10 because it can vara bit
print ("Something wrong happened here, you better check this one out")
#return np.nan, np.nan, np.nan, np.nan
# to remove trials that ended in the same frame as they started
tokeep = np.ones(trial_start.shape[0]).astype('bool')
for ind in np.arange(trial_start.shape[0]):
if (trial_end[ind] - trial_start[ind]) == 0 :
tokeep[ind]=False
hitloc = np.where(trial_end[ind]==hits)[0]
misloc = np.where(trial_end[ind]==miss)[0]
if len(hitloc) > 0:
hits[hitloc[0]] = np.nan
if len(misloc) > 0:
miss[misloc[0]] = np.nan
hits = hits[~np.isnan(hits)]
miss = miss[~np.isnan(miss)]
trial_end = trial_end[tokeep]
trial_start = trial_start[tokeep]
elif trial_end.shape[0] < trial_start.shape[0]:
trial_start = trial_start[:trial_end.shape[0]]
# preparing the arrays (number of trial for hits/miss)
array_t1 = np.zeros(hits.shape[0], dtype=int)
array_miss = np.zeros(miss.shape[0], dtype=int)
for hh, hit in enumerate(hits): array_t1[hh] = np.where(trial_end==hit)[0][0]
for mm, mi in enumerate(miss): array_miss[mm] = np.where(trial_end==mi)[0][0]
return trial_end, trial_start, array_t1, array_miss, hits, miss
def red_channel(red, neuron_plane, nerden, Afull, new_com, all_red_im, all_base_im, fanal, number_planes=4, maxdist=4, toplot=True):
"""
Function to identify red neurons with components returned by caiman
red(array-int): mask of red neurons position for each frame
nerden(array-bool): array of bool labelling as true components identified as neurons.
neuron_plane: list number_of neurons for each plane
new_com(array): position of the neurons
all_red_im(array): matrix MxNxplanes: Image of the red channel
all_base_im(array): matrix MxNxplanes: Image of the green channel
fanal(str): folder where to store the analysis sanity check
number_planes(int): number of planes that carry information
maxdist(int): spatial tolerance to assign red label to a caiman component
returns
redlabel(array-bool): boolean vector labelling as True the components that are red neurons
"""
#function to identify red neurons
all_red = []
ind_neur = 0
if len(red) < number_planes:
number_planes = len(red)
for plane in np.arange(number_planes):
maskred = copy.deepcopy(np.transpose(red[plane]))
mm = np.sum(maskred,1)
maskred = maskred[~np.isnan(mm),:].astype('float32')
# for some reason the motion correction sometimes works one way but not the other
_, _, shift, _ = motion_correct_iteration(all_base_im[:,:,plane].astype('float32'), all_red_im[:,:,plane].astype('float32'),1)
if np.nansum(abs(np.asarray(shift))) < 20: # hopefully this motion correction worked
maskred[:,0] -= shift[1].astype('float32')
maskred[:,1] -= shift[0].astype('float32')
# creates a new image with the shifts found
M = np.float32([[1, 0, -shift[1]], [0, 1, -shift[0]]])
min_, max_ = np.min(all_red_im[:,:,plane]), np.max(all_red_im[:,:,plane])
new_img = np.clip(cv2.warpAffine(all_red_im[:,:,plane], M, (all_red_im.shape[0], all_red_im.shape[1]), flags=cv2.INTER_CUBIC, borderMode=cv2.BORDER_REFLECT), min_, max_)
else:
print ('Trying other way since shift was: ' + str(shift))
# do the motion correctio the other way arround
new_img, _, shift, _ = motion_correct_iteration(all_red_im[:,:,plane].astype('float32'), all_base_im[:,:,plane].astype('float32'),1)
if np.nansum(abs(np.asarray(shift))) < 20:
maskred[:,0] += shift[1].astype('float32')
maskred[:,1] += shift[0].astype('float32')
else:
print ('didnt work with shift: ' + str(shift))
new_img = all_red_im[:,:,plane]
#somehow it didn't work either way
print ('There was an issue with the motion correction during red_channel comparison. Please check plane: ' + str(plane))
# find distances
neur_plane = neuron_plane[plane].astype('int')
aux_nc = np.zeros(neur_plane)
aux_nc = new_com[ind_neur:neur_plane+ind_neur, :2]
aux_nerden = nerden[ind_neur:neur_plane+ind_neur]
redlabel = np.zeros(neuron_plane.shape[0]).astype('bool')
dists = np.zeros((neur_plane,maskred.shape[0]))
dists = scipy.spatial.distance.cdist(aux_nc, maskred)
# identfify neurons based on distance
redlabel = np.zeros(neur_plane).astype('bool')
aux_redlabel = np.zeros(neur_plane)
iden_pairs = [] # to debug
for neur in np.arange(neur_plane):
if aux_nerden[neur]:
aux_redlabel[neur] = np.sum(dists[neur,:]<maxdist)
redn = np.where(dists[neur,:]<maxdist)[0]
if len(redn):
iden_pairs.append([neur, redn[0]]) # to debug
redlabel[aux_redlabel>0]=True
all_red.append(redlabel)
auxtoplot = aux_nc[redlabel,:]
if toplot:
imgtoplot = np.zeros((new_img.shape[0], new_img.shape[1]))
fig1 = plt.figure()
ax1 = fig1.add_subplot(1,2,1)
for ind in np.arange(maskred.shape[0]):
auxlocx = maskred[ind,1].astype('int')
auxlocy = maskred[ind,0].astype('int')
imgtoplot[auxlocx-1:auxlocx+1,auxlocy-1:auxlocy+1] = np.nanmax(new_img)
ax1.imshow(new_img + imgtoplot, vmax=np.nanmax(new_img))
imgtoplot = np.zeros((new_img.shape[0], new_img.shape[1]))
ax2 = fig1.add_subplot(1,2,2)
for ind in np.arange(auxtoplot.shape[0]):
auxlocx = auxtoplot[ind,1].astype('int')
auxlocy = auxtoplot[ind,0].astype('int')
imgtoplot[auxlocx-1:auxlocx+1,auxlocy-1:auxlocy+1] = np.nanmax(new_img)
ax2.imshow(new_img + imgtoplot, vmax=np.nanmax(new_img))
plt.savefig(fanal + str(plane) + '/redneurmask.png', bbox_inches="tight")
fig2 = plt.figure()
R = new_img
auxA = np.unique(np.arange(Afull.shape[2])[ind_neur:neur_plane+ind_neur]*redlabel)
G = np.transpose(np.nansum(Afull[:,:,auxA],2))
B = np.zeros((R.shape))
R = R/np.nanmax(R)
G = G/np.nanmax(G)
RGB = np.dstack((R,G,B))
plt.imshow(RGB)
plt.savefig(fanal + str(plane) + '/redneurmask_RG.png', bbox_inches="tight")
plt.close("all")
ind_neur += neur_plane
all_red = np.concatenate(all_red)
return all_red
def obtain_real_com(fanal, Afull, all_com, nerden, toplot=True, img_size = 20, thres=0.1):
"""
Function to obtain the "real" position of the neuron regarding the spatial filter
fanal(str): folder where the plots will be stored
Afull(array): matrix with all the spatial components
all_com(array): matrix with the position in xyz given by caiman
nerden(array-bool): array of bool labelling as true components identified as neurons
toplot(bool): flag to plot and save
thres(int): tolerance to identify the soma of the spatial filter
minsize(int): minimum size of a neuron. Should be change for types of neurons / zoom / spatial resolution
Returns
new_com(array): matrix with new position of the neurons
"""
#function to obtain the real values of com
faplot = fanal + 'Aplot/'
if not os.path.exists(faplot):
os.makedirs(faplot)
new_com = np.zeros((Afull.shape[2], 3))
for neur in np.arange(Afull.shape[2]):
center_mass = scipy.ndimage.measurements.center_of_mass(Afull[:,:,neur]>thres)
if np.nansum(center_mass)==0 :
center_mass = scipy.ndimage.measurements.center_of_mass(Afull[:,:,neur])
new_com[neur,:] = [center_mass[0], center_mass[1], all_com[neur,2]]
if (center_mass[0] + img_size) > Afull.shape[0]:
x2 = Afull.shape[0]
else:
x2 = int(center_mass[0]+img_size)
if (center_mass[0] - img_size) < 0:
x1 = 0
else:
x1 = int(center_mass[0]-img_size)
if (center_mass[1] + img_size) > Afull.shape[1]:
y2 = Afull.shape[1]
else:
y2 = int(center_mass[1]+img_size)
if (center_mass[1] - img_size) < 0:
y1 = 0
else:
y1 = int(center_mass[1]-img_size)
if toplot:
img = Afull[x1:x2,y1:y2,neur]
fig1 = plt.figure()
ax1 = fig1.add_subplot(121)
ax1.imshow(np.transpose(img))
ax1.set_xlabel('nd: ' + str(nerden[neur]))
ax2 = fig1.add_subplot(122)
ax2.imshow(np.transpose(img>thres))
ax2.set_xlabel('neuron: ' + str(neur))
plt.savefig(faplot + str(neur) + '.png', bbox_inches="tight")
plt.close('all')
return new_com
def detect_ensemble_neurons(fanal, all_C, online_data, units, com, metadata, neuron_plane, number_planes_total, len_base, auxtol=10, cormin=0.6):
"""
Function to identify the ensemble neurons across all components
fanal(str): folder where the plots will be stored
dff(array): Dff values of the components. given by caiman
online_data(array): activity of the ensemble neurons registered on the online bmi experiment
units (int): number of neurons in the ensembles
com(array): position of the neurons
mask(array): position of the ensemble neurons as given by the experiment
number_planes_total(int): number of planes given back by the recording system, it may differ from number_planes
len_base(int): lenght of the baseline
auxtol (int): max difference distance for ensemble neurons
cormin (int): minimum correlation between neuronal activity from caiman DFF and online recording
returns
final_neur(array): index of the ensemble neurons"""
# initialize vars
neurcor = np.ones((units, all_C.shape[0])) * np.nan
finalcorr = np.zeros(units)
finalneur = np.zeros(units)
finaldist = np.zeros(units)
pmask = np.zeros((metadata['FrameData']['SI.hRoiManager.pixelsPerLine'], metadata['FrameData']['SI.hRoiManager.linesPerFrame']))
iter = 40
ind_neuron_plane = np.cumsum(neuron_plane).astype('int')
#extract reference from metadata
a = metadata['RoiGroups']['imagingRoiGroup']['rois']['scanfields']['pixelToRefTransform'][0][0]
b = metadata['RoiGroups']['imagingRoiGroup']['rois']['scanfields']['pixelToRefTransform'][0][2]
all_zs = metadata['FrameData']['SI.hStackManager.zs']
for un in np.arange(units):
print(['finding neuron: ' + str(un)])
tol = copy.deepcopy(auxtol)
tempcormin = copy.deepcopy(cormin)
for npro in np.arange(all_C.shape[0]):
ens = (online_data.keys())[2+un]
frames = (np.asarray(online_data['frameNumber']) / number_planes_total).astype('int') + len_base
auxonline = (np.asarray(online_data[ens]) - np.nanmean(online_data[ens]))/np.nanmean(online_data[ens])
auxC = all_C[npro,frames]/10000
neurcor[un, npro] = pd.DataFrame(np.transpose([auxC[~np.isnan(auxonline)], auxonline[~np.isnan(auxonline)]])).corr()[0][1]
neurcor[neurcor<tempcormin] = np.nan
auxneur = copy.deepcopy(neurcor)
# extract position from metadata
relativepos = metadata['RoiGroups']['integrationRoiGroup']['rois'][un]['scanfields']['centerXY']
centermass = np.reshape(((np.asarray(relativepos) - b)/a).astype('int'), [1,2])
zs = metadata['RoiGroups']['integrationRoiGroup']['rois'][un]['zs']
plane = all_zs.index(zs)
if plane == 0:
neurcor[un, ind_neuron_plane[0]:] = np.nan
else:
neurcor[un, :ind_neuron_plane[plane-1]] = np.nan
neurcor[un, ind_neuron_plane[plane]:] = np.nan
not_good_enough = True
while not_good_enough:
if np.nansum(neurcor[un,:]) != 0:
if np.nansum(np.abs(neurcor[un, :])) > 0 :
maxcor = np.nanmax(neurcor[un, :])
indx = np.where(neurcor[un, :]==maxcor)[0][0]
auxcom = np.reshape(com[indx,:2], [1,2])
dist = scipy.spatial.distance.cdist(centermass, auxcom)[0][0]
not_good_enough = dist > tol
finalcorr[un] = neurcor[un, indx]
finalneur[un] = indx
finaldist[un] = dist
neurcor[un, indx] = np.nan
else:
print('Error couldnt find neuron' + str(un) + ' with this tolerance. Reducing tolerance')
if iter > 0:
neurcor = auxneur
tol *= 1.5
iter -= 1
else:
print ('where are my neurons??')
break
else:
print('trying new thresholds')
if tempcormin > 0:
tempcormin-= 0.1
tol-= 2
not_good_enough = True
else:
print('No luck, finding neurons by distance')
auxcom = com[:,:2]
dist = scipy.spatial.distance.cdist(centermass, auxcom)[0]
if plane == 0:
dist[ind_neuron_plane[0]:] = np.nan
else:
dist[:ind_neuron_plane[plane-1]] = np.nan
dist[ind_neuron_plane[plane]:] = np.nan
indx = np.where(dist==np.nanmin(dist))[0][0]
finalcorr[un] = np.nan
if np.nanmin(dist) < auxtol:
finaldist[un] = np.nanmin(dist)
finalneur[un] = indx
else:
finalneur[un] = np.nan
finaldist[un] = np.nan
not_good_enough = False
print('tol value at: ', str(tol), 'correlation thres at: ', str(tempcormin))
print('Correlated with value: ', str(finalcorr[un]), ' with a distance: ', str(finaldist[un]))
if ~np.isnan(finalneur[un]):
auxp = com[finalneur[un].astype(int),:2].astype(int)
pmask[auxp[1], auxp[0]] = 2 #to detect
pmask[centermass[0,1], centermass[0,0]] = 1 #to detect
plt.figure()
plt.imshow(pmask)
plt.savefig(fanal + 'ens_masks.png', bbox_inches="tight")
fig1 = plt.figure(figsize=(16,6))
for un in np.arange(units):
ax = fig1.add_subplot(units, 1, un + 1)
ens = (online_data.keys())[2+un]
frames = (np.asarray(online_data['frameNumber']) / number_planes_total).astype('int') + len_base
auxonline = (np.asarray(online_data[ens]) - np.nanmean(online_data[ens]))/np.nanmean(online_data[ens])
auxonline[np.isnan(auxonline)] = 0
ax.plot(zscore(auxonline[-5000:]))
if ~np.isnan(finalneur[un]):
auxC = all_C[finalneur[un].astype('int'), frames]
ax.plot(zscore(auxC[-5000:]))