Store ElectricalSeries in uV for spyglass compatability

src/jdb_to_nwb/convert_raw_ephys.py

@@ -21,15 +21,15 @@
 )
 from pynwb import NWBFile
 from pynwb.ecephys import ElectricalSeries
+from hdmf.data_utils import DataChunkIterator
 from spikeinterface.extractors import OpenEphysBinaryRecordingExtractor
 
 from .utils import get_logger_directory
 from .timestamps_alignment import align_via_interpolation
 from .plotting.plot_ephys import plot_channel_map, plot_channel_impedances, plot_neuropixels
 from ndx_franklab_novela import AssociatedFiles, Probe, NwbElectrodeGroup, Shank, ShanksElectrode
 
-MICROVOLTS_PER_VOLT = 1e6
-VOLTS_PER_MICROVOLT = 1 / MICROVOLTS_PER_VOLT
+VOLTS_PER_MICROVOLT = 1e-6
 
 MIN_IMPEDANCE_OHMS = 1e5
 MAX_IMPEDANCE_OHMS = 3e6
@@ -248,7 +248,7 @@ def get_raw_ephys_data(
         traces_as_iterator (SpikeInterfaceRecordingDataChunkIterator):
             To be used as the data argument in pynwb.ecephys.ElectricalSeries.
         channel_conversion_factor (float):
-            The conversion factor from the raw data to volts.
+            The conversion factor from the raw data to uV.
         original_timestamps (np.ndarray):
             Array that could be used as the timestamps argument in pynwb.ecephys.ElectricalSeries
             or may need to be time aligned with the other data streams in the NWB file.
@@ -298,8 +298,9 @@ def get_raw_ephys_data(
             "The channel conversion factors are not the same for all channels. "
             "This is unexpected and may indicate a problem with the conversion factors."
         )
-    channel_conversion_factor_v = channel_conversion_factors_uv[0] * VOLTS_PER_MICROVOLT
-    logger.debug(f"Channel conversion factor in V: {channel_conversion_factor_v}")
+    # Just grab the first one, because it should be the same for all channels
+    channel_conversion_factor_uv = channel_conversion_factors_uv[0]
+    logger.debug(f"Channel conversion factor in uV: {channel_conversion_factor_uv}")
 
     # NOTE channel offsets should be 0 for all channels in openephys data
     channel_conversion_offsets = recording_sliced.get_channel_offsets()
@@ -317,7 +318,7 @@ def get_raw_ephys_data(
 
     return (
         traces_as_iterator,
-        channel_conversion_factor_v,
+        channel_conversion_factor_uv,
         original_timestamps,
     )
 
@@ -1129,7 +1130,7 @@ def add_raw_ephys(
     # Get raw ephys data
     (
         traces_as_iterator,
-        channel_conversion_factor_v,
+        channel_conversion_factor_uv,
         original_timestamps,
     ) = get_raw_ephys_data(folder_path=openephys_folder_path, logger=logger, exclude_channels=exclude_channel_names)
     num_samples, num_channels = traces_as_iterator.maxshape
@@ -1165,13 +1166,26 @@ def add_raw_ephys(
         description="Electrodes used in raw ElectricalSeries recording",
     )
 
+    # Convert to uV without loading the whole thing at once
+    def traces_in_microvolts_iterator(traces_as_iterator, conversion_factor_uv):
+        for chunk in traces_as_iterator:
+            yield (chunk * conversion_factor_uv).astype("int16")
+
+    # Wrap iterator in DataChunkIterator for H5DataIO
+    data_iterator = DataChunkIterator(
+        traces_in_microvolts_iterator(traces_as_iterator, channel_conversion_factor_uv),
+        buffer_size=1,  # number of chunks to keep in memory
+        maxshape=(num_samples, num_channels),
+        dtype=np.dtype("int16"),
+    )
+
     # A chunk of shape (81920, 64) and dtype int16 (2 bytes) is ~10 MB, which is the recommended chunk size
     # by the NWB team.
     # We could also add compression here. zstd/blosc-zstd are recommended by the NWB team, but
     # they require the hdf5plugin library to be installed. gzip is available by default.
     # Use gzip for now, but consider zstd/blosc-zstd in the future.
     data_data_io = H5DataIO(
-        traces_as_iterator,
+        data_iterator,
         chunks=(min(num_samples, 81920), min(num_channels, 64)),
         compression="gzip",
     )
@@ -1200,11 +1214,11 @@ def add_raw_ephys(
     # For now, we do not chunk or compress the timestamps, which are relatively small
     eseries = ElectricalSeries(
         name="ElectricalSeries",
-        description="Raw ephys data from OpenEphys recording (multiply by conversion factor to get data in volts).",
+        description="Raw ephys data from OpenEphys recording, in uV (multiply by conversion factor to get data in V).",
         data=data_data_io,
         timestamps=ephys_timestamps,
         electrodes=electrode_table_region,
-        conversion=channel_conversion_factor_v,
+        conversion=VOLTS_PER_MICROVOLT, # conversion from uV to V
     )
 
     # Add the ElectricalSeries to the NWBFile

tests/test_convert_raw_ephys.py

@@ -209,7 +209,7 @@ def test_get_raw_ephys_data(dummy_logger):
     folder_path = "tests/test_data/raw_ephys/2022-07-25_15-30-00"
     traces_as_iterator, channel_conversion_factor, original_timestamps = get_raw_ephys_data(folder_path, dummy_logger)
     assert traces_as_iterator.maxshape == (3_000, 256)
-    np.testing.assert_allclose(channel_conversion_factor, [0.19499999284744263 * 1e-6] * 256)
+    np.testing.assert_allclose(channel_conversion_factor, [0.19499999284744263] * 256)
     assert len(original_timestamps) == 3_000
 
 
@@ -255,7 +255,7 @@ def test_add_raw_ephys(dummy_logger):
     metadata["ephys"]["targeted_z"] = -2.0  # DV in mm
     metadata["ephys"]["probe"] = ["256-ch Silicon Probe, 3mm length, 66um pitch"]
 
-    # Mapping from probe electrode to channel number
+    # Mapping from probe electrode to channel number for 256-ch silicon probe
     intan_channel_numbers = np.array([
         191,190,189,188,187,186,185,184,183,182,181,180,179,178,177,176,175,174,173,172,171,170,169,128,129,130,
         131,132,133,134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151,152,153,154,155,156,
@@ -278,13 +278,13 @@ def test_add_raw_ephys(dummy_logger):
     assert "ElectricalSeries" in nwbfile.acquisition
     es = nwbfile.acquisition["ElectricalSeries"]
     assert es.description == (
-        "Raw ephys data from OpenEphys recording (multiply by conversion factor to get data in volts)."
+        "Raw ephys data from OpenEphys recording, in uV (multiply by conversion factor to get data in V)."
     )
     assert es.data.maxshape == (3_000, 256)
     assert es.data.dtype == np.int16
     assert es.electrodes.data == expected_electrode_data_mapping
     assert es.timestamps.shape == (3_000,)
-    assert es.conversion == 0.19499999284744263 * 1e-6
+    assert es.conversion == 1e-6
 
     # Test that the nwbfile has the expected associated files
     assert "associated_files" in nwbfile.processing
@@ -374,6 +374,21 @@ def test_add_raw_ephys_complete_data():
     metadata["ephys"]["targeted_z"] = -2.0  # DV in mm
     metadata["ephys"]["probe"] = ["256-ch Silicon Probe, 3mm length, 66um pitch"]
 
+    # Mapping from probe electrode to channel number for 256-ch silicon probe
+    intan_channel_numbers = np.array([
+        191,190,189,188,187,186,185,184,183,182,181,180,179,178,177,176,175,174,173,172,171,170,169,128,129,130,
+        131,132,133,134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151,152,153,154,155,156,
+        157,158,159,160,161,162,163,164,165,166,167,168,192,193,194,195,196,197,198,199,200,201,202,203,204,205,
+        206,207,208,209,210,211,212,213,214,215,216,217,218,219,220,221,222,223,224,225,226,227,228,229,230,231,
+        232,233,234,235,236,237,238,239,240,241,242,243,244,245,246,247,248,249,250,251,252,253,254,255,64,65,66,
+        67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,
+        101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,
+        127,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,
+        56,57,58,59,60,61,62,63,22,21,20,19,18,17,16,15,14,13,12,11,10,9,8,7,6,5,4,3,2,1,0
+    ])
+    # Reverse mapping from channel (aka row in ElectricalSeries) to row in electrode table
+    expected_electrode_data_mapping = list(np.argsort(intan_channel_numbers))
+
     ephys_data_dict = add_raw_ephys(nwbfile=nwbfile, metadata=metadata)
 
     assert len(nwbfile.electrodes) == 256
@@ -382,13 +397,13 @@ def test_add_raw_ephys_complete_data():
     assert "ElectricalSeries" in nwbfile.acquisition
     es = nwbfile.acquisition["ElectricalSeries"]
     assert es.description == (
-        "Raw ephys data from OpenEphys recording (multiply by conversion factor to get data in volts)."
+        "Raw ephys data from OpenEphys recording, in uV (multiply by conversion factor to get data in V)."
     )
     assert es.data.maxshape == (157_733_308, 256)
     assert es.data.dtype == np.int16
-    assert es.electrodes.data == list(range(256))
+    assert es.electrodes.data == expected_electrode_data_mapping
     assert es.timestamps.shape == (157_733_308,)
-    assert es.conversion == 0.19499999284744263 * 1e-6
+    assert es.conversion == 1e-6
 
     # Test that the nwbfile has the expected associated files
     assert "associated_files" in nwbfile.processing