Anesthesia_CA1/validation/sleep_data_calcium_imaging/Validation_stability.py

#!/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("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_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 = ["animal",'number.neurons','traces.median','traces.skewness','decay','median.stability']
labels = ["Animal \n ID",'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=(8, 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


        #else:
        sns.boxplot(x='recording', y=param, data=df_estimators[(df_estimators.recording>=rmin)&(df_estimators.recording<=rmax)], width=0.9, dodge=False, showfliers = False,ax=axes[i],linewidth=lw)
        if (i == 0):
            param = "animal_cat"
            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 > 1:
            axes[i].set_ylim([0.0,2000.0])
        if i > 2:
            axes[i].set_ylim([0.0,10.0])
        if i > 3:
            axes[i].set_ylim([0.0,1.0])
        if i > 4:
            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("Recording", 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.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()