'script2_linelengths_of_snippets_Fig2.py' ändern

Linelengths_of_snippets_Fig2.py

@@ -1,270 +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')
-
+# -*- 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')
+