Merge branch 'main' of https://whale.am28.uni-tuebingen.de/git/awendt/oephys2nix

oephys2nix/main.py

@@ -129,6 +129,7 @@ def convert(
 
         log.info("Starting with stimulus recreation")
         stim = StimulusToNix(open_ephys, str(relacs), str(nix_path))
+        stim.open()
         stim.create_repros_automatically()
         stim.print_table()
         stim.close()

oephys2nix/sorting.py

@@ -3,11 +3,12 @@ import pathlib
 
 import nixio
 import numpy as np
-import spikeinterface.core as si
-import spikeinterface.extractors as extractors
+import spikeinterface.core as sc
+import spikeinterface.extractors as se
 from IPython import embed
 from nixio.exceptions import DuplicateName
 from rich.console import Console
+from scipy.signal import butter, sosfiltfilt
 
 log = logging.getLogger(__name__)
 console = Console()
@@ -30,29 +31,82 @@ class AppendSorting:
     def __init__(self, sorter_path: pathlib.Path, recording_path: pathlib.Path):
         self.sorter_path = sorter_path
         self.recording_path = recording_path
-        self.sorting = extractors.read_phy(self.sorter_path)
+        self.sorting = se.read_phy(self.sorter_path)
         self.nixfile = nixio.File.open(str(self.recording_path), nixio.FileMode.ReadWrite)
         self.block = self.nixfile.blocks[0]
         self.das = self.block.data_arrays
-        self.channel_ids = si.get_template_extremum_channel(
-            self.sorting, mode="extremum", peak_sign="neg", outputs="index"
-        )
+        self.channel_ids = self.sorting.get_property("ch")
         self.data = self.block.data_arrays["data"]
+        self.analyzer = sc.load_sorting_analyzer(self.recording_path.parent / "sorting_analyzser")
         self._clean()
+        self.append_processed_recording()
         self.append_sorting_to_recording()
 
+    def bandpass_filter(
+        self, data: np.ndarray, fs: float, lowcut: float, highcut: float, order: int = 5
+    ) -> np.ndarray:
+        coeffs = butter(
+            order,
+            [lowcut, highcut],
+            btype="band",
+            analog=False,
+            fs=fs,
+            output="sos",
+        )
+
+        return sosfiltfilt(sos=coeffs, x=data, axis=0).astype(np.float32)
+
+    def common_median_reference(self, data: np.ndarray) -> np.ndarray:
+        return (data - np.median(data, axis=1, keepdims=True)).astype(np.float32)
+
+    def append_processed_recording(self):
+
+        log.debug("Appending raw data")
+        processed_data = np.zeros_like(self.data, dtype=np.float32)
+        data_ref = self.common_median_reference(self.data)
+        for ch in range(data_ref.shape[1]):
+            processed_data[:, ch] = self.bandpass_filter(
+                data_ref[:, ch], self.sorting.sampling_frequency, 250, 5000
+            )
+
+        gr = self.block.groups["neuronal-data"]
+        try:
+            nix_data_array = self.block.create_data_array(
+                name="data_preprocessed",
+                array_type="open-ephys.data.sampled",
+                dtype=nixio.DataType.Float,
+                data=processed_data,
+                unit="uV",
+            )
+            nix_data_array.append_sampled_dimension(
+                self.data.dimensions[0].sampling_interval, label="time", unit="s"
+            )
+            nix_data_array.append_sampled_dimension(1, label="channel")
+        except DuplicateName:
+            del self.block.data_arrays["data_preprocessed"]
+            nix_data_array = self.block.create_data_array(
+                name="data_preprocessed",
+                array_type="open-ephys.data.sampled",
+                dtype=nixio.DataType.Float,
+                data=processed_data,
+                unit="uV",
+            )
+            nix_data_array.append_sampled_dimension(
+                self.data.dimensions[0].sampling_interval, label="time", unit="s"
+            )
+            nix_data_array.append_sampled_dimension(1, label="channel")
+
+        gr.data_arrays.append(nix_data_array)
+
     def append_sorting_to_recording(self):
         try:
             gr = self.block.create_group("units", "sorting.group")
         except DuplicateName:
             del self.block.groups["units"]
             gr = self.block.create_group("units", "sorting.group")
-        for unit in self.sorting.unit_ids:
-            spike_times = self.sorting.sorting.get_unit_spike_train_in_seconds(
-                unit, segment_index=0
-            )
-            unit_channel = self.channel_ids[unit]
-            channel_tag = np.repeat(unit_channel, spike_times.shape[0])
+        for i, unit in enumerate(self.sorting.unit_ids):
+            spike_times = self.sorting.get_unit_spike_train_in_seconds(unit, segment_index=0)
+            channel_tag = np.repeat(self.channel_ids[i], spike_times.shape[0])
             multi_tag_positions = np.column_stack((spike_times, channel_tag))
 
             try:
@@ -73,6 +127,7 @@ class AppendSorting:
                     f"unit-{unit}", "sorting.spike_index", multi_tag_positions
                 )
                 multi_tag.references.append(self.data)
