gin commit from andrey-UX510UXK

New files: 160
Modified files: 9
Deleted files: 3

validation/sleep_data_calcium_imaging/Figure2_compact.py

@@ -0,0 +1,261 @@
+#!/usr/bin/env python
+# coding: utf-8
+
+# ### Link to the file with meta information on recordings
+
+# In[14]:
+
+#import matplotlib.pyplot as plt
+#plt.rcParams["figure.figsize"] = (20,3)
+
+
+database_path = '/media/andrey/My Passport/GIN/backup_Anesthesia_CA1/meta_data/meta_recordings_sleep.xlsx'
+
+
+# ### Select the range of recordings for the analysis (see "Number" row in the meta data file)
+
+# In[4]:
+
+rec = [x for x in range(0,74+1)]
+#rec = [1,2]
+
+# In[1]:
+
+
+import numpy as np
+
+import numpy.ma as ma
+
+import matplotlib.pyplot as plt
+import matplotlib.ticker as ticker
+
+import pandas as pd
+import seaborn as sns
+import pickle
+import os
+
+sns.set()
+sns.set_style("whitegrid")
+
+from scipy.signal import medfilt 
+
+from scipy.stats import skew, kurtosis, zscore
+
+from scipy import signal
+
+from sklearn.linear_model import LinearRegression, TheilSenRegressor
+
+plt.rcParams['figure.figsize'] = [16, 8]
+
+color_awake = (0,191/255,255/255)
+color_mmf = (245/255,143/255,32/255)
+color_keta = (181./255,34./255,48./255)
+color_iso = (143./255,39./255,143./255)
+
+custom_palette ={'keta':color_keta, 'iso':color_iso,'fenta':color_mmf,'awa':color_awake}
+
+
+# In[2]:
+
+
+from capipeline import *
+
+
+
+# ### Run the analysis
+# /media/andrey/My Passport/GIN/Anesthesia_CA1/validation/calcium_imaging
+# It creates a data frame *df_estimators* that contains basic information regarding stability of the recordings, such as 
+# 
+# - total number of identified neurons,
+# - traces and neuropils median inntensities for each ROI
+# - their standard deviation
+# - skewness of the signal
+# - estamation of their baseline (defined as a bottom quartile of signal intensities) 
+# - their temporal stability (defined as the ratio between median signals of all ROIs in the first and the second parts of the recording)  
+
+# In[5]:
+
+'''
+df_estimators = pd.DataFrame()
+
+
+for r in rec:
+    
+    Traces, Npils, n_accepted_and_rejected = traces_and_npils(r, database_path, concatenation=False)
+
+    print(Traces.shape[0], "vs", n_accepted_and_rejected)
+    
+    animal = get_animal_from_recording(r, database_path)
+        
+    #condition = get_condition(r, database_path)
+    
+    #print("#" +  str(r) + " " + str(animal) + " " + str(condition) + " ")
+    
+    Traces_median = ma.median(Traces, axis=1) 
+    Npils_median = ma.median(Npils, axis=1)
+    
+    Traces_std = ma.std(Npils, axis=1)    
+    Npils_std = ma.std(Npils, axis=1)
+    
+    Traces_skewness = skew(Traces,axis=1)
+    Npils_skewness = skew(Npils,axis=1)
+    
+    baseline = np.quantile(Traces,0.25,axis=1)
+
+
+    num_cells = np.shape(Traces)[0]
+    decay_isol = np.zeros((num_cells))
+     
+    fs = 30 
+        
+    for neuron in np.arange(num_cells):
+            
+        if np.all(np.isnan(Traces[neuron])):
+            decay_isol[neuron] = np.nan
+        else:
+            _, _, _, decay_neuron_isolated10, _ = deconvolve(np.double(Traces[neuron, ] + 100000),
+                                                                 penalty = 0, optimize_g = 10)
+            decay_isol[neuron] = - 1 / (fs * np.log(decay_neuron_isolated10))
+
+    
+    recording_length = int(Traces.shape[1])
+    
+    half = int(recording_length/2)
+    
+    print(recording_length)
+    m1 = ma.median(Traces[:,:half])
+    m2 = ma.median(Traces[:,half:])
+
+    print("Stability:",m2/m1*100)
+
+    norm_9000 = 9000/recording_length   # normalize to 9000 frames (5 min recording)
+    
+    traces_median_half_vs_half = norm_9000*(m2-m1)*100/m1 + 100
+
+    print("Stability (9000 frames normalization)",traces_median_half_vs_half)
+
+    df_e = pd.DataFrame({ "animal":animal,
+                        "recording":r,
+                        #"condition":condition,
+                        "number.neurons":Traces.shape[0],
+                        "traces.median":Traces_median,
+                        "npils.median":Npils_median,
+                        "traces.std":Traces_std,
+                        "npils.std":Npils_std,
+                        
+                        "traces.skewness":Traces_skewness,
+                        "npils.skewness":Npils_skewness,
+     
+                        "baseline.quantile.25":baseline,
+
+                        "decay":decay_isol,
+                         
+                        "median.stability":traces_median_half_vs_half # in percent
+
+                      })
+        
+    df_estimators = pd.concat([df_estimators,df_e])
+
+    print("*****")
+   
+
+
+# ### Save the result of the analysis
+
+# In[7]:
+
+
+df_estimators.to_pickle("./sleep_data_calcium_imaging_stability_validation.pkl") 
+
+'''
+# ### Load the result of the analysis
+
+# In[8]:
+
+
+df_estimators = pd.read_pickle("./sleep_data_calcium_imaging_stability_validation.pkl") 
+df_estimators['neuronID'] = df_estimators.index
+
+df_estimators["animal_cat"] = df_estimators["animal"].astype("category")
+
+
+# ### Plot 
+
+# In[9]:
+
+parameters = ['number.neurons','traces.median','traces.skewness','decay','median.stability']
+labels = ['Extracted \n ROIs','Median, \n A.U.','Skewness','Decay time, \n s','1st/2nd \n ratio, %']
+number_subplots = len(parameters)
+
+recordings_ranges = [[0,74]]
+
+for rmin,rmax in recordings_ranges:
+
+    f, axes = plt.subplots(number_subplots, 1, figsize=(1, 4.5)) # sharex=Truerex=True
+    #plt.subplots_adjust(left=None, bottom=0.1, right=None, top=0.9, wspace=None, hspace=0.2)
+    #f.tight_layout()
+    sns.despine(left=True)
+
+    for i, param in enumerate(parameters):
+        
+        lw = 0.8
+
+        df_nneurons = df_estimators.groupby(['recording','animal'], as_index=False)['number.neurons'].median()
+
+        if (i == 0):
+            sns.scatterplot(x='animal', y='number.neurons', data=df_nneurons[(df_nneurons.recording>=rmin)&(df_nneurons.recording<=rmax)],  ax=axes[i], marker = '.',edgecolor="black", linewidth = lw*0.5)  #
+        else:
+            sns.boxplot(x='animal', y=param, data=df_estimators[(df_estimators.recording>=rmin)&(df_estimators.recording<=rmax)], color = (15./255.,117./255.,188./255.),showfliers = False,ax=axes[i],linewidth=lw)
+
+        if i == 0:
+            axes[i].set_ylim([0.0,1200.0])       
+        if i > 0:
+            axes[i].set_ylim([0.0,1000.0])
+        if i > 1:
+            axes[i].set_ylim([0.0,10.0])
+        if i > 2:
+            axes[i].set_ylim([0.0,1.0])
+        if i > 3:
+            axes[i].set_ylim([80,120])
+            axes[i].get_xaxis().set_visible(True)
+        else:
+            axes[i].get_xaxis().set_visible(False)
+           
+        if i < number_subplots-1:
+            axes[i].xaxis.label.set_visible(False)
+        #if i==0:
+        #    axes[i].set_title("Validation: stability check (recordings #%d-#%d)" % (rmin,rmax), fontsize=9, pad=30) #45
+        axes[i].set_ylabel(labels[i], fontsize=8,labelpad=5) #40
+        axes[i].set_xlabel("Animal", fontsize=8,labelpad=5) #40
+        #axes[i].axis('off')
+
+        axes[i].xaxis.set_tick_params(labelsize=6) #35
+        axes[i].yaxis.set_tick_params(labelsize=6) #30
+
+        #axes[i].get_legend().remove()
+        #axes[i].xaxis.set_major_locator(ticker.MultipleLocator(10))
+        #axes[i].xaxis.set_major_formatter(ticker.ScalarFormatter())
+
+    #plt.legend(loc='upper right',fontsize=35)
+    #plt.legend(bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0.,fontsize=25)
+    plt.xticks(rotation=45)
+    plt.savefig("Validation_stability_check_figure2.png",dpi=400)
+    plt.savefig("Validation_stability_check_figure2.svg")
+    #plt.show()
+
+
+# In[13]:
+
+
+sns.displot(data=df_estimators, x="median.stability", kind="kde", hue = "animal")
+plt.xlim([80,120])
+plt.xlabel("Stability, %", fontsize = 15)
+plt.title("Validation: summary on stability (recordings #%d-#%d)" % (min(rec),max(rec)), fontsize = 20, pad=20)
+plt.grid(False)
+plt.savefig("Validation_summary_stability_recordings_#%d-#%d)" % (min(rec),max(rec)))
+#plt.show()
+
+
+# In[62]:
+
+

