1
1

example.py 12 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280
  1. # -*- coding: utf-8 -*-
  2. """
  3. Example code for loading and processing of a recording of the reach-
  4. to-grasp experiments conducted at the Institute de Neurosciences de la Timone
  5. by Thomas Brochier and Alexa Riehle.
  6. Authors: Julia Sprenger, Lyuba Zehl, Michael Denker
  7. Copyright (c) 2017, Institute of Neuroscience and Medicine (INM-6),
  8. Forschungszentrum Juelich, Germany
  9. All rights reserved.
  10. Redistribution and use in source and binary forms, with or without
  11. modification, are permitted provided that the following conditions are met:
  12. * Redistributions of source code must retain the above copyright notice, this
  13. list of conditions and the following disclaimer.
  14. * Redistributions in binary form must reproduce the above copyright notice,
  15. this list of conditions and the following disclaimer in the documentation
  16. and/or other materials provided with the distribution.
  17. * Neither the names of the copyright holders nor the names of the contributors
  18. may be used to endorse or promote products derived from this software without
  19. specific prior written permission.
  20. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
  21. ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
  22. WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
  23. DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
  24. FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  25. DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
  26. SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
  27. CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
  28. OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
  29. OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  30. """
  31. # This loads the Neo and odML libraries shipped with this code. For production
  32. # use, please use the newest releases of odML and Neo.
  33. import load_local_neo_odml_elephant
  34. import os
  35. import numpy as np
  36. import matplotlib.pyplot as plt
  37. import quantities as pq
  38. from neo import Block, Segment
  39. from elephant.signal_processing import butter
  40. from reachgraspio import reachgraspio
  41. from neo.utils import cut_segment_by_epoch, add_epoch, get_events
  42. from neo_utils import load_segment
  43. # =============================================================================
  44. # Load data
  45. #
  46. # As a first step, we partially load the data file into memory as a Neo object.
  47. # =============================================================================
  48. # Specify the path to the recording session to load, eg,
  49. # '/home/user/l101210-001'
  50. session_name = os.path.join('..', 'datasets', 'i140703-001')
  51. # session_name = os.path.join('..', 'datasets', 'l101210-001')
  52. odml_dir = os.path.join('..', 'datasets')
  53. # Open the session for reading
  54. session = reachgraspio.ReachGraspIO(session_name, odml_directory=odml_dir)
  55. # Read a the complete dataset in lazy mode generating all neo objects,
  56. # but not loading data into memory. The lazy neo structure will contain objects
  57. # to capture all recorded data types (time series at 1000Hz (ns2) and 30kHz (ns6)
  58. # scaled to units of voltage, sorted spike trains, spike waveforms and events)
  59. # of the recording session and return it as a Neo Block. The
  60. # time shift of the ns2 signal (LFP) induced by the online filter is
  61. # automatically corrected for by a heuristic factor stored in the metadata
  62. # (correct_filter_shifts=True).
  63. block = session.read_block(lazy=True, correct_filter_shifts=True)
  64. # Validate there is only a single Segment present in the block
  65. assert len(block.segments) == 1
  66. # loading data content of all data objects during the first 300 seconds
  67. data_segment = load_segment(block.segments[0], time_range=(None, 300*pq.s))
  68. # =============================================================================
  69. # Create offline filtered LFP
  70. #
  71. # Here, we construct one offline filtered LFP from each ns5 (monkey L) or ns6
  72. # (monkey N) raw recording trace. For monkey N, this filtered LFP can be
  73. # compared to the LFPs in the ns2 file (note that monkey L contains only
  74. # behavioral signals in the ns2 file). Also, we assign telling names to each
  75. # Neo AnalogSignal, which is used for plotting later on in this script.
  76. # =============================================================================
  77. target_channel_id = 62
  78. nsx_to_anasig_name = {2: 'LFP signal (online filtered)',
  79. 5: 'raw signal',
  80. 6: 'raw signal'}
  81. filtered_anasig = None
  82. raw_anasig = None
  83. # identify neuronal signals and provide labels for plotting
  84. for anasig in data_segment.analogsignals:
  85. # skip non-neuronal signals
  86. if not anasig.annotations['neural_signal']:
  87. continue
  88. # identify nsx source of signals in this AnalogSignal object
  89. nsx = np.unique(anasig.array_annotations['nsx'])
  90. assert len(nsx) == 1, 'Different nsx sources in AnalogSignal'
  91. nsx = nsx[0]
  92. if nsx == 2:
  93. # AnalogSignal is LFP from ns2
  94. anasig.name = f'LFP (online filter, ns2)'
  95. filtered_anasig = anasig
  96. elif nsx in [5, 6]:
  97. # AnalogSignal is raw signal from ns5 or ns6
  98. anasig.name = f'raw (ns{nsx})'
  99. raw_anasig = anasig
  100. # Create LFP signal by filtering raw signal if not present already
  101. if filtered_anasig is None:
  102. # Use the Elephant library to filter the signal, filter only target channel
  103. target_channel_index = np.where(target_channel_id == raw_anasig.array_annotations['channel_ids'])[0]
  104. raw_signal = raw_anasig[:, target_channel_index]
  105. f_anasig = butter(
  106. raw_signal,
  107. highpass_freq=None,
  108. lowpass_freq=250 * pq.Hz,
  109. order=4)
  110. f_anasig.name = f'LFP (offline filtered ns{anasig.array_annotations["nsx"][0]})'
  111. # Attach all offline filtered LFPs to the segment of data
  112. data_segment.analogsignals.append(f_anasig)
  113. # =============================================================================
  114. # Construct analysis epochs
  115. #
  116. # In this step we extract and cut the data into time segments (termed analysis
  117. # epochs) that we wish to analyze. We contrast these analysis epochs to the
  118. # behavioral trials that are defined by the experiment as occurrence of a Trial
  119. # Start (TS-ON) event in the experiment. Concretely, here our analysis epochs
  120. # are constructed as a cutout of 25ms of data around the TS-ON event of all
  121. # successful behavioral trials.
  122. # =============================================================================
  123. # Get Trial Start (TS-ON) events of all successful behavioral trials
  124. # (corresponds to performance code 255, which is accessed for convenience and
  125. # better legibility in the dictionary attribute performance_codes of the
  126. # ReachGraspIO class).
  127. #
  128. # To this end, we filter all event objects of the loaded data to match the name
  129. # "TrialEvents", which is the Event object containing all Events available (see
  130. # documentation of ReachGraspIO). From this Event object we extract only events
  131. # matching "TS-ON" and the desired trial performance code (which are
  132. # annotations of the Event object).
  133. start_events = get_events(
  134. data_segment,
  135. name='TrialEvents',
  136. trial_event_labels='TS-ON',
  137. performance_in_trial=session.performance_codes['correct_trial'])
  138. print('got start events.')
  139. # Extract single Neo Event object containing all TS-ON triggers
  140. assert len(start_events) == 1
  141. start_event = start_events[0]
  142. # Construct analysis epochs from 10ms before the TS-ON of a successful
  143. # behavioral trial to 15ms after TS-ON. The name "analysis_epochs" is given to
  144. # the resulting Neo Epoch object. The object is not attached to the Neo
  145. # Segment. The parameter event2 of add_epoch() is left empty, since we are
  146. # cutting around a single event, as opposed to cutting between two events.
  147. pre = -10 * pq.ms
  148. post = 15 * pq.ms
  149. epoch = add_epoch(
  150. data_segment,
  151. event1=start_event, event2=None,
  152. pre=pre, post=post,
  153. attach_result=False,
  154. name='analysis_epochs',
  155. array_annotations=start_event.array_annotations)
  156. print('added epoch.')
  157. # Create new segments of data cut according to the analysis epochs of the
  158. # 'analysis_epochs' Neo Epoch object. The time axes of all segments are aligned
  159. # such that each segment starts at time 0 (parameter reset_times); annotations
  160. # describing the analysis epoch are carried over to the segments. A new Neo
  161. # Block named "data_cut_to_analysis_epochs" is created to capture all cut
  162. # analysis epochs. For execution time reason, we are only considering the
  163. # first 10 epochs here.
  164. cut_trial_block = Block(name="data_cut_to_analysis_epochs")
  165. cut_trial_block.segments = cut_segment_by_epoch(
  166. data_segment, epoch[:10], reset_time=True)
  167. # =============================================================================
  168. # Plot data
  169. # =============================================================================
  170. # Determine the first existing trial ID i from the Event object containing all
  171. # start events. Then, by calling the filter() function of the Neo Block
  172. # "data_cut_to_analysis_epochs" containing the data cut into the analysis
  173. # epochs, we ask to return all Segments annotated by the behavioral trial ID i.
  174. # In this case this call should return one matching analysis epoch around TS-ON
  175. # belonging to behavioral trial ID i. For monkey N, this is trial ID 1, for
  176. # monkey L this is trial ID 2 since trial ID 1 is not a correct trial.
  177. trial_id = int(np.min(start_event.array_annotations['trial_id']))
  178. trial_segments = cut_trial_block.filter(
  179. targdict={"trial_id": trial_id}, objects=Segment)
  180. assert len(trial_segments) == 1
  181. trial_segment = trial_segments[0]
  182. # Create figure
  183. fig = plt.figure(facecolor='w')
  184. time_unit = pq.CompoundUnit('1./30000*s')
  185. amplitude_unit = pq.microvolt
  186. nsx_colors = {2: 'k', 5: 'r', 6: 'b'}
  187. # Loop through all AnalogSignal objects and plot the signal of the target channel
  188. # in a color corresponding to its sampling frequency (i.e., originating from the ns2/ns5 or ns2/ns6).
  189. for i, anasig in enumerate(trial_segment.analogsignals):
  190. # only visualize neural data
  191. if anasig.annotations['neural_signal']:
  192. nsx = anasig.array_annotations['nsx'][0]
  193. target_channel_index = np.where(anasig.array_annotations['channel_ids'] == target_channel_id)[0]
  194. target_signal = anasig[:, target_channel_index]
  195. plt.plot(
  196. target_signal.times.rescale(time_unit),
  197. target_signal.squeeze().rescale(amplitude_unit),
  198. label=target_signal.name,
  199. color=nsx_colors[nsx])
  200. # Loop through all spike trains and plot the spike time, and overlapping the
  201. # wave form of the spike used for spike sorting stored separately in the nev
  202. # file.
  203. for st in trial_segment.spiketrains:
  204. color = np.random.rand(3,)
  205. if st.annotations['channel_id'] == target_channel_id:
  206. for spike_id, spike in enumerate(st):
  207. # Plot spike times
  208. plt.axvline(
  209. spike.rescale(time_unit).magnitude,
  210. color=color,
  211. label='Unit ID %i' % st.annotations['unit_id'])
  212. # Plot waveforms
  213. waveform = st.waveforms[spike_id, 0, :]
  214. waveform_times = np.arange(len(waveform))*time_unit + spike
  215. plt.plot(
  216. waveform_times.rescale(time_unit).magnitude,
  217. waveform.rescale(amplitude_unit),
  218. '--',
  219. linewidth=2,
  220. color=color,
  221. zorder=0)
  222. # Loop through all events
  223. for event in trial_segment.events:
  224. if event.name == 'TrialEvents':
  225. for ev_id, ev in enumerate(event):
  226. plt.axvline(
  227. ev.rescale(time_unit),
  228. alpha=0.2,
  229. linewidth=3,
  230. linestyle='dashed',
  231. label=f'event {event.array_annotations["trial_event_labels"][ev_id]}')
  232. # Finishing touches on the plot
  233. plt.autoscale(enable=True, axis='x', tight=True)
  234. plt.xlabel(time_unit.name)
  235. plt.ylabel(amplitude_unit.name)
  236. plt.legend(loc=4, fontsize=10)
  237. # Save plot
  238. fname = 'example_plot'
  239. for file_format in ['eps', 'png', 'pdf']:
  240. fig.savefig(fname + '.%s' % file_format, dpi=400, format=file_format)