+                multi_tag.references.append(self.block.data_arrays["data_preprocessed"])
             except DuplicateName:
                 del self.block.multi_tags[f"unit-{unit}"]
                 del self.das[f"unit-{unit}-positions"]
@@ -80,18 +135,18 @@ class AppendSorting:
                     f"unit-{unit}", "sorting.spike_index", multi_tag_positions
                 )
                 multi_tag.references.append(self.data)
-            gr.mulit_tags.append(positions)
+            gr.multi_tags.append(multi_tag)
 
     def _clean(self):
         try:
             gr = self.block.groups["units"]
-        except IndexError:
+        except KeyError:
             return
 
         for das in gr.data_arrays:
             del self.das[das.name]
         for mtag in gr.multi_tags:
-            del self.block.mulit_tags[mtag.name]
+            del self.block.multi_tags[mtag.name]
             del self.das[mtag.name + "-positions"]
 
     def close(self):

oephys2nix/stimulus_recreation.py

@@ -63,26 +63,32 @@ class StimulusToNix:
     """
 
     def __init__(self, open_ephys_path: pathlib.Path, relacs_nix_path: str, nix_file: str):
+        self.open_ephys_path = open_ephys_path
         self.relacs_nix_path = relacs_nix_path
         self.nix_file_path = nix_file
 
-        self.dataset = rlx.Dataset(relacs_nix_path)
-
-        self.relacs_nix_file = nixio.File.open(relacs_nix_path, nixio.FileMode.ReadOnly)
-        self.relacs_block = self.relacs_nix_file.blocks[0]
-        self.relacs_sections = self.relacs_nix_file.sections
-
-        self.neo_data = OpenEphysBinaryIO(open_ephys_path).read(lazy=True)
-        self.fs = self.neo_data[0].segments[0].analogsignals[0].sampling_rate.magnitude
-
-        self.nix_file = nixio.File(nix_file, nixio.FileMode.ReadWrite)
-        self.block = self.nix_file.blocks[0]
-
         # Threshold for TTL peak
         # constants
         self.threshold = 2
         self.new_start_jiggle = 0.1
 
+    def open(self):
+        self.relacs_nix_file = nixio.File.open(self.relacs_nix_path, nixio.FileMode.ReadOnly)
+        self.relacs_block = self.relacs_nix_file.blocks[0]
+        self.relacs_sections = self.relacs_nix_file.sections
+
+        self.neo_data = OpenEphysBinaryIO(self.open_ephys_path).read(lazy=True)
+        self.fs = self.neo_data[0].segments[0].analogsignals[0].sampling_rate.magnitude
+
+        self.nix_file = nixio.File(self.nix_file_path, nixio.FileMode.ReadWrite)
+        self.block = self.nix_file.blocks[0]
+
+        self.dataset = rlx.Dataset(self.relacs_nix_path)
+
+    def open_nix_file(self):
+        self.nix_file = nixio.File(self.nix_file_path, nixio.FileMode.ReadWrite)
+        self.block = self.nix_file.blocks[0]
+
     def _append_relacs_tag_mtags(self) -> None:
         """Append relacs tags and multi tags to new nix file."""
         for t in self.relacs_block.tags:
@@ -456,6 +462,9 @@ class StimulusToNix:
                     last_repro_position.reshape(-1, 1),
                     (current_position - last_repro_position).reshape(-1, 1),
                 )
+            self.close_nix_file()
+            log.debug("Flush the nix file")
+            self.open_nix_file()
 
     def create_repros_from_config_file(self) -> None:
         """Creates repros form a config file.
@@ -719,6 +728,9 @@ class StimulusToNix:
         plt.legend(loc="upper right")
         plt.show()
 
+    def close_nix_file(self):
+        self.nix_file.close()
+
     def close(self):
         self.dataset.close()
         self.nix_file.close()