validation/sleep_data_calcium_imaging/ROI_adjustment_check.py

@@ -0,0 +1,176 @@
+#!/usr/bin/env python
+# coding: utf-8
+
+# ### Define paths and animals for the analysis
+
+# In[1]:
+
+
+path = '/media/andrey/My Passport/GIN/backup_Anesthesia_CA1/meta_data/meta_recordings_sleep.xlsx'
+path4results = '/media/andrey/My Passport/GIN/Anesthesia_CA1/validation/sleep_data_calcium_imaging/' #To store transformation matrix
+save_plots_path = '/media/andrey/My Passport/GIN/Anesthesia_CA1/validation/sleep_data_calcium_imaging/'
+log_file_path = save_plots_path + 'registration_logs.txt'
+
+animals_for_analysis = [8235,8237,8238]
+
+# ### Align FOV's for all recordings
+
+# In[2]:
+
+repeat_calc = 1
+silent_mode = False
+
+#######################
+
+import pandas as pd
+import numpy as np
+import os
+
+import matplotlib.pyplot as plt
+
+import sys
+np.set_printoptions(threshold=sys.maxsize)
+
+from pystackreg import StackReg
+
+# Sobel filter (not used)
+#from scipy import ndimage #https://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.sobel.html
+
+meta_data = pd.read_excel(path, engine='openpyxl')
+
+#%% compute transformations matrices between recordings
+
+recordings = meta_data['Number']
+
+animals = animals_for_analysis
+
+#print("Recordings: ", recordings)
+
+#ALIGNMENT ALGORITHM
+
+sr = StackReg(StackReg.AFFINE)
+
+
+# ROI adjustment check for all recordings
+
+from matplotlib import colors
+Amap = colors.ListedColormap(['white','tab:blue', 'tab:red','purple'])
+
+
+if not os.path.exists(save_plots_path + 'ROIAdjustmentCheck/'):
+    os.makedirs(save_plots_path + 'ROIAdjustmentCheck/')
+'''   
+for animal in animals:
+    if not os.path.exists(save_plots_path + 'ROIAdjustmentCheck/' + str(animal) + "/"):
+        os.makedirs(save_plots_path + 'ROIAdjustmentCheck/' + str(animal) + "/")  
+    
+    tmats_loaded = np.load(path4results + 'StackReg/' + str(animal) + "_best_tmats" + '.npy')
+    meta_animal = meta_data[meta_data['Mouse'] == animal]
+    recordings = meta_animal['Number'] #,'blue', 'green','cyan'])
+    print(recordings)
+    images = np.zeros((512, 512, np.shape(meta_animal)[0]))
+
+    for idx, recording in enumerate(recordings):
+        stats = np.load(meta_data['Folder'][recording] +
+                        str(meta_data['Subfolder'][recording]) +
+                        str(int(meta_data['Recording idx'][recording])) +
+                        '/suite2p/plane0/stat.npy',
+                        allow_pickle=True)
+        iscell = np.load(meta_data['Folder'][recording] +
+                         str(meta_data['Subfolder'][recording]) +
+                         str(int(meta_data['Recording idx'][recording])) +
+                         '/suite2p/plane0/iscell.npy',
+                         allow_pickle=True)
+        stats = stats[iscell[:, 0].astype(bool)]
+
+        for stats_neuron in stats:
+                images[stats_neuron['ypix'], stats_neuron['xpix'], idx] += 1
+        images[:,:,idx][images[:,:,idx]>1] = 1
+        
+    roi_corresp = np.zeros((512, 512))
+ 
+    #for idx0 in range(np.shape(images)[2]):
+        #for idx1 in range(idx0, np.shape(images)[2]):
+    for idx0, recording0 in enumerate(recordings):
+        for idx1, recording1 in enumerate(recordings):
+           # roi_corresp = images[:,:,idx0] + 2*images[:,:,idx1]
+            roi_corresp = images[:,:,idx0] + 2*sr.transform(images[:, :, idx1], tmats_loaded[idx0, idx1, :, :])
+            file_title = str(idx0) + "_" + str(idx1) + "_" + meta_data['Condition'][recording0] + "_" + meta_data['Condition'][recording1]
+            print(file_title)
+            plt.imshow(roi_corresp[:,:],cmap=Amap) #'tab10'
+            #plt.annotate('blue - reference / green - corrected / cyan - intersection', xy=(0, 0), xytext=(.8, 0), fontsize=12)
+            plt.savefig(save_plots_path + "ROIAdjustmentCheck/" + str(animal) + '/' + file_title + "_ROI_correspondence.svg")
+            #plt.show()
+'''
+
+#####################################################################################
+############################ Draw ROI in all recordings #############################
+#####################################################################################
+
+# - check if the iscell is active
+# - check if the alignment is still active 
+
+if not os.path.exists(save_plots_path + 'ROIs/'):
+    os.makedirs(save_plots_path + 'ROIs/')
+   
+for animal in animals:
+    if not os.path.exists(save_plots_path + 'ROIs/' + str(animal) + "/"):
+        os.makedirs(save_plots_path + 'ROIs/' + str(animal) + "/")  
+    
+    tmats_loaded = np.load(path4results + 'StackReg/' + str(animal) + "_best_tmats" + '.npy')
+    meta_animal = meta_data[meta_data['Mouse'] == animal]
+    recordings = meta_animal['Number'] #,'blue', 'green','cyan'])
+    print(recordings)
+    images = np.zeros((512, 512, np.shape(meta_animal)[0]))
+
+    for idx, recording in enumerate(recordings):
+        stats = np.load(meta_data['Folder'][recording] +
+                        str(meta_data['Subfolder'][recording]) +
+                        #str(int(meta_data['Recording idx'][recording])) +
+                        '/suite2p/plane0/stat.npy',
+                        allow_pickle=True)
+        iscell = np.load(meta_data['Folder'][recording] +
+                         str(meta_data['Subfolder'][recording]) +
+                         #str(int(meta_data['Recording idx'][recording])) +
+                         '/suite2p/plane0/iscell.npy',
+                         allow_pickle=True)
+
+        #stats = stats[iscell[:, 0].astype(bool)]   ### all ROI or only classified as neurons
+
+        for stats_neuron in stats:
+                images[stats_neuron['ypix'], stats_neuron['xpix'], idx] += 1
+        images[:,:,idx][images[:,:,idx]>1] = 1
+        
+    roi_corresp = np.zeros((512, 512))
+
+    for idx0, recording0 in enumerate(recordings):
+    
+            #roi_corresp = sr.transform(images[:, :, idx0], tmats_loaded[0, idx0, :, :])  ### Alignment with respect to 0 recording
+            roi_corresp = images[:, :, idx0] # No alignment
+            #roi_corresp [roi_corresp > 0.5] = 1
+            file_title = str(idx0)
+            print(file_title)
+            plt.imshow(roi_corresp[:,:],cmap='binary') #'tab10'
+
+            plt.savefig(save_plots_path + "ROIs/" + str(animal) + '/' + file_title + "_ROIs.svg")
+
+
+    
+    
+    
+    
+    
+    
+    
+    
+    
+    
+    
+    
+    
+    
+
+
+
+
+