LightDisco/runme.py

import pylab as pl
import numpy as np
from os import listdir, makedirs
from os.path import isfile, join, isdir
import seaborn as sns
import scipy.optimize as opt


from Helpers.file_helpers import data_dictionary_di

experiments = [2,3,4,5,6,8,9,10,11,12,13,15,19,21,24,25]

pathtofiles = '/Users/ms81/Desktop/Disco_plots'

datapath = pathtofiles + '/Lightdisco_raw/'
figure_path = pathtofiles + '/img/'


days = []
disco_flag = []
culture_flag = []
cvs = []
info = []
rates = []
waveforms = []
for ex in experiments:
    print('Collecting...' + str(ex))
    
    di = data_dictionary_di(ex, verbose = 'True', infofile = '/Users/ms81/Desktop/Disco_plots/data.txt')
    age = di[ex]['DIV']
    discoflag = di[ex]['disco']
    good_channels_file = di[ex]['channels']
    culture = di[ex]['culture']
    
    saveresults = datapath + 'rawdata_' + str(ex)
    a = np.load(saveresults + '.npy', encoding='latin1')
    final_results = a.item()
    elenumber = final_results['data_dict']['elenumber']
    data = final_results['electrode_data']
    sorted_dict = final_results['sorted_dict']
    mean_waveforms = sorted_dict['mean_waveforms_good']

    for ele in range(len(elenumber)):
        if elenumber[ele] in good_channels_file:
            spikes = data['Electrode_' + str(elenumber[ele])]
            mu = np.mean(np.diff(spikes))
            sigma = np.std(np.diff(spikes))
            cv = sigma/mu
            rate = len(spikes)/np.max(spikes)   #number of spike, until the last one
            print( str(elenumber[ele]) + '  ' + str(cv) + '  ' + str(rate) )
            cvs.append(cv)
            rates.append(rate)
            days.append(age)
            info.append(elenumber[ele])
            disco_flag.append(discoflag)
            culture_flag.append(culture)
            waveforms.append(mean_waveforms[ele])
print('Collection done, now printing.')


#~ sns.set_context("talk")
sns.set_style("whitegrid")
colors = sns.color_palette("Set3",n_colors=7)

days_ref = []
for s in days:
    if s ==12 or s == 13 or s == 15 or s == 16:
        days_ref.append(14)
    if s == 20 or s == 23:
        days_ref.append(21)
    if s == 9:
        days_ref.append(7)


pl.figure(1, figsize = (15,5))
pl.subplot(1,2,1)
pl.title('CVs, Lightdisco vs. controls')
sns.boxplot(x=days_ref, y=cvs, hue = disco_flag, palette=colors)
pl.ylim([0,6])
pl.xlabel('Age [DIV]')
pl.ylabel('Coefficient of Variation')
pl.subplot(1,2,2)
pl.title('Firing rates, Lightdisco vs. controls')
sns.boxplot(x=days_ref, y=rates, hue = disco_flag, palette="Set3")
pl.ylim([0,5])
pl.xlabel('Age [DIV]')
pl.ylabel('Rates')
pl.savefig(figure_path + 'fig1.pdf', bbox_inches='tight')
pl.close()




pl.figure(2, figsize = (15,5))
pl.subplot(1,2,1)
for i in range(len(days)):
    if disco_flag[i] == 'yes':
        pl.plot(days[i],cvs[i],'.', markerfacecolor = colors[int(culture_flag[i][0])])
    else:
        pl.plot(days[i]+.2,cvs[i],'b.')
    #~ pl.text(days[i]+.5,cvs[i],info[i])
pl.xlabel('DIV')
pl.ylabel('Coefficient of Variation')
pl.ylim([0,6])
pl.subplot(1,2,2)
for i in range(len(days)):
    if disco_flag[i] == 'yes':
        pl.plot(days[i],rates[i],'.', markerfacecolor = colors[int(culture_flag[i][0])])
    else:
        pl.plot(days[i]+.2,rates[i],'b.')
    #~ pl.text(days[i]+.5,cvs[i],info[i])
pl.xlabel('DIV')
pl.ylabel('Firing Rate')
pl.ylim([0,5])

pl.savefig(figure_path + 'fig2.pdf', bbox_inches='tight')
pl.close()


pl.figure(3, figsize = (15,5))
for i in range(len(days)):
    if disco_flag[i] == 'yes':
        pl.plot(cvs[i],rates[i],'r.')
    else:
        pl.plot(cvs[i],rates[i],'b.')
