Changes to Ephys Schema

1) map channel to electrode function
2) Fix bug in trace extraction (currently doesn't account for file boundaries
3) Add impedance measurement tables
4) Add STTC tables

Author: judewerth

element_array_ephys/ephys_no_curation.py

@@ -9,6 +9,9 @@
 from element_interface.utils import dict_to_uuid, find_full_path, find_root_directory
 from scipy import signal
 import intanrhdreader
+import neo
+import quantities as pq
+from elephant.spike_train_correlation import spike_time_tiling_coefficient
 
 from . import ephys_report, probe
 from .readers import kilosort, openephys, spikeglx
@@ -115,7 +118,36 @@ def get_processed_root_data_dir() -> str:
     else:
         return get_ephys_root_data_dir()[0]
 
+def map_channel_to_electrode(probe_type="A1x32-6mm-100-177-H32_21mm", input_indices=None, electrode_to_channel=False):
+    """
+    Maps channel indices from recording controller to specific probe geometry as defined in probe.ElectrodeConfig.Electrode.
+
+    Args:
+        probe_type (str): Name of the probe used in the recording session. See probe.ProbeType() for inserted probes.
+        electrode_to_channel (bool): If True, maps from electrode indices to channel indices. If False, maps from channel indices to electrode indices. Default is False.
+    
+    Returns:
+        electrodes (array-like): Array of electrode indices corresponding to the input channel indices. 
+            If electrode_to_channel is False, the output will be electrode indices. If electrode_to_channel is True, the output will be channel indices.
+    """
 
+    # get electrode and channel info
+    num_electrodes = len(probe.ProbeType.Electrode & f"probe_type='{probe_type}'")
+    electrode_mapping, channel_mapping = probe.ElectrodeConfig.Electrode.fetch("electrode", "channel_idx")
+
+    # create lookup to convert
+    lookup = np.empty(num_electrodes, dtype=int)
+    if electrode_to_channel:
+        lookup[electrode_mapping] = channel_mapping
+    else:
+        lookup[channel_mapping] = electrode_mapping
+
+    # correctly map electrode indices
+    if input_indices is None:
+        input_indices = np.arange(num_electrodes)
+
+    electrode_ids = lookup[input_indices]
+    return electrode_ids
 # ----------------------------- Table declarations ----------------------
 
 
@@ -152,7 +184,6 @@ class EphysRawFile(dj.Manual):
     filename_prefix   : varchar(64)  #  filename prefix, if any, excluding the datetime information
     """
 
-
 @schema
 class EphysSession(dj.Manual):
     definition = """ # User defined ephys session for downstream analysis.
@@ -224,7 +255,6 @@ def make(self, key):
             ]
         )
 
