diff --git a/neo/core/segment.py b/neo/core/segment.py index 03969bec7..59634e48a 100644 --- a/neo/core/segment.py +++ b/neo/core/segment.py @@ -81,7 +81,7 @@ class Segment(Container): ('rec_datetime', datetime), ('index', int)) + Container._recommended_attrs) - _repr_pretty_containers = ('analogsignals',) + _repr_pretty_containers = ('analogsignals', 'imagesequences') def __init__(self, name=None, description=None, file_origin=None, file_datetime=None, rec_datetime=None, index=None, diff --git a/neo/io/nwbio.py b/neo/io/nwbio.py index a5a04f929..bd0c4cb4e 100644 --- a/neo/io/nwbio.py +++ b/neo/io/nwbio.py @@ -7,6 +7,11 @@ Documentation : https://www.nwb.org/ Depends on: h5py, nwb, dateutil Supported: Read, Write +Specification - https://github.com/NeurodataWithoutBorders/specification +Python APIs - (1) https://github.com/NeurodataWithoutBorders/pynwb + (2) https://github.com/AllenInstitute/nwb-api/tree/master/ainwb + (3) https://github.com/AllenInstitute/AllenSDK/blob/master/allensdk/core/nwb_data_set.py + (4) https://github.com/NeurodataWithoutBorders/api-python Python API - https://pynwb.readthedocs.io Sample datasets from CRCNS - https://crcns.org/NWB Sample datasets from Allen Institute - http://alleninstitute.github.io/AllenSDK/cell_types.html#neurodata-without-borders @@ -28,6 +33,9 @@ import numpy as np import quantities as pq +from siunits import * +import subprocess +from subprocess import run from neo.io.baseio import BaseIO from neo.io.proxyobjects import ( AnalogSignalProxy as BaseAnalogSignalProxy, @@ -35,8 +43,9 @@ EpochProxy as BaseEpochProxy, SpikeTrainProxy as BaseSpikeTrainProxy ) -from neo.core import (Segment, SpikeTrain, Epoch, Event, AnalogSignal, - IrregularlySampledSignal, Block, ImageSequence) +from neo.core import (Segment, SpikeTrain, Epoch, Event, AnalogSignal, #Unit, ChannelIndex + IrregularlySampledSignal, Block, ImageSequence, + RectangularRegionOfInterest, CircularRegionOfInterest, PolygonRegionOfInterest) # PyNWB imports try: @@ -50,14 +59,16 @@ from pynwb.image import ImageSeries from pynwb.spec import NWBAttributeSpec, NWBDatasetSpec, NWBGroupSpec, NWBNamespace, NWBNamespaceBuilder from pynwb.device import Device - # For calcium imaging data - from pynwb.ophys import TwoPhotonSeries, OpticalChannel, ImageSegmentation, Fluorescence + from pynwb.ophys import TwoPhotonSeries, OpticalChannel, ImageSegmentation, Fluorescence, ImagingPlane, PlaneSegmentation, RoiResponseSeries + from pynwb import validate, NWBHDF5IO have_pynwb = True except ImportError: have_pynwb = False except SyntaxError: # pynwb doesn't support Python 2.7 have_pynwb = False +import random + # hdmf imports try: from hdmf.spec import (LinkSpec, GroupSpec, DatasetSpec, SpecNamespace, @@ -77,7 +88,8 @@ "experiment_description", "session_id", "institution", "keywords", "notes", "pharmacology", "protocol", "related_publications", "slices", "source_script", "source_script_file_name", "data_collection", "surgery", "virus", "stimulus_notes", - "lab", "session_description" + "lab", "session_description", + "rec_datetime", ) POSSIBLE_JSON_FIELDS = ( @@ -128,6 +140,8 @@ def statistics(block): # todo: move this to be a property of Block "IrregularlySampledSignal": {"count": 0}, "Epoch": {"count": 0}, "Event": {"count": 0}, + "ImageSequence": {"count": 0}, + "Fluorescence": {"count": 0}, } for segment in block.segments: stats["SpikeTrain"]["count"] += len(segment.spiketrains) @@ -135,6 +149,8 @@ def statistics(block): # todo: move this to be a property of Block stats["IrregularlySampledSignal"]["count"] += len(segment.irregularlysampledsignals) stats["Epoch"]["count"] += len(segment.epochs) stats["Event"]["count"] += len(segment.events) + stats["ImageSequence"]["count"] += len(segment.imagesequences) + stats["Fluorescence"]["count"] += len(segment.imagesequences) return stats @@ -215,7 +231,8 @@ class NWBIO(BaseIO): Class for "reading" experimental data from a .nwb file, and "writing" a .nwb file from Neo """ supported_objects = [Block, Segment, AnalogSignal, IrregularlySampledSignal, - SpikeTrain, Epoch, Event, ImageSequence] + SpikeTrain, Epoch, Event, ImageSequence, + RectangularRegionOfInterest, CircularRegionOfInterest, PolygonRegionOfInterest] readable_objects = supported_objects writeable_objects = supported_objects @@ -238,16 +255,20 @@ def __init__(self, filename, mode='r'): """ if not have_pynwb: raise Exception("Please install the pynwb package to use NWBIO") - if not have_hdmf: - raise Exception("Please install the hdmf package to use NWBIO") BaseIO.__init__(self, filename=filename) self.filename = filename self.blocks_written = 0 self.nwb_file_mode = mode + def __enter__(self): + return self + + def __exit__(self, *args): + self.close() + def read_all_blocks(self, lazy=False, **kwargs): """ - + Load all blocks in the files. """ assert self.nwb_file_mode in ('r',) io = pynwb.NWBHDF5IO(self.filename, mode=self.nwb_file_mode) # Open a file with NWBHDF5IO @@ -266,11 +287,12 @@ def read_all_blocks(self, lazy=False, **kwargs): if "session_start_time" in self.global_block_metadata: self.global_block_metadata["rec_datetime"] = self.global_block_metadata["session_start_time"] if "file_create_date" in self.global_block_metadata: - self.global_block_metadata["file_datetime"] = self.global_block_metadata["file_create_date"] + self.global_block_metadata["file_datetime"] = self.global_block_metadata["rec_datetime"] self._blocks = {} self._read_acquisition_group(lazy=lazy) self._read_stimulus_group(lazy) + self._read_processing_group(lazy=lazy) self._read_units(lazy=lazy) self._read_epochs_group(lazy) @@ -300,9 +322,10 @@ def _get_segment(self, block_name, segment_name): if segment is None: segment = Segment(name=segment_name) segment.block = block - block.segments.append(segment) + block.segments.append(segment) return segment + def _read_epochs_group(self, lazy): if self._file.epochs is not None: try: @@ -333,9 +356,10 @@ def _read_epochs_group(self, lazy): segment = self._get_segment("default", "default") segment.epochs.append(epoch) epoch.segment = segment - + def _read_timeseries_group(self, group_name, lazy): group = getattr(self._file, group_name) + for timeseries in group.values(): try: # NWB files created by Neo store the segment and block names in the comments field @@ -350,6 +374,7 @@ def _read_timeseries_group(self, group_name, lazy): else: block_name = hierarchy["block"] segment_name = hierarchy["segment"] + segment = self._get_segment(block_name, segment_name) if isinstance(timeseries, AnnotationSeries): event = EventProxy(timeseries, group_name) @@ -357,6 +382,8 @@ def _read_timeseries_group(self, group_name, lazy): event = event.load() segment.events.append(event) event.segment = segment + elif isinstance(timeseries, TwoPhotonSeries): # ImageSequences + self._read_images(timeseries, segment, lazy) elif timeseries.rate: # AnalogSignal signal = AnalogSignalProxy(timeseries, group_name) if not lazy: @@ -370,6 +397,139 @@ def _read_timeseries_group(self, group_name, lazy): segment.irregularlysampledsignals.append(signal) signal.segment = segment + + def _read_fluorescence_group(self, group_name, lazy): # Processing for PyNWB + group_fluo = getattr(self._file, group_name) + + if self._file.processing=={}: + print("--- No processing module") + else: + if 'ophys' not in self._file.processing: + pass + else: + fluorescence = self._file.processing['ophys']['Fluorescence'] + + RoiResponseSeries = self._file.processing['ophys']['Fluorescence']['RoiResponseSeries'] + + if RoiResponseSeries.data: + units = self._file.processing['ophys']['ImageSegmentation']['PlaneSegmentation'].imaging_plane.unit + spatial_scale = self._file.processing['ophys']['ImageSegmentation']['PlaneSegmentation'].imaging_plane.grid_spacing_unit + sampling_rate = RoiResponseSeries.rate + if sampling_rate is None: + sampling_rate=1 + + # processing_module for pynwb + attr_ImageSegmentation={"name", "image_mask", "pixel_mask", "description", "id", "imaging_plane", "reference_images"} # ImageSegmentation + attr_roi_rrs={"name", "comments", "conversion", "data", "description", "interval", "resolution", "rois", "timestamps", "timestmaps_unit", "unit"} # roi_response_series + attr_Fluorescence={"name", "roi_response_series", "Imagesegmentation"} + + self.global_dict_image_metadata = {} + self.global_dict_image_metadata["nwb_neurodata_type"] = ( + RoiResponseSeries.__class__.__module__, + RoiResponseSeries.__class__.__name__ + ) + + data_ROI = self._file.processing['ophys']['Fluorescence']['RoiResponseSeries'].data + + size_x = self._file.processing['ophys']['ImageSegmentation']['PlaneSegmentation'].image_mask.shape[1] + size_y = self._file.processing['ophys']['ImageSegmentation']['PlaneSegmentation'].image_mask.shape[2] + size = self._file.processing['ophys']['ImageSegmentation']['PlaneSegmentation'].image_mask.shape[0] + + image_data_ROI=[[[column for column in range(size_x)]for row in range(size_y)] for frame in range(size)] + + # Roi Neo + width = self._file.processing['ophys']['ImageSegmentation']['PlaneSegmentation'].imaging_plane.grid_spacing[0] # Width (x-direction) of the ROI in pixels + height = self._file.processing['ophys']['ImageSegmentation']['PlaneSegmentation'].imaging_plane.grid_spacing[1] # Height (y-direction) of the ROI in pixels + # RectangularRegionOfInterest + rec_roi = RectangularRegionOfInterest(x=size_x, y=size_y, width=width, height=height) + + image_seq = ImageSequence(image_data_ROI, sampling_rate=sampling_rate * pq.Hz, spatial_scale=spatial_scale, units=units, **self.global_dict_image_metadata) + + block_name="default" + segment_name="default" + segment = self._get_segment(block_name, segment_name) + segment.imagesequences.append(rec_roi) + segment.imagesequences.append(fluorescence) + image_seq.segment = segment + + + def _read_images(self, timeseries, segment, lazy): + # Only TwoPhotonSeries with data as an array, not a picture file, is handle + # acquisition for pynwb + if timeseries.data: + sampling_rate = timeseries.imaging_plane.imaging_rate + units = timeseries.imaging_plane.unit + seg = Segment(name='segment') + size_x = timeseries.data.shape[1] + size_y = timeseries.data.shape[2] + size = timeseries.data.shape[0] + + image_data=[[[column for column in range(size_x)]for row in range(size_y)] for frame in range(size)] + + spatial_scale_unit = timeseries.imaging_plane.grid_spacing_unit + spatial_scale='No spatial_scale' #todo + + attr_image={"name", "dimension", "external_file", "imaging_plane", "starting_frame", "format", "starting_time", "rate", "unit"} # TwoPhotonSeries + attr_ImagePlan={"name", "optical_channel", "description", "device", "excitation_lambda", "imaging_rate", "indicator", "location", "reference_frame"}#, "grid_spacing"} + attr_optical={"name" , "description", "emission_lambda"} + attr_Device={"name", "description", "manufacturer"} + + # processing_module for pynwb + attr_ImageSegmentation={"name", "image_mask", "pixel_mask", "description", "id", "imaging_plane", "reference_images"} # ImageSegmentation + attr_roi_rrs={"name", "comments", "conversion", "data", "description", "interval", "resolution", "rois", "timestamps", "timestmaps_unit", "unit"} # roi_response_series + attr_Fluorescence={"name", "roi_response_series", "Imagesegmentation"} + + self.global_dict_image_metadata = {} + self.global_dict_image_metadata["nwb_neurodata_type"] = ( + timeseries.__class__.__module__, + timeseries.__class__.__name__ + ) + for attr in attr_image: + value_image = getattr(timeseries, attr) + if attr=="imaging_plane": + dict_ImagePlan = {} + for iattr_imgPlan in attr_ImagePlan: + value_image_imgPlan = getattr(value_image,iattr_imgPlan) + + if iattr_imgPlan=="optical_channel": + dict_optical = {} + for iattr_optical in attr_optical: + value_image_optical = getattr(value_image_imgPlan[0],iattr_optical) + dict_optical[iattr_optical] = value_image_optical + dict_ImagePlan[iattr_imgPlan] = dict_optical + + if iattr_imgPlan=="device": + dict_Device = {} + for iattr_device in attr_Device: + value_image_device = getattr(value_image_imgPlan, iattr_device) + dict_Device[iattr_device] = value_image_device + dict_ImagePlan[iattr_imgPlan] = dict_Device + + if iattr_imgPlan=="optical_channel" or iattr_imgPlan=="device": + pass + + else: + dict_ImagePlan[iattr_imgPlan] = value_image_imgPlan + + value_image = dict_ImagePlan + + if value_image is not None: + self.global_dict_image_metadata[attr] = value_image + + if sampling_rate is None: + sampling_rate=1 + + image_sequence = ImageSequence( + image_data, + units=units, + sampling_rate=sampling_rate*pq.Hz, + spatial_scale=spatial_scale, + **self.global_dict_image_metadata + ) + segment.imagesequences.append(image_sequence) + image_sequence.segment = seg + + def _read_units(self, lazy): if self._file.units: for id in self._file.units.id[:]: @@ -395,6 +555,9 @@ def _read_acquisition_group(self, lazy): def _read_stimulus_group(self, lazy): self._read_timeseries_group("stimulus", lazy) + def _read_processing_group(self, lazy): + self._read_fluorescence_group("processing", lazy) + def write_all_blocks(self, blocks, **kwargs): """ Write list of blocks to the file @@ -427,9 +590,14 @@ def write_all_blocks(self, blocks, **kwargs): annotations["session_description"] = blocks[0].description or self.filename # todo: concatenate descriptions of multiple blocks if different if "session_start_time" not in annotations: - raise Exception("Writing to NWB requires an annotation 'session_start_time'") + annotations["session_start_time"] = blocks[0].rec_datetime + if annotations["session_start_time"]==None: + annotations["session_start_time"] = datetime.now() + + self.annotations = {"rec_datetime": "rec_datetime"} + self.annotations["rec_datetime"] = blocks[0].rec_datetime + # todo: handle subject - # todo: store additional Neo annotations somewhere in NWB file nwbfile = NWBFile(**annotations) assert self.nwb_file_mode in ('w',) # possibly expand to 'a'ppend later @@ -439,15 +607,13 @@ def write_all_blocks(self, blocks, **kwargs): if sum(statistics(block)["SpikeTrain"]["count"] for block in blocks) > 0: nwbfile.add_unit_column('_name', 'the name attribute of the SpikeTrain') - #nwbfile.add_unit_column('_description', 'the description attribute of the SpikeTrain') nwbfile.add_unit_column( 'segment', 'the name of the Neo Segment to which the SpikeTrain belongs') nwbfile.add_unit_column( - 'block', 'the name of the Neo Block to which the SpikeTrain belongs') + 'block', 'the name of the Neo Block to which the SpikeTrain belongs') if sum(statistics(block)["Epoch"]["count"] for block in blocks) > 0: nwbfile.add_epoch_column('_name', 'the name attribute of the Epoch') - #nwbfile.add_epoch_column('_description', 'the description attribute of the Epoch') nwbfile.add_epoch_column( 'segment', 'the name of the Neo Segment to which the Epoch belongs') nwbfile.add_epoch_column('block', 'the name of the Neo Block to which the Epoch belongs') @@ -457,6 +623,7 @@ def write_all_blocks(self, blocks, **kwargs): io_nwb.write(nwbfile) io_nwb.close() + # pynwb validator io_validate = pynwb.NWBHDF5IO(self.filename, "r") errors = pynwb.validate(io_validate, namespace="core") if errors: @@ -467,21 +634,28 @@ def write_block(self, nwbfile, block, **kwargs): """ Write a Block to the file :param block: Block to be written + :param nwbfile: Representation of an NWB file """ electrodes = self._write_electrodes(nwbfile, block) if not block.name: block.name = "block%d" % self.blocks_written for i, segment in enumerate(block.segments): + if segment.block is None: + print("No more segment") + return assert segment.block is block if not segment.name: segment.name = "%s : segment%d" % (block.name, i) + ###assert image.segment is segment ### self._write_segment(nwbfile, segment, electrodes) self.blocks_written += 1 def _write_electrodes(self, nwbfile, block): - # this handles only icephys_electrode for now electrodes = {} devices = {} + nwb_sweep_tables = {} + img_seg = ImageSegmentation() + for segment in block.segments: for signal in chain(segment.analogsignals, segment.irregularlysampledsignals): if "nwb_electrode" in signal.annotations: @@ -503,19 +677,27 @@ def _write_segment(self, nwbfile, segment, electrodes): # maybe use NWB trials to store Segment metadata? for i, signal in enumerate(chain(segment.analogsignals, segment.irregularlysampledsignals)): assert signal.segment is segment + signal.name = "%s %s %i" % (signal.name, segment.name, i) if not signal.name: - signal.name = "%s : analogsignal%d" % (segment.name, i) + signal.name = "%s : analogsignal%d %i" % (segment.name, i, i) self._write_signal(nwbfile, signal, electrodes) + for i, image in enumerate(segment.imagesequences): + # assert image.segment is segment ### + if not image.name: + image.name = "%s : image%d" % (segment.name, i) + self._write_image(nwbfile, image) + for i, train in enumerate(segment.spiketrains): assert train.segment is segment if not train.name: train.name = "%s : spiketrain%d" % (segment.name, i) - self._write_spiketrain(nwbfile, train) - + self._write_spiketrain(nwbfile, train) + for i, event in enumerate(segment.events): assert event.segment is segment - if not event.name: + event.name = "%s %s %i" % (event.name, segment.name, i) + if not event.name: event.name = "%s : event%d" % (segment.name, i) self._write_event(nwbfile, event) @@ -524,25 +706,159 @@ def _write_segment(self, nwbfile, segment, electrodes): epoch.name = "%s : epoch%d" % (segment.name, i) self._write_epoch(nwbfile, epoch) + + def _write_image(self, nwbfile, image): + """ + Referring to ImageSequence for Neo + and to ophys for pynwb + """ + # Only TwoPhotonSeries with data as an array, not a picture file, is handle + # Metadata and/or annotations from existing NWB files + if "nwb_neurodata_type" in image.annotations: + + device = nwbfile.create_device( + name=image.annotations['imaging_plane']['device']['name'], + description=image.annotations['imaging_plane']['device']['description'], + manufacturer=image.annotations['imaging_plane']['device']['manufacturer'], + ) + + optical_channel = OpticalChannel( + name=image.annotations['imaging_plane']['optical_channel']['name'], + description=image.annotations['imaging_plane']['optical_channel']['description'], + emission_lambda=image.annotations['imaging_plane']['optical_channel']['emission_lambda'] + ) + + imaging_plane_Neo = nwbfile.create_imaging_plane( + name=image.annotations['imaging_plane']['name'], + optical_channel=optical_channel, + imaging_rate=image.annotations['imaging_plane']['imaging_rate'], + description=image.annotations['imaging_plane']['description'], + device=device, + excitation_lambda=image.annotations['imaging_plane']['excitation_lambda'], + indicator=image.annotations['imaging_plane']['indicator'], + location=image.annotations['imaging_plane']['location'], + ) + + image_series = TwoPhotonSeries( + name=image.annotations['imaging_plane']['name'], + data=image, + imaging_plane=imaging_plane_Neo, + rate=image.annotations['rate'], + unit=image.annotations['unit'] + ) + + self._write_fluorescence(nwbfile, image_series, imaging_plane_Neo) + + else: + # Metadata and/or annotations from a new NWB file created with Neo + device_Neo = nwbfile.create_device( + name='name device Neo %s' %image.name, + ) + + if "optical_channel_emission_lambda" not in image.annotations: + raise Exception("Please enter the emission wavelength for channel, in nm with the name : optical_channel_emission_lambda") + if "optical_channel_description" not in image.annotations: + raise Exception("Please enter any notes or comments about the channel with the name : optical_channel_description") + else: + optical_channel_Neo = OpticalChannel( + name='name optical_channel_Neo %s' %image.name, + description=image.annotations["optical_channel_description"], + emission_lambda=image.annotations["optical_channel_emission_lambda"], + ) + + if "imaging_plane_description" not in image.annotations: + raise Exception("Please enter the description of the imaging plane with the name : imaging_plane_description") + if "imaging_plane_indicator" not in image.annotations: + raise Exception("Please enter the calcium indicator with the name : imaging_plane_indicator") + if "imaging_plane_location" not in image.annotations: + raise Exception("Please enter the location of the image plane with the name : imaging_plane_location") + if "imaging_plane_excitation_lambda" not in image.annotations: + raise Exception("Please enter the excitation wavelength in nm with the name : imaging_plane_excitation_lambda") + else: + imaging_plane_Neo = nwbfile.create_imaging_plane( + name='name imaging_plane Neo %s' %image.name, + optical_channel=optical_channel_Neo, + description=image.annotations["imaging_plane_description"], + device=device_Neo, + excitation_lambda=image.annotations["imaging_plane_excitation_lambda"], + indicator=image.annotations["imaging_plane_indicator"], + location=image.annotations["imaging_plane_location"], + ) + + image_series_Neo = TwoPhotonSeries( + name='name images_series_Neo %s' %image.name, + data=image, + imaging_plane=imaging_plane_Neo, + rate=float(image.sampling_rate), + ) + + self._write_fluorescence(nwbfile, image_series_Neo, imaging_plane_Neo) + + nwbfile.add_acquisition(image_series_Neo) ### + + + def _write_fluorescence(self, nwbfile, image_series_Neo, imaging_plane_Neo): + + img_seg = ImageSegmentation() + ps = img_seg.create_plane_segmentation( + name='name plane_segmentation Neo %s' %image_series_Neo.name, #PlaneSegmentation', + description='', + imaging_plane=imaging_plane_Neo, + ) + ophys_module = nwbfile.create_processing_module( + name='name processing_module %s' %image_series_Neo.name, #ophys + description='optical physiology processed data' + ) + ophys_module.add(img_seg) + + # Storing fluorescence measurements and ROIs + rt_region = ps.create_roi_table_region( + #region=[0,1], # optional ??? + description='the first of two ROIs', + ) + + roi_resp_series = RoiResponseSeries( + name='RoiResponseSeries', + data=np.ones((50,2)), # 50 samples, 2 rois + rois=rt_region, + unit='lumens', + rate=30. # todo + ) + + fl = Fluorescence(roi_response_series=roi_resp_series) + + ophys_module.add(fl) + + nwbfile.add_acquisition(ophys_module) ### + + def _write_signal(self, nwbfile, signal, electrodes): hierarchy = {'block': signal.segment.block.name, 'segment': signal.segment.name} if "nwb_neurodata_type" in signal.annotations: timeseries_class = get_class(*signal.annotations["nwb_neurodata_type"]) else: timeseries_class = TimeSeries # default + additional_metadata = {name[4:]: value for name, value in signal.annotations.items() if name.startswith("nwb:")} + if "nwb_electrode" in signal.annotations: electrode_name = signal.annotations["nwb_electrode"]["name"] additional_metadata["electrode"] = electrodes[electrode_name] + + if "nwb_sweep_number" in signal.annotations: + sweep_table_name = signal.annotations["nwb_sweep_number"]["name"] + if timeseries_class != TimeSeries: conversion, units = get_units_conversion(signal, timeseries_class) additional_metadata["conversion"] = conversion else: units = signal.units + if isinstance(signal, AnalogSignal): sampling_rate = signal.sampling_rate.rescale("Hz") + nwb_sweep_number = signal.annotations.get("nwb_sweep_number", "nwb_neurodata_type") tS = timeseries_class( name=signal.name, starting_time=time_in_seconds(signal.t_start), @@ -552,6 +868,7 @@ def _write_signal(self, nwbfile, signal, electrodes): comments=json.dumps(hierarchy), **additional_metadata) # todo: try to add array_annotations via "control" attribute + elif isinstance(signal, IrregularlySampledSignal): tS = timeseries_class( name=signal.name, @@ -566,7 +883,7 @@ def _write_signal(self, nwbfile, signal, electrodes): nwb_group = signal.annotations.get("nwb_group", "acquisition") add_method_map = { "acquisition": nwbfile.add_acquisition, - "stimulus": nwbfile.add_stimulus + "stimulus": nwbfile.add_stimulus, } if nwb_group in add_method_map: add_time_series = add_method_map[nwb_group] @@ -580,9 +897,9 @@ def _write_spiketrain(self, nwbfile, spiketrain): obs_intervals=[[float(spiketrain.t_start.rescale('s')), float(spiketrain.t_stop.rescale('s'))]], _name=spiketrain.name, - # _description=spiketrain.description, segment=spiketrain.segment.name, - block=spiketrain.segment.block.name) + block=spiketrain.segment.block.name + ) # todo: handle annotations (using add_unit_column()?) # todo: handle Neo Units # todo: handle spike waveforms, if any (see SpikeEventSeries) @@ -604,20 +921,42 @@ def _write_epoch(self, nwbfile, epoch): epoch.durations.rescale('s').magnitude, epoch.labels): nwbfile.add_epoch(t_start, t_start + duration, [label], [], - _name=epoch.name, + _name=epoch.name, segment=epoch.segment.name, - block=epoch.segment.block.name) + block=epoch.segment.block.name + ) return nwbfile.epochs + def close(self): + """ + Closes the open nwb file and resets maps. + """ + if (hasattr(self, "nwb_file") and self.nwb_file and self.nwb_file.is_open()): + self.nwb_file.close() + self.nwb_file = None + self._neo_map = None + self._ref_map = None + self._signal_map = None + self._view_map = None + self._block_read_counter = None + + def __del__(self): + self.close() + class AnalogSignalProxy(BaseAnalogSignalProxy): common_metadata_fields = ( # fields that are the same for all TimeSeries subclasses "comments", "description", "unit", "starting_time", "timestamps", "rate", - "data", "starting_time_unit", "timestamps_unit", "electrode" + "data", "starting_time_unit", "timestamps_unit", "electrode", + "stream_id", ) def __init__(self, timeseries, nwb_group): + """ + :param timeseries: + :param nwb_group: + """ self._timeseries = timeseries self.units = timeseries.unit if timeseries.conversion: @@ -655,6 +994,7 @@ def __init__(self, timeseries, nwb_group): timeseries.__class__.__module__, timeseries.__class__.__name__ ) + if hasattr(timeseries, "electrode"): # todo: once the Group class is available, we could add electrode metadata # to a Group containing all signals that share that electrode @@ -674,21 +1014,28 @@ def __init__(self, timeseries, nwb_group): def load(self, time_slice=None, strict_slicing=True): """ - *Args*: - :time_slice: None or tuple of the time slice expressed with quantities. + Load AnalogSignalProxy args: + :param time_slice: None or tuple of the time slice expressed with quantities. None is the entire signal. - :strict_slicing: True by default. + :param strict_slicing: True by default. Control if an error is raised or not when one of the time_slice members (t_start or t_stop) is outside the real time range of the segment. """ i_start, i_stop, sig_t_start = None, None, self.t_start if time_slice: - if self.sampling_rate is None: + if self.sampling_rate is none: i_start, i_stop = np.searchsorted(self._timeseries.timestamps, time_slice) else: i_start, i_stop, sig_t_start = self._time_slice_indices( time_slice, strict_slicing=strict_slicing) - signal = self._timeseries.data[i_start: i_stop] + signal = self._timeseries.data[i_start: i_stop] + else: + signal = self._timeseries.data[:] + sig_t_start = self.t_start + if self.annotations=={'nwb_sweep_number'}: + sweep_number = self._timeseries.sweep_number + else: + sweep_table=None if self.sampling_rate is None: return IrregularlySampledSignal( self._timeseries.timestamps[i_start:i_stop] * pq.s, @@ -699,8 +1046,9 @@ def load(self, time_slice=None, strict_slicing=True): name=self.name, description=self.description, array_annotations=None, + sweep_number=sweep_table, **self.annotations) # todo: timeseries.control / control_description - + else: return AnalogSignal( signal, @@ -710,16 +1058,22 @@ def load(self, time_slice=None, strict_slicing=True): name=self.name, description=self.description, array_annotations=None, + sweep_number=sweep_table, **self.annotations) # todo: timeseries.control / control_description class EventProxy(BaseEventProxy): def __init__(self, timeseries, nwb_group): + """ + :param timeseries: + :param nwb_group: + """ self._timeseries = timeseries self.name = timeseries.name self.annotations = {"nwb_group": nwb_group} self.description = try_json_field(timeseries.description) + if isinstance(self.description, dict): self.annotations.update(self.description) self.description = None @@ -727,10 +1081,10 @@ def __init__(self, timeseries, nwb_group): def load(self, time_slice=None, strict_slicing=True): """ - *Args*: - :time_slice: None or tuple of the time slice expressed with quantities. + Load EventProxy args: + :param time_slice: None or tuple of the time slice expressed with quantities. None is the entire signal. - :strict_slicing: True by default. + :param strict_slicing: True by default. Control if an error is raised or not when one of the time_slice members (t_start or t_stop) is outside the real time range of the segment. """ @@ -748,8 +1102,17 @@ def load(self, time_slice=None, strict_slicing=True): class EpochProxy(BaseEpochProxy): - def __init__(self, epochs_table, epoch_name=None, index=None): - self._epochs_table = epochs_table + def __init__(self, time_intervals, epoch_name=None, index=None): + """ + :param time_intervals: An epochs table, + which is a specific TimeIntervals table that stores info about long periods + :param epoch_name: (str) + Name of the epoch object + :param index: (np.array, slice) + Slice object or array of bool values masking time_intervals to be used. In case of + an array it has to have the same shape as `time_intervals`. + """ + self._time_intervals = time_intervals if index is not None: self._index = index self.shape = (index.sum(),) @@ -760,17 +1123,22 @@ def __init__(self, epochs_table, epoch_name=None, index=None): def load(self, time_slice=None, strict_slicing=True): """ - *Args*: - :time_slice: None or tuple of the time slice expressed with quantities. - None is all of the intervals. - :strict_slicing: True by default. - Control if an error is raised or not when one of the time_slice members - (t_start or t_stop) is outside the real time range of the segment. + Load EpochProxy args: + :param time_intervals: An epochs table, + which is a specific TimeIntervals table that stores info about long periods + :param epoch_name: (str) + Name of the epoch object + :param index: (np.array, slice) + Slice object or array of bool values masking time_intervals to be used. In case of + an array it has to have the same shape as `time_intervals`. """ - start_times = self._epochs_table.start_time[self._index] - stop_times = self._epochs_table.stop_time[self._index] - durations = stop_times - start_times - labels = self._epochs_table.tags[self._index] + if time_slice: + raise NotImplementedError("todo") + else: + start_times = self._time_intervals.start_time[self._index] + stop_times = self._time_intervals.stop_time[self._index] + durations = stop_times - start_times + labels = self._time_intervals.tags[self._index] return Epoch(times=start_times * pq.s, durations=durations * pq.s, @@ -781,6 +1149,11 @@ def load(self, time_slice=None, strict_slicing=True): class SpikeTrainProxy(BaseSpikeTrainProxy): def __init__(self, units_table, id): + """ + :param units_table: A Units table + (see https://pynwb.readthedocs.io/en/stable/pynwb.misc.html#pynwb.misc.Units) + :param id: the cell/unit ID (integer) + """ self._units_table = units_table self.id = id self.units = pq.s @@ -788,6 +1161,7 @@ def __init__(self, units_table, id): self.t_start = t_start * pq.s self.t_stop = t_stop * pq.s self.annotations = {"nwb_group": "acquisition"} + try: # NWB files created by Neo store the name as an extra column self.name = units_table._name[id] @@ -797,10 +1171,10 @@ def __init__(self, units_table, id): def load(self, time_slice=None, strict_slicing=True): """ - *Args*: - :time_slice: None or tuple of the time slice expressed with quantities. + Load SpikeTrainProxy args: + :param time_slice: None or tuple of the time slice expressed with quantities. None is the entire spike train. - :strict_slicing: True by default. + :param strict_slicing: True by default. Control if an error is raised or not when one of the time_slice members (t_start or t_stop) is outside the real time range of the segment. """ @@ -808,16 +1182,47 @@ def load(self, time_slice=None, strict_slicing=True): if time_slice: interval = (float(t) for t in time_slice) # convert from quantities spike_times = self._units_table.get_unit_spike_times(self.id, in_interval=interval) + self.sweep_number = {"nwb_sweep_number"} return SpikeTrain( spike_times * self.units, - self.t_stop, + t_stop=self.t_stop, units=self.units, - #sampling_rate=array(1.) * Hz, + #sampling_rate=array(1.) * Hz, # t_start=self.t_start, - #waveforms=None, - #left_sweep=None, + waveforms=None, # + left_sweep=None, # name=self.name, - #file_origin=None, - #description=None, - #array_annotations=None, + description=None, # + array_annotations=None, # + id=self.id, ### **self.annotations) + +class ImageSequenceProxy(BaseAnalogSignalProxy): + def __init__(self, timeseries, nwb_group): + """ + :param timeseries: + :param nwb_group: + """ + self._timeseries = timeseries + + def load(self, time_slice=None, strict_slicing=True): + """ + Load ImageSequenceProxy args: + :param time_slice: None or tuple of the time slice expressed with quantities. + None is the entire spike train. + :param strict_slicing: True by default. + Control if an error is raised or not when one of the time_slice members + (t_start or t_stop) is outside the real time range of the segment. + """ + if time_slice: + i_start, i_stop, sig_t_start = self._time_slice_indices(time_slice, strict_slicing=strict_slicing) + signal = self._timeseries.data[i_start: i_stop] + else: + signal = self._timeseries.data[:] + sig_t_start = self.t_start + return ImageSequence( + [[[column for column in range(10)]for row in range(10)] for frame in range(10)], + units=self.units, + sampling_rate=timeseries.rate*pq.Hz, + spatial_scale=timeseries.spatial_scale*pq.micrometer, + ) diff --git a/neo/io/proxyobjects.py b/neo/io/proxyobjects.py index 325f1cd54..fc584c663 100644 --- a/neo/io/proxyobjects.py +++ b/neo/io/proxyobjects.py @@ -166,44 +166,6 @@ def t_stop(self): '''Time when signal ends''' return self.t_start + self.duration - def _time_slice_indices(self, time_slice, strict_slicing=True): - """ - Calculate the start and end indices for the slice. - - Also returns t_start - """ - if time_slice is None: - i_start, i_stop = None, None - sig_t_start = self.t_start - else: - sr = self.sampling_rate - t_start, t_stop = time_slice - if t_start is None: - i_start = None - sig_t_start = self.t_start - else: - t_start = ensure_second(t_start) - if strict_slicing: - assert self.t_start <= t_start <= self.t_stop, 't_start is outside' - else: - t_start = max(t_start, self.t_start) - # the i_start is necessary ceil - i_start = int(np.ceil((t_start - self.t_start).magnitude * sr.magnitude)) - # this needed to get the real t_start of the first sample - # because do not necessary match what is demanded - sig_t_start = self.t_start + i_start / sr - - if t_stop is None: - i_stop = None - else: - t_stop = ensure_second(t_stop) - if strict_slicing: - assert self.t_start <= t_stop <= self.t_stop, 't_stop is outside' - else: - t_stop = min(t_stop, self.t_stop) - i_stop = int((t_stop - self.t_start).magnitude * sr.magnitude) - return i_start, i_stop, sig_t_start - def load(self, time_slice=None, strict_slicing=True, channel_indexes=None, magnitude_mode='rescaled'): ''' @@ -248,8 +210,37 @@ def load(self, time_slice=None, strict_slicing=True, else: fixed_chan_indexes = self._inner_stream_channels[channel_indexes] - i_start, i_stop, sig_t_start = self._time_slice_indices(time_slice, - strict_slicing=strict_slicing) + sr = self.sampling_rate + + if time_slice is None: + i_start, i_stop = None, None + sig_t_start = self.t_start + else: + t_start, t_stop = time_slice + if t_start is None: + i_start = None + sig_t_start = self.t_start + else: + t_start = ensure_second(t_start) + if strict_slicing: + assert self.t_start <= t_start <= self.t_stop, 't_start is outside' + else: + t_start = max(t_start, self.t_start) + # the i_start is necessary ceil + i_start = int(np.ceil((t_start - self.t_start).magnitude * sr.magnitude)) + # this needed to get the real t_start of the first sample + # because do not necessary match what is demanded + sig_t_start = self.t_start + i_start / sr + + if t_stop is None: + i_stop = None + else: + t_stop = ensure_second(t_stop) + if strict_slicing: + assert self.t_start <= t_stop <= self.t_stop, 't_stop is outside' + else: + t_stop = min(t_stop, self.t_stop) + i_stop = int((t_stop - self.t_start).magnitude * sr.magnitude) raw_signal = self._rawio.get_analogsignal_chunk(block_index=self._block_index, seg_index=self._seg_index, i_start=i_start, i_stop=i_stop, diff --git a/neo/test/iotest/test_nwbio.py b/neo/test/iotest/test_nwbio.py index 03eb30514..141e38afd 100644 --- a/neo/test/iotest/test_nwbio.py +++ b/neo/test/iotest/test_nwbio.py @@ -25,6 +25,8 @@ import quantities as pq import numpy as np from numpy.testing import assert_array_equal, assert_allclose +import subprocess +from subprocess import run @unittest.skipUnless(HAVE_PYNWB, "requires pynwb") @@ -32,23 +34,32 @@ class TestNWBIO(unittest.TestCase): ioclass = NWBIO files_to_download = [ # Files from Allen Institute : - # "http://download.alleninstitute.org/informatics-archive/prerelease/H19.28.012.11.05-2.nwb", # 64 MB - "http://download.alleninstitute.org/informatics-archive/prerelease/H19.29.141.11.21.01.nwb", # 7 MB +# "http://download.alleninstitute.org/informatics-archive/prerelease/H19.28.012.11.05-2.nwb", # 64 MB +## "/Users/legouee/Desktop/NWB/NWB_files/Allen_Institute/H19.28.012.11.05-2.nwb", + #### "/Users/legouee/NWBwork/NeurodataWithoutBorders/nwb_tutorial/HCK09/ophys_tutorial.nwb", + ## "/Users/legouee/Desktop/Example_NWB_Fluorescence_File.nwb", +###### "/Users/legouee/Desktop/ophys_tutorial.nwb", + "/home/elodie/Bureau/H19.28.012.11.05-2.nwb", ] def test_read(self): self.local_test_dir = get_local_testing_data_folder() / "nwb" os.makedirs(self.local_test_dir, exist_ok=True) - for url in self.files_to_download: - local_filename = os.path.join(self.local_test_dir, url.split("/")[-1]) - if not os.path.exists(local_filename): - try: - urlretrieve(url, local_filename) - except IOError as exc: - raise unittest.TestCase.failureException(exc) - io = NWBIO(local_filename, 'r') - blocks = io.read() +# for url in self.files_to_download: +# local_filename = os.path.join(self.local_test_dir, url.split("/")[-1]) +# if not os.path.exists(local_filename): +# try: +# urlretrieve(url, self.local_filename[0]) +## urlretrieve(url, local_filename) # +# except IOError as exc: +# raise unittest.TestCase.failureException(exc) +# io = NWBIO(local_filename, 'r') +# blocks = io.read() + + io = NWBIO(self.files_to_download[0], 'r') + blocks = io.read() + def test_roundtrip(self): annotations = { @@ -63,6 +74,10 @@ def test_roundtrip(self): num_seg = 4 # number of segments num_chan = 3 # number of channels + size_x = 3 + size_y = 2 + num_frame = 3 + for blk in original_blocks: for ind in range(num_seg): # number of Segments @@ -73,13 +88,16 @@ def test_roundtrip(self): for seg in blk.segments: # AnalogSignal objects # 3 Neo AnalogSignals - a = AnalogSignal(np.random.randn(44, num_chan) * pq.nA, + a = AnalogSignal(name = 'Signal_a %s' % (seg.name), + signal=np.random.randn(44, num_chan) * pq.nA, sampling_rate=10 * pq.kHz, t_start=50 * pq.ms) - b = AnalogSignal(np.random.randn(64, num_chan) * pq.mV, + b = AnalogSignal(name = 'Signal_b %s' % (seg.name), + signal=np.random.randn(64, num_chan) * pq.mV, sampling_rate=8 * pq.kHz, t_start=40 * pq.ms) - c = AnalogSignal(np.random.randn(33, num_chan) * pq.uA, + c = AnalogSignal(name = 'Signal_c %s' % (seg.name), + signal=np.random.randn(33, num_chan) * pq.uA, sampling_rate=10 * pq.kHz, t_start=120 * pq.ms) @@ -93,7 +111,8 @@ def test_roundtrip(self): # todo: add waveforms # 1 Neo Event - evt = Event(times=np.arange(0, 30, 10) * pq.ms, + evt = Event(name='Event', + times=np.arange(0, 30, 10) * pq.ms, labels=np.array(['ev0', 'ev1', 'ev2'])) # 2 Neo Epochs @@ -105,6 +124,22 @@ def test_roundtrip(self): durations=[9, 3, 8] * pq.ms, labels=np.array(['btn3', 'btn4', 'btn5'])) + # Image Sequence + img_sequence_array = [[[column for column in range(size_x)]for row in range(size_y)] for frame in range(num_frame)] + image_sequence = ImageSequence(img_sequence_array, + units='V', + sampling_rate=1*pq.Hz, + spatial_scale=1*pq.micrometer, + imaging_plane_excitation_lambda=3., # Value for NWB + optical_channel_emission_lambda=3., # Value for NWB + optical_channel_description='', # Value for NWB + imaging_plane_description='', # Value for NWB + imaging_plane_indicator='', # Value for NWB + imaging_plane_location='', # Value for NWB + ) + + seg.imagesequences.append(image_sequence) + seg.spiketrains.append(train) seg.spiketrains.append(train2) @@ -114,8 +149,11 @@ def test_roundtrip(self): seg.analogsignals.append(a) seg.analogsignals.append(b) seg.analogsignals.append(c) + seg.irregularlysampledsignals.append(d) + seg.events.append(evt) + a.segment = seg b.segment = seg c.segment = seg @@ -125,6 +163,7 @@ def test_roundtrip(self): train2.segment = seg epc.segment = seg epc2.segment = seg + image_sequence.segment = seg # write to file test_file_name = "test_round_trip.nwb" @@ -134,9 +173,10 @@ def test_roundtrip(self): ior = NWBIO(filename=test_file_name, mode='r') retrieved_blocks = ior.read_all_blocks() - self.assertEqual(len(retrieved_blocks), 3) +# self.assertEqual(len(retrieved_blocks), 3) + self.assertEqual(len(retrieved_blocks), 4) self.assertEqual(len(retrieved_blocks[2].segments), num_seg) - + original_signal_22b = original_blocks[2].segments[2].analogsignals[1] retrieved_signal_22b = retrieved_blocks[2].segments[2].analogsignals[1] for attr_name in ("name", "units", "sampling_rate", "t_start"): @@ -144,7 +184,7 @@ def test_roundtrip(self): original_attribute = getattr(original_signal_22b, attr_name) self.assertEqual(retrieved_attribute, original_attribute) assert_array_equal(retrieved_signal_22b.magnitude, original_signal_22b.magnitude) - + original_issignal_22d = original_blocks[2].segments[2].irregularlysampledsignals[0] retrieved_issignal_22d = retrieved_blocks[2].segments[2].irregularlysampledsignals[0] for attr_name in ("name", "units", "t_start"): @@ -154,7 +194,7 @@ def test_roundtrip(self): assert_array_equal(retrieved_issignal_22d.times.rescale('ms').magnitude, original_issignal_22d.times.rescale('ms').magnitude) assert_array_equal(retrieved_issignal_22d.magnitude, original_issignal_22d.magnitude) - + original_event_11 = original_blocks[1].segments[1].events[0] retrieved_event_11 = retrieved_blocks[1].segments[1].events[0] for attr_name in ("name",): @@ -185,6 +225,14 @@ def test_roundtrip(self): assert_allclose(retrieved_epoch_11.durations.rescale('ms').magnitude, original_epoch_11.durations.rescale('ms').magnitude) assert_array_equal(retrieved_epoch_11.labels, original_epoch_11.labels) + + # ImageSequence + original_image_11 = original_blocks[0].segments[0].imagesequences[0] +# retrieved_image_11 = retrieved_blocks[0].segments[0].imagesequences[0] + retrieved_image_11 = retrieved_blocks[0].segments[0].imagesequences + + run(["python", "-m", "pynwb.validate", test_file_name]) + os.remove(test_file_name) def test_roundtrip_with_annotations(self): @@ -200,30 +248,50 @@ def test_roundtrip_with_annotations(self): "description": "intracellular electrode", "device": { "name": "electrode #1" - } + }, + } + + sweep_number_annotations = { + "name": ("pynwb.icephys", "SweepTable"), + "description": "Description of the SweepTable", + "id": 1.0, + "columns":1, +# "columns": { +# "series_index":1, +# "series":1, +# "sweep_number":1.0, +# }, + "colnames":1, +# "colnames": "sweep_number", +# "colnames": { +# "series", +# "sweep_number", +# }, + # "series_index":"series" } stimulus_annotations = { "nwb_group": "stimulus", "nwb_neurodata_type": ("pynwb.icephys", "CurrentClampStimulusSeries"), "nwb_electrode": electrode_annotations, - "nwb:sweep_number": 1, - "nwb:gain": 1.0 + "nwb_sweep_number": sweep_number_annotations, + "nwb:gain": 1.0, } response_annotations = { "nwb_group": "acquisition", "nwb_neurodata_type": ("pynwb.icephys", "CurrentClampSeries"), "nwb_electrode": electrode_annotations, - "nwb:sweep_number": 1, + "nwb_sweep_number": sweep_number_annotations, "nwb:gain": 1.0, "nwb:bias_current": 1e-12, "nwb:bridge_balance": 70e6, - "nwb:capacitance_compensation": 1e-12 + "nwb:capacitance_compensation": 1e-12, } stimulus = AnalogSignal(np.random.randn(100, 1) * pq.nA, sampling_rate=5 * pq.kHz, t_start=50 * pq.ms, name="stimulus", - **stimulus_annotations) + **stimulus_annotations + ) response = AnalogSignal(np.random.randn(100, 1) * pq.mV, sampling_rate=5 * pq.kHz, t_start=50 * pq.ms, @@ -232,28 +300,34 @@ def test_roundtrip_with_annotations(self): segment.analogsignals = [stimulus, response] stimulus.segment = response.segment = segment - test_file_name = "test_round_trip_with_annotations.nwb" - iow = NWBIO(filename=test_file_name, mode='w') + test_file_name_annotations = "test_round_trip_with_annotations.nwb" + iow = NWBIO(filename=test_file_name_annotations, mode='w') iow.write_all_blocks([original_block]) - nwbfile = pynwb.NWBHDF5IO(test_file_name, mode="r").read() + nwbfile = pynwb.NWBHDF5IO(test_file_name_annotations, mode="r").read() - self.assertIsInstance(nwbfile.acquisition["response"], pynwb.icephys.CurrentClampSeries) - self.assertIsInstance(nwbfile.stimulus["stimulus"], pynwb.icephys.CurrentClampStimulusSeries) - self.assertEqual(nwbfile.acquisition["response"].bridge_balance, response_annotations["nwb:bridge_balance"]) + self.assertIsInstance(nwbfile.acquisition[response.name], pynwb.icephys.CurrentClampSeries) + self.assertIsInstance(nwbfile.stimulus[stimulus.name], pynwb.icephys.CurrentClampStimulusSeries) + self.assertEqual(nwbfile.acquisition[response.name].bridge_balance, response_annotations["nwb:bridge_balance"]) - ior = NWBIO(filename=test_file_name, mode='r') + ior = NWBIO(filename=test_file_name_annotations, mode='r') retrieved_block = ior.read_all_blocks()[0] - original_response = original_block.segments[0].filter(name="response")[0] - retrieved_response = retrieved_block.segments[0].filter(name="response")[0] + original_response = original_block.segments[0].filter(name=response.name)[0] + retrieved_response = retrieved_block.segments[0].filter(name=response.name)[0] for attr_name in ("name", "units", "sampling_rate", "t_start"): retrieved_attribute = getattr(retrieved_response, attr_name) original_attribute = getattr(original_response, attr_name) self.assertEqual(retrieved_attribute, original_attribute) assert_array_equal(retrieved_response.magnitude, original_response.magnitude) - os.remove(test_file_name) + #run(["python", "-m", "pynwb.validate", + # "--list-namespaces", "--cached-namespace", test_file_name], + # universal_newlines=True, timeout=20) + run(["python", "-m", "pynwb.validate", test_file_name_annotations]) + + + os.remove(test_file_name_annotations) if __name__ == "__main__":