PACK_in-vivo/script2_linelengths_of_snippets_Fig2.py

# -*- coding: utf-8 -*-
"""
Created on Wed Oct 14 11:18:44 2020

@author: kleisp
This is a script for getting matrices of snippets around the stimulation pulses in PACK-expressing mice, 
which can later be used for plotting the snippets on top of each other and for line length
calculations.
Note that recordings with 0.1 Hz and 0.05 Hz illumination need to be treated separately
One needs to load dictionaries containing recIDs with corresponding smr filenames for using raw data
Abbreviations:
dict- dictionary
stim- light pulse
pre- before light pulse
post-after light pulse
snipmat- matrix containing snippets

"""

import os
os.chdir('C:\EAfiles\MyCode')
from glob import glob
import numpy as np
import pickle
import pandas as pd
from scipy import stats

os.chdir('C:\EAfiles\EAdetection')

import core.ea_management as eam
reload(eam)

#%% (1) define functions to extract snippets before and after the stim
''' 'pre' or 'post' are the snippet length in seconds, here 2 is used, sr=sampling rate (here 500)'''

# get the spike trace before stimulation (pre):
#stimtime is one element from stimtimes array
def get_spiketrace_pre(stimtimes, datatrace, pre, sr):
    cutout_pre = np.int(pre*sr)
    stim_index= np.int(stimtimes*sr) 
    spiketrace = datatrace[stim_index-cutout_pre:stim_index]
    return spiketrace
  
# get the spike trace after(post) stimulation:

def get_spiketrace_post(stimtime, datatrace, post, sr):
    cutout_post = np.int(post*sr)
    stim_index= np.int(stimtime*sr)
    spiketrace = datatrace[stim_index:stim_index+cutout_post]
    return spiketrace

#get snippet before each stimulation and insert it into a matrix, called snipmat:
  
def get_snipmat_pre(datatrace, recID, pre, stimtimes,  sr):
    cutout_pre = np.int(pre*sr)
    N_stim = len(stimtimes) #number of stimuli
    snippts = cutout_pre
    snipmat = np.zeros((N_stim,snippts))
    for tt,stimtime in enumerate(stimtimes):
        spiketrace = get_spiketrace_pre(stimtime, datatrace, pre, sr) # instead of datatrace you can use zTrace
        snipmat[tt] = spiketrace
    return snipmat  
  
#get snippet after each stimulation and insert it into a matrix, called snipmat:
def get_snipmat_post(datatrace, recID, post, stimtimes,  sr):
    cutout_post = np.int(post*sr)
    N_stim = len(stimtimes) #number of stimuli
    snippts = cutout_post
    snipmat = np.zeros((N_stim,snippts))
    for tt,stimtime in enumerate(stimtimes):
        spiketrace = get_spiketrace_post(stimtime, datatrace, post, sr) # stimtime is just one timepoint from stimtimes
        snipmat[tt] = spiketrace
    return snipmat   
  
#%% (2) defining the function to get line length of each snippet(snipmat line length)

def get_snipmat_linelength(snipmat_pre, snipmat_post, linelength_dict_pre, linelength_dict_post, recID):  
  
        linelength_list=[]
        linelength_s_list=[]
        rms_list=[]
        count=0
        for i in xrange(snipmat_pre.shape[0]):
            data=snipmat_pre[i]
            N_datapoints=len(data)
            duration=N_datapoints/500 #sr=500
            linelength = np.sum(np.abs(np.diff(data)))
            linelength_list.append(linelength)
            linelength_s=linelength/duration
            linelength_s_list.append(linelength_s)
            count=count+1
            
            linelength_dict_pre[recID]={}
            
        linelength_dict_pre[recID]['sem pre-pulse linelength/second']=stats.sem(linelength_s_list) 
        linelength_dict_pre[recID]['mean pre-pulse linelength/second']=np.mean(linelength_s_list) 
        linelength_dict_pre[recID]['mean pre-pulse linelength']=np.mean(linelength_list)   
        linelength_dict_pre[recID]['pulse number']=count
        
            
            
        linelength_list=[]
        linelength_s_list=[]
        rms_list=[]
        count=0
        for i in xrange(snipmat_post.shape[0]):
            data=snipmat_post[i]
            N_datapoints=len(data)
            duration=N_datapoints/500
            linelength = np.sum(np.abs(np.diff(data)))
            linelength_list.append(linelength)
            linelength_s=linelength/duration
            linelength_s_list.append(linelength_s)
            count=count+1
            linelength_dict_post[recID]={}

        
        linelength_dict_post[recID]['sem post-pulse linelength/second']=stats.sem(linelength_s_list) 
        linelength_dict_post[recID]['mean post-pulse linelength/second']=np.mean(linelength_s_list) 
        linelength_dict_post[recID]['mean post-pulse linelength']=np.mean(linelength_list)  
        linelength_dict_post[recID]['pulse number']=count
        
