#!/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=(0.75, 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.swarmplot(x='animal', y='number.neurons', data=df_nneurons[(df_nneurons.recording>=rmin)&(df_nneurons.recording<=rmax)], ax=axes[i], marker = '.',size=1,edgecolor="black", linewidth = lw*0.5) # #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([-100.0,1500.0]) if i > 1: axes[i].set_ylim([-3.0,15.0]) if i > 2: axes[i].set_ylim([-0.5,1.5]) if i > 3: axes[i].set_ylim([90,110]) 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_sleep_figure2.png",dpi=400) plt.savefig("Validation_stability_check_sleep_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]: