data_overview_2.py 15 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385
  1. # -*- coding: utf-8 -*-
  2. """
  3. Code for generating the second data figure in the manuscript.
  4. Authors: Julia Sprenger, Lyuba Zehl, Michael Denker
  5. Copyright (c) 2017, Institute of Neuroscience and Medicine (INM-6),
  6. Forschungszentrum Juelich, Germany
  7. All rights reserved.
  8. Redistribution and use in source and binary forms, with or without
  9. modification, are permitted provided that the following conditions are met:
  10. * Redistributions of source code must retain the above copyright notice, this
  11. list of conditions and the following disclaimer.
  12. * Redistributions in binary form must reproduce the above copyright notice,
  13. this list of conditions and the following disclaimer in the documentation
  14. and/or other materials provided with the distribution.
  15. * Neither the names of the copyright holders nor the names of the contributors
  16. may be used to endorse or promote products derived from this software without
  17. specific prior written permission.
  18. THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
  19. ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
  20. WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
  21. DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
  22. FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
  23. DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
  24. SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
  25. CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
  26. OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
  27. OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  28. """
  29. import os
  30. import matplotlib.pyplot as plt
  31. from matplotlib import gridspec, transforms
  32. import numpy as np
  33. import quantities as pq
  34. from neo import SpikeTrain
  35. from neo.utils import *
  36. from reachgraspio import reachgraspio
  37. from neo_utils import load_segment
  38. # =============================================================================
  39. # Define data and metadata directories and general settings
  40. # =============================================================================
  41. def get_monkey_datafile(monkey):
  42. if monkey == "Lilou":
  43. return "l101210-001" # ns2 (behavior) and ns5 present
  44. elif monkey == "Nikos2":
  45. return "i140703-001" # ns2 and ns6 present
  46. else:
  47. return ""
  48. # Enter your dataset directory here
  49. datasetdir = os.path.join('..', 'datasets')
  50. chosen_els = {'Lilou': range(3, 97, 7), 'Nikos2': range(1, 97, 7)}
  51. chosen_el = {
  52. 'Lilou': chosen_els['Lilou'][0],
  53. 'Nikos2': chosen_els['Nikos2'][0]}
  54. chosen_unit = 1
  55. trial_indexes = range(14)
  56. trial_index = trial_indexes[0]
  57. chosen_events = ['TS-ON', 'WS-ON', 'CUE-ON', 'CUE-OFF', 'GO-ON', 'SR-ON',
  58. 'RW-ON', 'WS-OFF'] # , 'RW-OFF'
  59. # =============================================================================
  60. # Load data and metadata for a monkey
  61. # =============================================================================
  62. # CHANGE this parameter to load data of the different monkeys
  63. monkey = 'Nikos2'
  64. # monkey = 'Lilou'
  65. datafile = get_monkey_datafile(monkey)
  66. session = reachgraspio.ReachGraspIO(
  67. filename=os.path.join(datasetdir, datafile),
  68. odml_directory=datasetdir,
  69. verbose=False)
  70. bl = session.read_block(lazy=True)
  71. seg = bl.segments[0]
  72. # get start and stop events of trials
  73. start_events = get_events(
  74. seg, **{
  75. 'name': 'TrialEvents',
  76. 'trial_event_labels': 'TS-ON',
  77. 'performance_in_trial': session.performance_codes['correct_trial']})
  78. stop_events = get_events(
  79. seg, **{
  80. 'name': 'TrialEvents',
  81. 'trial_event_labels': 'RW-ON',
  82. 'performance_in_trial': session.performance_codes['correct_trial']})
  83. # there should only be one event object for these conditions
  84. assert len(start_events) == 1
  85. assert len(stop_events) == 1
  86. # insert epochs between 10ms before TS to 50ms after RW corresponding to trails
  87. ep = add_epoch(seg,
  88. start_events[0],
  89. stop_events[0],
  90. pre=-250 * pq.ms,
  91. post=500 * pq.ms,
  92. segment_type='complete_trials')
  93. ep.array_annotate(trialtype=start_events[0].array_annotations['belongs_to_trialtype'])
  94. # access single epoch of this data_segment
  95. epochs = get_epochs(seg, **{'segment_type': 'complete_trials'})
  96. assert len(epochs) == 1
  97. # remove spiketrains not belonging to chosen_electrode
  98. seg.spiketrains = seg.filter(targdict={'unit_id': chosen_unit},
  99. recursive=True, objects='SpikeTrainProxy')
  100. # remove all non-neural signals
  101. seg.analogsignals = seg.filter(targdict={'neural_signal': True},
  102. objects='AnalogSignalProxy')
  103. # use prefiltered data if multiple versions are present
  104. raw_signal = seg.analogsignals[0]
  105. for sig in seg.analogsignals:
  106. if sig.sampling_rate < raw_signal.sampling_rate:
  107. raw_signal = sig
  108. seg.analogsignals = [raw_signal]
  109. # replacing the segment with a new segment containing all data
  110. # to speed up cutting of segments
  111. seg = load_segment(seg, load_wavefroms=True)
  112. # only keep the chosen electrode signal in the AnalogSignal object
  113. mask = np.isin(np.asarray(seg.analogsignals[0].array_annotations['channel_ids'], dtype=int),
  114. chosen_els[monkey])
  115. seg.analogsignals[0] = seg.analogsignals[0][:, mask]
  116. # cut segments according to inserted 'complete_trials' epochs and reset trial
  117. # times
  118. cut_segments = cut_segment_by_epoch(seg, epochs[0], reset_time=True)
  119. # explicitly adding trial type annotations to cut segments
  120. for i, cut_seg in enumerate(cut_segments):
  121. cut_seg.annotate(trialtype=epochs[0].array_annotations['trialtype'][i])
  122. # =============================================================================
  123. # Define figure and subplot axis for first data overview
  124. # =============================================================================
  125. fig = plt.figure(facecolor='w')
  126. fig.set_size_inches(7.0, 9.9) # (w, h) in inches
  127. # #(7.0, 9.9) corresponds to A4 portrait ratio
  128. gs = gridspec.GridSpec(
  129. nrows=2,
  130. ncols=2,
  131. left=0.1,
  132. bottom=0.05,
  133. right=0.9,
  134. top=0.975,
  135. wspace=0.1,
  136. hspace=0.1,
  137. width_ratios=None,
  138. height_ratios=[2, 1])
  139. ax1 = plt.subplot(gs[0, 0]) # top left
  140. ax2 = plt.subplot(gs[0, 1], sharex=ax1) # top right
  141. ax3 = plt.subplot(gs[1, 0], sharex=ax1) # bottom left
  142. ax4 = plt.subplot(gs[1, 1], sharex=ax1) # bottom right
  143. fontdict_titles = {'fontsize': 9, 'fontweight': 'bold'}
  144. fontdict_axis = {'fontsize': 10, 'fontweight': 'bold'}
  145. # the x coords of the event labels are data, and the y coord are axes
  146. event_label_transform = transforms.blended_transform_factory(ax1.transData,
  147. ax1.transAxes)
  148. trialtype_colors = {
  149. 'SGHF': 'MediumBlue', 'SGLF': 'Turquoise',
  150. 'PGHF': 'DarkGreen', 'PGLF': 'YellowGreen',
  151. 'LFSG': 'Orange', 'LFPG': 'Yellow',
  152. 'HFSG': 'DarkRed', 'HFPG': 'OrangeRed',
  153. 'SGSG': 'SteelBlue', 'PGPG': 'LimeGreen',
  154. None: 'black'}
  155. event_colors = {
  156. 'TS-ON': 'indigo', 'TS-OFF': 'indigo',
  157. 'WS-ON': 'purple', 'WS-OFF': 'purple',
  158. 'CUE-ON': 'crimson', 'CUE-OFF': 'crimson',
  159. 'GO-ON': 'orangered', 'GO-OFF': 'orangered',
  160. 'SR-ON': 'darkorange',
  161. 'RW-ON': 'orange', 'RW-OFF': 'orange'}
  162. electrode_cmap = plt.get_cmap('bone')
  163. electrode_colors = [electrode_cmap(x) for x in
  164. np.tile(np.array([0.3, 0.7]), int(len(chosen_els[monkey]) / 2))]
  165. time_unit = 'ms'
  166. lfp_unit = 'uV'
  167. # define scaling factors for analogsignals
  168. anasig_std = np.mean(np.std(cut_segments[trial_index].analogsignals[0].rescale(lfp_unit), axis=0))
  169. anasig_offset = 3 * anasig_std
  170. # =============================================================================
  171. # SUPPLEMENTARY PLOTTING functions
  172. # =============================================================================
  173. def add_scalebar(ax, std):
  174. # the x coords of the scale bar are axis, and the y coord are data
  175. scalebar_transform = transforms.blended_transform_factory(ax.transAxes,
  176. ax.transData)
  177. # adding scalebar
  178. yscalebar = max(int(std.rescale(lfp_unit)), 1) * getattr(pq, lfp_unit) * 2
  179. scalebar_offset = -2 * std
  180. ax.vlines(x=0.4,
  181. ymin=(scalebar_offset - yscalebar).magnitude,
  182. ymax=scalebar_offset.magnitude,
  183. color='k',
  184. linewidth=4,
  185. transform=scalebar_transform)
  186. ax.text(0.4, (scalebar_offset - 0.5 * yscalebar).magnitude,
  187. ' %i %s' % (yscalebar.magnitude, lfp_unit),
  188. ha="left", va="center", rotation=0, color='k',
  189. size=8, transform=scalebar_transform)
  190. # =============================================================================
  191. # PLOT DATA OF SINGLE TRIAL (left plots)
  192. # =============================================================================
  193. # get data of selected trial
  194. selected_trial = cut_segments[trial_index]
  195. # PLOT DATA FOR EACH CHOSEN ELECTRODE
  196. for el_idx, electrode_id in enumerate(chosen_els[monkey]):
  197. # PLOT ANALOGSIGNALS in upper plot
  198. chids = np.asarray(cut_segments[0].analogsignals[0].array_annotations['channel_ids'], dtype=int)
  199. chosen_el_idx = np.where(chids == electrode_id)[0][0]
  200. anasig = selected_trial.analogsignals[0][:, chosen_el_idx]
  201. ax1.plot(anasig.times.rescale(time_unit),
  202. np.asarray(anasig.rescale(lfp_unit))
  203. + anasig_offset.magnitude * el_idx,
  204. color=electrode_colors[el_idx])
  205. # PLOT SPIKETRAINS in lower plot
  206. spiketrains = selected_trial.filter(
  207. channel_id=electrode_id, objects=SpikeTrain)
  208. for spiketrain in spiketrains:
  209. ax3.plot(spiketrain.times.rescale(time_unit),
  210. np.zeros(len(spiketrain.times)) + el_idx, 'k|')
  211. # PLOT EVENTS in both plots
  212. for event_type in chosen_events:
  213. # get events of each chosen event type
  214. event_data = get_events(selected_trial,
  215. **{'trial_event_labels': event_type})
  216. for event in event_data:
  217. event_color = event_colors[event.array_annotations['trial_event_labels'][0]]
  218. # adding lines
  219. for ax in [ax1, ax3]:
  220. ax.axvline(event.times.rescale(time_unit),
  221. color=event_color,
  222. zorder=0.5)
  223. # adding labels
  224. ax1.text(event.times.rescale(time_unit), 0,
  225. event.array_annotations['trial_event_labels'][0],
  226. ha="center", va="top", rotation=45, color=event_color,
  227. size=8, transform=event_label_transform)
  228. # SUBPLOT ADJUSTMENTS
  229. ax1.set_title('single trial', fontdict=fontdict_titles)
  230. ax1.set_ylabel('electrode id', fontdict=fontdict_axis)
  231. ax1.set_yticks(np.arange(len(chosen_els[monkey])) * anasig_offset)
  232. ax1.set_yticklabels(chosen_els[monkey])
  233. ax1.autoscale(enable=True, axis='y')
  234. plt.setp(ax1.get_xticklabels(), visible=False) # show no xticklabels
  235. ax3.set_ylabel('electrode id', fontdict=fontdict_axis)
  236. ax3.set_yticks(range(0, len(chosen_els[monkey])))
  237. ax3.set_yticklabels(np.asarray(chosen_els[monkey]))
  238. ax3.set_ylim(-1, len(chosen_els[monkey]))
  239. ax3.set_xlabel('time [%s]' % time_unit, fontdict=fontdict_axis)
  240. # ax3.autoscale(axis='y')
  241. # =============================================================================
  242. # PLOT DATA OF SINGLE ELECTRODE
  243. # =============================================================================
  244. # plot data for each chosen trial
  245. chids = np.asarray(cut_segments[0].analogsignals[0].array_annotations['channel_ids'], dtype=int)
  246. chosen_el_idx = np.where(chids == chosen_el[monkey])[0][0]
  247. for trial_idx, trial_id in enumerate(trial_indexes):
  248. trial_spikes = cut_segments[trial_id].filter(channel_id=chosen_el[monkey], objects='SpikeTrain')
  249. trial_type = cut_segments[trial_id].annotations['trialtype']
  250. trial_color = trialtype_colors[trial_type]
  251. t_signal = cut_segments[trial_id].analogsignals[0][:, chosen_el_idx]
  252. # PLOT ANALOGSIGNALS in upper plot
  253. ax2.plot(t_signal.times.rescale(time_unit),
  254. np.asarray(t_signal.rescale(lfp_unit))
  255. + anasig_offset.magnitude * trial_idx,
  256. color=trial_color, zorder=1)
  257. for t_data in trial_spikes:
  258. # PLOT SPIKETRAINS in lower plot
  259. ax4.plot(t_data.times.rescale(time_unit),
  260. np.ones(len(t_data.times)) + trial_idx, 'k|')
  261. # PLOT EVENTS in both plots
  262. for event_type in chosen_events:
  263. # get events of each chosen event type
  264. event_data = get_events(cut_segments[trial_id], **{'trial_event_labels': event_type})
  265. for event in event_data:
  266. color = event_colors[event.array_annotations['trial_event_labels'][0]]
  267. ax2.vlines(x=event.times.rescale(time_unit),
  268. ymin=(trial_idx - 0.5) * anasig_offset,
  269. ymax=(trial_idx + 0.5) * anasig_offset,
  270. color=color,
  271. zorder=2)
  272. ax4.vlines(x=event.times.rescale(time_unit),
  273. ymin=trial_idx + 1 - 0.4,
  274. ymax=trial_idx + 1 + 0.4,
  275. color=color,
  276. zorder=0.5)
  277. # SUBPLOT ADJUSTMENTS
  278. ax2.set_title('single electrode', fontdict=fontdict_titles)
  279. ax2.set_ylabel('trial id', fontdict=fontdict_axis)
  280. ax2.set_yticks(np.asarray(trial_indexes) * anasig_offset)
  281. ax2.set_yticklabels(
  282. [epochs[0].array_annotations['trial_id'][_] for _ in trial_indexes])
  283. ax2.yaxis.set_label_position("right")
  284. ax2.tick_params(direction='in', length=3, labelleft='off', labelright='on')
  285. ax2.autoscale(enable=True, axis='y')
  286. add_scalebar(ax2, anasig_std)
  287. plt.setp(ax2.get_xticklabels(), visible=False) # show no xticklabels
  288. ax4.set_ylabel('trial id', fontdict=fontdict_axis)
  289. ax4.set_xlabel('time [%s]' % time_unit, fontdict=fontdict_axis)
  290. start, end = ax4.get_xlim()
  291. ax4.xaxis.set_ticks(np.arange(start, end, 1000))
  292. ax4.xaxis.set_ticks(np.arange(start, end, 500), minor=True)
  293. ax4.set_yticks(range(1, len(trial_indexes) + 1))
  294. ax4.set_yticklabels(np.asarray(
  295. [epochs[0].array_annotations['trial_id'][_] for _ in trial_indexes]))
  296. ax4.yaxis.set_label_position("right")
  297. ax4.tick_params(direction='in', length=3, labelleft='off', labelright='on')
  298. ax4.autoscale(enable=True, axis='y')
  299. # GENERAL PLOT ADJUSTMENTS
  300. # adjust font sizes of ticks
  301. for ax in [ax4.yaxis, ax4.xaxis, ax3.xaxis, ax3.yaxis]:
  302. for tick in ax.get_major_ticks():
  303. tick.label.set_fontsize(10)
  304. # adjust time range on x axis
  305. t_min = np.min([cut_segments[tid].t_start.rescale(time_unit)
  306. for tid in trial_indexes])
  307. t_max = np.max([cut_segments[tid].t_stop.rescale(time_unit)
  308. for tid in trial_indexes])
  309. ax1.set_xlim(t_min, t_max)
  310. add_scalebar(ax1, anasig_std)
  311. # =============================================================================
  312. # SAVE FIGURE
  313. # =============================================================================
  314. fname = 'data_overview_2_%s' % monkey
  315. for file_format in ['eps', 'pdf', 'png']:
  316. fig.savefig(fname + '.%s' % file_format, dpi=400, format=file_format)