#%% in order to exclude unwanted recordings
      
pathtobadrecIDs = 'C:\EAfiles\ID_Dateien\\badrecIDs.txt' 
badrecIDs = []
with open (pathtobadrecIDs,'r') as f:
    badIDs = f.read().splitlines()
    badrecIDs = badrecIDs + badIDs
    

#%% (3) Example of using the script to extract the snipmats
#Saline PACK mice 0.05Hz: get snipmat dicts (for before and after 0.05 Hz pulses) (downsampled data)

os.chdir('C:\EAfiles\EAdetection')
ymlfiles = glob('C:\EAfiles\yml_setup_PK43\*.yml') #takes all yml files from this folder
ymlfile = ymlfiles
lines = np.arange(0, len(ymlfile)) #as many lines as there are ymlfiles
snipmat_005Hz_pre_controls={}
snipmat_005Hz_post_controls={}
sr= 500
pre, post = 2, 2 #2 seconds before and after pulse    
for i in lines:
    recID = os.path.basename(ymlfile[i]).split('__')[0] #getting the IDs
    print (i)
    print (recID)
  
    if recID in badrecIDs: 
        pass
    
    else:
        excelpath = 'C:\EAfiles\Excel\\startstop_controls_HCc.xlsx'
        df = pd.read_excel(excelpath) #reading an excel file with start and stop times of each recording
        df.set_index('recID', inplace=True) #recID
        stimstop = int(df.stop[recID]) #stop time
        if recID=='PK16_005Hz1_HCc': #adding an exception (can skip this step if you have no exceptions)
            stimstart=int(df.start[recID]) 
        else:
            stimstart=10 #the stim usually starts at 10s
        stimstep=20 #speriod is 20 s
        stimtimes=np.arange(stimstart,stimstop,stimstep)
        
        aRec = eam.Rec(init_ymlpath=ymlfile[i])
        snipmat_005Hz_pre_controls[recID]={}
        data=aRec.raw_data 
        #z_data=(data-np.mean(data))/np.std(data)
        snipmat_pre=get_snipmat_pre(data, recID, pre, stimtimes, sr=500)
        snipmat_005Hz_pre_controls[recID]=snipmat_pre

        snipmat_005Hz_post_controls[recID]={}
        snipmat_post=get_snipmat_post(data, recID, post, stimtimes, sr=500)
        snipmat_005Hz_post_controls[recID]=snipmat_post

        
# save snipmat_dict 

filename = 'C:\EAfiles\MyDicts\snipmat_005Hz_pre_controls'
outfile = open(filename, 'wb') #write binary: mit Erlaubnis zu (über)schreiben
pickle.dump(snipmat_005Hz_pre_controls, outfile)
outfile.close()

filename = 'C:\EAfiles\MyDicts\snipmat_005Hz_post_controls'
outfile = open(filename, 'wb') #write binary: mit Erlaubnis zu (über)schreiben
pickle.dump(snipmat_005Hz_post_controls, outfile)
outfile.close()

