#!/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/Anesthesia_CA1/meta_data/meta_recordings_transition_state.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,198+1)] #rec = [127,128] # 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,'mmf':color_mmf,'awake':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) print("Shape:" + str(Traces.shape[0]) + "N_accept_reject" + str(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: " + str(r) + "Recording length:" + str(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("./transition_state_calcium_imaging_stability_validation.pkl") ''' # ### Load the result of the analysis # In[8]: df_estimators = pd.read_pickle("./transition_state_calcium_imaging_stability_validation.pkl") df_estimators['neuronID'] = df_estimators.index #df_estimators["animal"] = df_estimators["animal"] print(np.unique(df_estimators["condition"])) df_estimators["CONDITION"] = df_estimators["condition"] df_estimators.loc[:,"CONDITION"] = 'awake' df_estimators.loc[(df_estimators.condition == 'iso'),"CONDITION"] = 'iso' df_estimators.loc[(df_estimators.condition == 'keta'),"CONDITION"] = 'keta' df_estimators.loc[(df_estimators.condition == 'mmf'),"CONDITION"] = 'mmf' print(np.unique(df_estimators["CONDITION"])) df_estimators["multihue"] = df_estimators["CONDITION"] + df_estimators["animal"].astype("string") print(np.unique(df_estimators["multihue"])) print(np.unique(df_estimators["CONDITION"])) # ### 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,198]] for rmin,rmax in recordings_ranges: f, axes = plt.subplots(number_subplots, 1, figsize=(8, 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 #else: sns.boxplot(x='multihue', y=param, data=df_estimators[(df_estimators.recording>=rmin)&(df_estimators.recording<=rmax)], width=0.9, hue = "CONDITION", palette = custom_palette, dodge=False, showfliers = False,ax=axes[i],linewidth=lw) #if (i == 0): # param = "animal" # print(np.unique(df_estimators[param])) # axes[i].set_yticks(np.unique(df_estimators[param])) # sns.swarmplot(x='recording', y=param, data=df_estimators[(df_estimators.recording>=rmin)&(df_estimators.recording<=rmax)&(df_estimators['neuronID'] == 0)],dodge=False, s=1, edgecolor='black', linewidth=1, ax=axes[i]) #ax.set(ylabel="") 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=9,labelpad=5) #40 #axes[i].set_xlabel("Recording", fontsize=5,labelpad=5) #40 #axes[i].axis('off') axes[i].xaxis.set_tick_params(labelsize=9) #35 axes[i].yaxis.set_tick_params(labelsize=9) #30 axes[i].get_legend().remove() #axes[i].xaxis.set_major_locator(ticker.MultipleLocator(10)) #axes[i].xaxis.set_major_formatter(ticker.ScalarFormatter()) plt.xlabel('xlabel', fontsize=6) plt.xticks(rotation=90) #plt.legend(bbox_to_anchor=(1.01, 1), loc=2, borderaxespad=0.,fontsize=25) plt.savefig("Validation_stability_check_rec_#%d-#%d).png" % (rmin,rmax),dpi=400) plt.savefig("Validation_stability_check_rec_#%d-#%d).svg" % (rmin,rmax)) #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]: '''