Dateien hochladen nach ''

Linelengths_of_snippets_Fig2.py

@@ -0,0 +1,270 @@
+# -*- 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, 
+this can later be used for plotting the snippets on top of each other (in this file below) and for line length calculations (in another file)
+0.1 Hz and 0.05 Hz need to be treated separately
+need to load dictionaries containing recIDs with corresponding smrs for using raw data
+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)
+
+#%% 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   
+  
+#%% 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
+    
+
+#%% 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()
+
+#%% 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')
+