-
 @schema
 class LFP(dj.Imported):
     definition = """ # Store pre-processed LFP traces per electrode. Only the LFPs collected from a pre-defined recording session.
@@ -384,6 +414,28 @@ def make_compute(
                 }
 
             lfps = data.pop("amplifier_data")[lfp_indices]
+
+            # account for boundaries
+            fs = header["sample_rate"]
+            if file_relpath == file_paths[0]:
+                file_start = datetime.strptime(
+                    "_".join(file_relpath.split("_")[3:5]).removesuffix(".rhd"),
+                    "%y%m%d_%H%M%S",
+                )
+                start_idx = int((key['start_time'] - file_start).total_seconds() * fs)
+
+                # trim lfps to start boundary
+                lfps = lfps[:, start_idx:]
+            elif file_relpath == file_paths[-1]:
+                file_start = datetime.strptime(
+                    "_".join(file_relpath.split("_")[3:5]).removesuffix(".rhd"),
+                    "%y%m%d_%H%M%S",
+                )
+                end_idx = int((key['end_time'] - file_start).total_seconds() * fs)
+
+                # trim lfps to end boundary
+                lfps = lfps[:, :end_idx]            
+
             lfp_concat.append(lfps)
 
         full_lfp = np.hstack(lfp_concat)
@@ -408,10 +460,11 @@ def make_compute(
             # Downsample the signal with `decimate`
             lfp = signal.decimate(lfp, downsample_factor, ftype="fir", zero_phase=True)
             all_lfps.append(lfp)
-
-        execution_duration = (
-            datetime.now(timezone.utc) - execution_time
-        ).total_seconds() / 3600
+            
+        execution_duration = ((
+                    datetime.now(timezone.utc) - execution_time
+                ).total_seconds()
+                / 3600)
         return (
             all_lfps,
             channels,
@@ -448,7 +501,89 @@ def make_insert(
                 }
             )
 
+@schema
+class ImpedanceFile(dj.Manual):
+    definition = """ # Insert files and organoid_id for impedance measurements
+    -> ephys.EphysRawFile
+    organoid_id        : varchar(4) # e.g. O17
+    """
+
+@schema
+class ImpedanceMeasurements(dj.Imported):
+    definition = """ # Store impedance measurements per channel
+    -> ImpedanceFile
+    ---
+    port_id: char(2)  # Port ID of the Intan acquisition system
+    """
+
+    class Channel(dj.Part):
+        definition = """
+        -> master
+        channel_idx: int  # channel index
+        ---
+        channel_id: varchar(64)  # channel id
+        impedance_magnitude: float  # in Ohms
+        impedance_phase: float  # in Degrees
+        """
 
+    def make(self, key):
+        # fetch file path from ephysrawfile entry
+        file_path = (EphysRawFile & key).fetch1("file_path")
+
+        # import file
+        file = find_full_path(get_ephys_root_data_dir(), file_path)
+        try:
+            data = intanrhdreader.load_file(file)
+        except OSError:
+            raise OSError(f"OS error occurred when loading file {file.name}")
+
+        # extract amplifier channels
+        amplifier_channels = data['header'].pop("amplifier_channels")
+
+        # Figure out `Port ID` from the existing EphysSessionProbe
+        port_id = set((EphysSessionProbe & key).fetch("port_id"))
+
+        # Figure out `Port ID` from the existing EphysSession
+        if not (EphysSessionProbe & key):
+            raise ValueError(
+                f"No EphysSessionProbe found for the {key} - cannot determine the port ID"
+            )
+
+        # Check if there are multiple port IDs for the same experiment, if so, it needs to be fixed in the EphysSessionProbe table
+        if len(port_id) > 1:
+            raise ValueError(
+                f"Multiple Port IDs found for the {key} - cannot determine the port ID"
+            )
+        port_id = port_id.pop()
+
+        # get channels for the correct port
+        port_channels = [channel for channel in amplifier_channels if channel['port_prefix'] == port_id]
+
+        # insert into master
+        self.insert1(
+            {
+                **key,
+                "port_id": port_id,
+            }
+        )
+
+        # loop through channels and insert impedance data
+        for channel in port_channels:
+
+            channel_idx = channel['custom_order']
+            channel_id = channel['custom_channel_name']
+            impedance_magnitude = channel['electrode_impedance_magnitude']
+            impedance_phase = channel['electrode_impedance_phase']
+
+            self.Channel.insert1(
+                {
+                    **key,
+                    "channel_idx": channel_idx,
+                    "channel_id": channel_id,
+                    "impedance_magnitude": impedance_magnitude,
+                    "impedance_phase": impedance_phase,
+                }
+            )
 # ------------ Clustering --------------
 
 
@@ -1101,3 +1236,66 @@ def make(self, key):
         self.insert1(key)
         self.Cluster.insert(metrics_list, ignore_extra_fields=True)
         self.Waveform.insert(metrics_list, ignore_extra_fields=True)
+
+"""
+Functional Connectivity (STTC)
+"""
+@schema
+class STTC(dj.Computed):
+    """
+    Spike Time Tiling Coefficient (STTC) between unit pairs. Automatically computed within ephys sessions (spike sorting).
+    Based on the method described in Sharf et al. (2022) Nature Communications.
+    """
+
+    definition = """
+    -> ephys.CuratedClustering
+    unit_a: int  # First unit in the pair
+    unit_b: int  # Second unit in the pair
+    ---
+    sttc: float  # STTC value between unit pairs
+    spike_time_latencies: longblob  # Latencies (ms) of spikes from unit A to nearest spike in unit B during (limited to +/- dt) 
+    """
+
+    def make(self, key):
+
+        # define parameters
+        dt = 20 # ms    
+
+        # fetch spike times for all units in the clustering
+        unit_ids, spike_times = (CuratedClustering.Unit & key).fetch('unit', 'spike_times', order_by='unit')
+
+        num_units = len(unit_ids)
+        t_stop = (key['end_time'] - key['start_time']) / timedelta(milliseconds=1)  # in ms
+
+        # REMOVE LATER
+        spike_times = np.array([st[st <= (key['end_time'] - key['start_time']).total_seconds()] for st in spike_times], dtype=object)
+
+        # loop through unit pairs and calculate STTC
+        for i in range(num_units - 1):
+            for j in range(i + 1, num_units):
+
+                # get spike times (convert from seconds to miliseconds)
+                spikes_A = (spike_times[i] * (timedelta(seconds=1) / timedelta(milliseconds=1))).astype(int)  
+                spikes_B = (spike_times[j] * (timedelta(seconds=1) / timedelta(milliseconds=1))).astype(int)  
+
+                # convert to spike trains (neo)
+                spiketrain_A = neo.SpikeTrain(spikes_A, units='ms', t_stop=t_stop)
+                spiketrain_B = neo.SpikeTrain(spikes_B, units='ms', t_stop=t_stop)
+
+                # calculate STTC
+                sttc = spike_time_tiling_coefficient(spiketrain_A, spiketrain_B, dt=dt*pq.ms)
+
+                # calculate spike time latencies
+                diff_matrix = np.abs(np.subtract.outer(spikes_A, spikes_B))
+                closest_spikes = np.min(diff_matrix, axis=1) # closest spike in B for each spike in A
+                spike_time_latencies = closest_spikes[closest_spikes <= dt]
+
+                self.insert1(
+                    {
+                        **key,
+                        'unit_a': unit_ids[i],
+                        'unit_b': unit_ids[j],
+                        'sttc': sttc,
+                        'spike_time_latencies': spike_time_latencies,
+                    }
+                )