pl.xlabel('CVs')
pl.ylabel('Rate')
pl.xlim([0,5])
pl.ylim([0,15])
pl.savefig(figure_path + 'fig3.pdf', bbox_inches='tight')
pl.close()





#This is the waveform figure


def make_gaussian(xx, x0, stdx, a):
    xx =  xx - x0
    gauss = a*np.exp(-(xx*xx/(2.*stdx*stdx)))
    return gauss


ref = []
fit = []
fit_error = []
discoshifted = []
controlshifted = []
c1, c2 = sns.color_palette("Set3", 2)


pl.figure(4, figsize = (8,10))
xx = np.linspace(0, 6, 150)    #0ms to 6 ms in 40mus steps
for i in range(len(rates)):
    wv = waveforms[i]
    normalized = wv*1e6 #-wv/np.min(wv)
    minwhere_ms = np.argmin(normalized)*6./150


    popt, pcov = opt.curve_fit(make_gaussian, xx, wv, p0=(2,0.2,-1e-4))
    perr = np.sqrt(np.diag(pcov))
    gaussian_fitted = make_gaussian(xx, *popt)
    fit.append(popt)
    fit_error.append(perr)

    if disco_flag[i] == 'yes':
        pl.subplot(3,2,1)        
        pl.plot(xx-minwhere_ms, normalized, color = c1, alpha = .1)
        pl.ylim([-130,70])
        ref.append(0)
        wvsh = np.array([0 for i in range(300)])
        start = int(minwhere_ms*150/6.)
        wvsh[start:(150+start)] = normalized
        discoshifted.append(wvsh)

    else:
        pl.subplot(3,2,2)
        pl.plot(xx-minwhere_ms, normalized, color = c2, alpha = .05)
        pl.ylim([-130,70])
        ref.append(1)
        wvsh = np.array([0 for i in range(300)])
        start = int(minwhere_ms*150/6.)
        wvsh[start:(150+start)] = normalized
        controlshifted.append(wvsh)
        
fit = np.array(fit)
pl.subplot(3,2,3)
sns.boxplot(x=ref, y=np.abs(fit[:,1]), palette="Set3")
pl.ylabel("width [ms]")
pl.ylim([0,0.35])

pl.subplot(3,2,4)
sns.boxplot(x=ref, y=np.abs(fit[:,2])*1e6, palette="Set3")
pl.ylabel("Amplitude [muV]")
pl.ylim([0,160])

pl.subplot(3,2,5)
pl.fill_between(x=np.arange(300)*6./150, y1=np.mean(discoshifted,axis=0) - np.std(discoshifted,axis=0), y2=np.mean(discoshifted,axis=0) + np.std(discoshifted,axis=0),alpha=0.5,facecolor = c1)
pl.plot(np.arange(300)*6./150,np.mean(discoshifted,axis=0),color = c1)
pl.fill_between(x=np.arange(300)*6./150, y1=np.mean(controlshifted,axis=0) - np.std(controlshifted,axis=0), y2=np.mean(controlshifted,axis=0) + np.std(controlshifted,axis=0), alpha=0.5, facecolor = c2)
pl.plot(np.arange(300)*6./150,np.mean(controlshifted,axis=0), color = c2)
pl.xlim([2,8])


pl.subplot(3,2,6)
controldistro = np.mean(controlshifted,axis=0)
dt = 6/150.
ref = []
dd = []
for i in range(len(rates)):
    wv = waveforms[i]
    normalized = wv*1e6 #-wv/np.min(wv)
    minwhere_ms = np.argmin(normalized)*6./150
    if disco_flag[i] == 'yes':
        ref.append(0)
        wvsh = np.array([0 for i in range(300)])
        start = int(minwhere_ms*150/6.)
        wvsh[start:(150+start)] = normalized
        dd.append(np.sqrt(dt*np.sum(controldistro - wvsh)**2))
    else:
        ref.append(1)
        wvsh = np.array([0 for i in range(300)])
        start = int(minwhere_ms*150/6.)
        wvsh[start:(150+start)] = normalized
        dd.append(np.sqrt(dt*np.sum(controldistro - wvsh)**2))


sns.boxplot(x=ref, y=dd, palette="Set3")
pl.ylabel("Deviation from control mean")
pl.savefig(figure_path + 'fig4.pdf', bbox_inches='tight')
pl.close()