Anesthesia_CA1/validation/calcium_imaging_transition_state/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/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]:
'''