#%% (4) Example of using the script to create dictionaries with line length before and after 0.05Hz pulses from snipmats above

os.chdir('C:\EAfiles\EAdetection')
ymlfiles = glob('C:\EAfiles\yml_setup_PK43\*.yml') #takes all yml files from this folder
ymlfile = ymlfiles
lines = np.arange(0, len(ymlfile)) #as many lines as there are ymlfiles
linelength_dict_005Hz_pre={}
linelength_dict_005Hz_post={}


filename = 'C:\EAfiles\MyDicts\snipmat_005Hz_pre_controls'
infile = open(filename,'rb') #read binary
snipmat_005Hz_pre_controls = pickle.load(infile)
infile.close()


filename = 'C:\EAfiles\MyDicts\snipmat_005Hz_post_controls'
infile = open(filename,'rb') #read binary
snipmat_005Hz_post_controls = pickle.load(infile)
infile.close()


for i in lines:
    recID = os.path.basename(ymlfile[i]).split('__')[0] #getting the IDs
    print (i)
    print (recID)
    
    if recID in badrecIDs: #template also needs to be excluded
        #badrecIDs = ['EPxy_xHz_HCi1','...']
        pass
    
    else:
        snipmat_pre=snipmat_005Hz_pre_controls[recID]
        snipmat_post=snipmat_005Hz_post_controls[recID]
        get_snipmat_linelength(snipmat_pre, snipmat_post, linelength_dict_005Hz_pre, linelength_dict_005Hz_post, recID)
        
        
            
filename = 'C:\EAfiles\MyDicts\linelength_dict_005Hz_post'
outfile = open(filename, 'wb') #write binary: mit Erlaubnis zu (über)schreiben
pickle.dump(linelength_dict_005Hz_post, outfile)
outfile.close()

filename = 'C:\EAfiles\MyDicts\linelength_dict_005Hz_pre'
outfile = open(filename, 'wb') #write binary: mit Erlaubnis zu (über)schreiben
pickle.dump(linelength_dict_005Hz_pre, outfile)
outfile.close()

#%% make an excel file of the pre- and post-pulse linelengths (mean and SEM)
recIDs = np.array(linelength_dict_005Hz_pre.keys())

flavour = 'mean pre-pulse linelength/second'
mean_pre_linelengths = np.array([linelength_dict_005Hz_pre[recID][flavour] for recID in linelength_dict_005Hz_pre.keys()])

flavour = 'sem pre-pulse linelength/second'
sem_pre_linelengths = np.array([linelength_dict_005Hz_pre[recID][flavour] for recID in linelength_dict_005Hz_pre.keys()])

flavour = 'mean post-pulse linelength/second'
mean_post_linelengths = np.array([linelength_dict_005Hz_post[recID][flavour] for recID in linelength_dict_005Hz_pre.keys()])

flavour = 'sem post-pulse linelength/second'
sem_post_linelengths = np.array([linelength_dict_005Hz_post[recID][flavour] for recID in linelength_dict_005Hz_pre.keys()])

flavour = 'pulse number'
n_pulses = np.array([linelength_dict_005Hz_pre[recID][flavour] for recID in linelength_dict_005Hz_pre.keys()])


os.chdir('C:\EAfiles\RESULTS')

df_sr=pd.DataFrame(columns=
    ['recID', 
    'mean pre-pulse [mv/s]',
    'sem pre-pulse [mv/s]',
    'mean post-pulse [mv/s]',
    'sem post-pulse [mv/s]',
    'pulse number'])

df_sr['recID'] = recIDs
df_sr['mean pre-pulse [mv/s]'] = mean_pre_linelengths 
df_sr['sem pre-pulse [mv/s]'] = sem_pre_linelengths 
df_sr['mean post-pulse [mv/s]'] = mean_post_linelengths 
df_sr['sem post-pulse [mv/s]'] = sem_post_linelengths 
df_sr['pulse number'] = n_pulses

df_sr.to_excel('Test_snipmat-005Hz-linelengths-control-PK43.xlsx')