Linelengths_of_snippets_Fig2.py 10 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270
  1. # -*- coding: utf-8 -*-
  2. """
  3. Created on Wed Oct 14 11:18:44 2020
  4. @author: kleisp
  5. This is a script for getting matrices of snippets around the stimulation pulses in PACK-expressing mice,
  6. 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)
  7. 0.1 Hz and 0.05 Hz need to be treated separately
  8. need to load dictionaries containing recIDs with corresponding smrs for using raw data
  9. dict- dictionary
  10. stim- light pulse
  11. pre- before light pulse
  12. post-after light pulse
  13. snipmat- matrix containing snippets
  14. """
  15. import os
  16. os.chdir('C:\EAfiles\MyCode')
  17. from glob import glob
  18. import numpy as np
  19. import pickle
  20. import pandas as pd
  21. from scipy import stats
  22. os.chdir('C:\EAfiles\EAdetection')
  23. import core.ea_management as eam
  24. reload(eam)
  25. #%% define functions to extract snippets before and after the stim
  26. ''' 'pre' or 'post' are the snippet length in seconds, here 2 is used, sr=sampling rate (here 500)'''
  27. # get the spike trace before stimulation (pre):
  28. #stimtime is one element from stimtimes array
  29. def get_spiketrace_pre(stimtimes, datatrace, pre, sr):
  30. cutout_pre = np.int(pre*sr)
  31. stim_index= np.int(stimtimes*sr)
  32. spiketrace = datatrace[stim_index-cutout_pre:stim_index]
  33. return spiketrace
  34. # get the spike trace after(post) stimulation:
  35. def get_spiketrace_post(stimtime, datatrace, post, sr):
  36. cutout_post = np.int(post*sr)
  37. stim_index= np.int(stimtime*sr)
  38. spiketrace = datatrace[stim_index:stim_index+cutout_post]
  39. return spiketrace
  40. #get snippet before each stimulation and insert it into a matrix, called snipmat:
  41. def get_snipmat_pre(datatrace, recID, pre, stimtimes, sr):
  42. cutout_pre = np.int(pre*sr)
  43. N_stim = len(stimtimes) #number of stimuli
  44. snippts = cutout_pre
  45. snipmat = np.zeros((N_stim,snippts))
  46. for tt,stimtime in enumerate(stimtimes):
  47. spiketrace = get_spiketrace_pre(stimtime, datatrace, pre, sr) # instead of datatrace you can use zTrace
  48. snipmat[tt] = spiketrace
  49. return snipmat
  50. #get snippet after each stimulation and insert it into a matrix, called snipmat:
  51. def get_snipmat_post(datatrace, recID, post, stimtimes, sr):
  52. cutout_post = np.int(post*sr)
  53. N_stim = len(stimtimes) #number of stimuli
  54. snippts = cutout_post
  55. snipmat = np.zeros((N_stim,snippts))
  56. for tt,stimtime in enumerate(stimtimes):
  57. spiketrace = get_spiketrace_post(stimtime, datatrace, post, sr) # stimtime is just one timepoint from stimtimes
  58. snipmat[tt] = spiketrace
  59. return snipmat
  60. #%% defining the function to get line length of each snippet(snipmat line length)
  61. def get_snipmat_linelength(snipmat_pre, snipmat_post, linelength_dict_pre, linelength_dict_post, recID):
  62. linelength_list=[]
  63. linelength_s_list=[]
  64. rms_list=[]
  65. count=0
  66. for i in xrange(snipmat_pre.shape[0]):
  67. data=snipmat_pre[i]
  68. N_datapoints=len(data)
  69. duration=N_datapoints/500 #sr=500
  70. linelength = np.sum(np.abs(np.diff(data)))
  71. linelength_list.append(linelength)
  72. linelength_s=linelength/duration
  73. linelength_s_list.append(linelength_s)
  74. count=count+1
  75. linelength_dict_pre[recID]={}
  76. linelength_dict_pre[recID]['sem pre-pulse linelength/second']=stats.sem(linelength_s_list)
  77. linelength_dict_pre[recID]['mean pre-pulse linelength/second']=np.mean(linelength_s_list)
  78. linelength_dict_pre[recID]['mean pre-pulse linelength']=np.mean(linelength_list)
  79. linelength_dict_pre[recID]['pulse number']=count
  80. linelength_list=[]
  81. linelength_s_list=[]
  82. rms_list=[]
  83. count=0
  84. for i in xrange(snipmat_post.shape[0]):
  85. data=snipmat_post[i]
  86. N_datapoints=len(data)
  87. duration=N_datapoints/500
  88. linelength = np.sum(np.abs(np.diff(data)))
  89. linelength_list.append(linelength)
  90. linelength_s=linelength/duration
  91. linelength_s_list.append(linelength_s)
  92. count=count+1
  93. linelength_dict_post[recID]={}
  94. linelength_dict_post[recID]['sem post-pulse linelength/second']=stats.sem(linelength_s_list)
  95. linelength_dict_post[recID]['mean post-pulse linelength/second']=np.mean(linelength_s_list)
  96. linelength_dict_post[recID]['mean post-pulse linelength']=np.mean(linelength_list)
  97. linelength_dict_post[recID]['pulse number']=count
  98. #%% in order to exclude unwanted recordings
  99. pathtobadrecIDs = 'C:\EAfiles\ID_Dateien\\badrecIDs.txt'
  100. badrecIDs = []
  101. with open (pathtobadrecIDs,'r') as f:
  102. badIDs = f.read().splitlines()
  103. badrecIDs = badrecIDs + badIDs
  104. #%% Example of using the script to extract the snipmats
  105. #Saline PACK mice 0.05Hz: get snipmat dicts (for before and after 0.05 Hz pulses) (downsampled data)
  106. os.chdir('C:\EAfiles\EAdetection')
  107. ymlfiles = glob('C:\EAfiles\yml_setup_PK43\*.yml') #takes all yml files from this folder
  108. ymlfile = ymlfiles
  109. lines = np.arange(0, len(ymlfile)) #as many lines as there are ymlfiles
  110. snipmat_005Hz_pre_controls={}
  111. snipmat_005Hz_post_controls={}
  112. sr= 500
  113. pre, post = 2, 2 #2 seconds before and after pulse
  114. for i in lines:
  115. recID = os.path.basename(ymlfile[i]).split('__')[0] #getting the IDs
  116. print (i)
  117. print (recID)
  118. if recID in badrecIDs:
  119. pass
  120. else:
  121. excelpath = 'C:\EAfiles\Excel\\startstop_controls_HCc.xlsx'
  122. df = pd.read_excel(excelpath) #reading an excel file with start and stop times of each recording
  123. df.set_index('recID', inplace=True) #recID
  124. stimstop = int(df.stop[recID]) #stop time
  125. if recID=='PK16_005Hz1_HCc': #adding an exception (can skip this step if you have no exceptions)
  126. stimstart=int(df.start[recID])
  127. else:
  128. stimstart=10 #the stim usually starts at 10s
  129. stimstep=20 #speriod is 20 s
  130. stimtimes=np.arange(stimstart,stimstop,stimstep)
  131. aRec = eam.Rec(init_ymlpath=ymlfile[i])
  132. snipmat_005Hz_pre_controls[recID]={}
  133. data=aRec.raw_data
  134. #z_data=(data-np.mean(data))/np.std(data)
  135. snipmat_pre=get_snipmat_pre(data, recID, pre, stimtimes, sr=500)
  136. snipmat_005Hz_pre_controls[recID]=snipmat_pre
  137. snipmat_005Hz_post_controls[recID]={}
  138. snipmat_post=get_snipmat_post(data, recID, post, stimtimes, sr=500)
  139. snipmat_005Hz_post_controls[recID]=snipmat_post
  140. # save snipmat_dict
  141. filename = 'C:\EAfiles\MyDicts\snipmat_005Hz_pre_controls'
  142. outfile = open(filename, 'wb') #write binary: mit Erlaubnis zu (über)schreiben
  143. pickle.dump(snipmat_005Hz_pre_controls, outfile)
  144. outfile.close()
  145. filename = 'C:\EAfiles\MyDicts\snipmat_005Hz_post_controls'
  146. outfile = open(filename, 'wb') #write binary: mit Erlaubnis zu (über)schreiben
  147. pickle.dump(snipmat_005Hz_post_controls, outfile)
  148. outfile.close()
  149. #%% Example of using the script to create dictionaries with line length before and after 0.05Hz pulses from snipmats above
  150. os.chdir('C:\EAfiles\EAdetection')
  151. ymlfiles = glob('C:\EAfiles\yml_setup_PK43\*.yml') #takes all yml files from this folder
  152. ymlfile = ymlfiles
  153. lines = np.arange(0, len(ymlfile)) #as many lines as there are ymlfiles
  154. linelength_dict_005Hz_pre={}
  155. linelength_dict_005Hz_post={}
  156. filename = 'C:\EAfiles\MyDicts\snipmat_005Hz_pre_controls'
  157. infile = open(filename,'rb') #read binary
  158. snipmat_005Hz_pre_controls = pickle.load(infile)
  159. infile.close()
  160. filename = 'C:\EAfiles\MyDicts\snipmat_005Hz_post_controls'
  161. infile = open(filename,'rb') #read binary
  162. snipmat_005Hz_post_controls = pickle.load(infile)
  163. infile.close()
  164. for i in lines:
  165. recID = os.path.basename(ymlfile[i]).split('__')[0] #getting the IDs
  166. print (i)
  167. print (recID)
  168. if recID in badrecIDs: #template also needs to be excluded
  169. #badrecIDs = ['EPxy_xHz_HCi1','...']
  170. pass
  171. else:
  172. snipmat_pre=snipmat_005Hz_pre_controls[recID]
  173. snipmat_post=snipmat_005Hz_post_controls[recID]
  174. get_snipmat_linelength(snipmat_pre, snipmat_post, linelength_dict_005Hz_pre, linelength_dict_005Hz_post, recID)
  175. filename = 'C:\EAfiles\MyDicts\linelength_dict_005Hz_post'
  176. outfile = open(filename, 'wb') #write binary: mit Erlaubnis zu (über)schreiben
  177. pickle.dump(linelength_dict_005Hz_post, outfile)
  178. outfile.close()
  179. filename = 'C:\EAfiles\MyDicts\linelength_dict_005Hz_pre'
  180. outfile = open(filename, 'wb') #write binary: mit Erlaubnis zu (über)schreiben
  181. pickle.dump(linelength_dict_005Hz_pre, outfile)
  182. outfile.close()
  183. #%% make an excel file of the pre- and post-pulse linelengths (mean and SEM)
  184. recIDs = np.array(linelength_dict_005Hz_pre.keys())
  185. flavour = 'mean pre-pulse linelength/second'
  186. mean_pre_linelengths = np.array([linelength_dict_005Hz_pre[recID][flavour] for recID in linelength_dict_005Hz_pre.keys()])
  187. flavour = 'sem pre-pulse linelength/second'
  188. sem_pre_linelengths = np.array([linelength_dict_005Hz_pre[recID][flavour] for recID in linelength_dict_005Hz_pre.keys()])
  189. flavour = 'mean post-pulse linelength/second'
  190. mean_post_linelengths = np.array([linelength_dict_005Hz_post[recID][flavour] for recID in linelength_dict_005Hz_pre.keys()])
  191. flavour = 'sem post-pulse linelength/second'
  192. sem_post_linelengths = np.array([linelength_dict_005Hz_post[recID][flavour] for recID in linelength_dict_005Hz_pre.keys()])
  193. flavour = 'pulse number'
  194. n_pulses = np.array([linelength_dict_005Hz_pre[recID][flavour] for recID in linelength_dict_005Hz_pre.keys()])
  195. os.chdir('C:\EAfiles\RESULTS')
  196. df_sr=pd.DataFrame(columns=
  197. ['recID',
  198. 'mean pre-pulse [mv/s]',
  199. 'sem pre-pulse [mv/s]',
  200. 'mean post-pulse [mv/s]',
  201. 'sem post-pulse [mv/s]',
  202. 'pulse number'])
  203. df_sr['recID'] = recIDs
  204. df_sr['mean pre-pulse [mv/s]'] = mean_pre_linelengths
  205. df_sr['sem pre-pulse [mv/s]'] = sem_pre_linelengths
  206. df_sr['mean post-pulse [mv/s]'] = mean_post_linelengths
  207. df_sr['sem post-pulse [mv/s]'] = sem_post_linelengths
  208. df_sr['pulse number'] = n_pulses
  209. df_sr.to_excel('Test_snipmat-005Hz-linelengths-control-PK43.xlsx')