From affca4975d44d0cdf5e8c8723899423d4312d0ab Mon Sep 17 00:00:00 2001 From: Jude Werth <143547832+judewerth@users.noreply.github.com> Date: Mon, 2 Mar 2026 11:55:37 -0700 Subject: [PATCH 1/2] 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 --- element_array_ephys/ephys_no_curation.py | 210 ++++++++++++++++++++++- 1 file changed, 204 insertions(+), 6 deletions(-) diff --git a/element_array_ephys/ephys_no_curation.py b/element_array_ephys/ephys_no_curation.py index 2afc7749..38845b65 100644 --- a/element_array_ephys/ephys_no_curation.py +++ b/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, + } + ) \ No newline at end of file From d2e2f086b230ff88298a83edfb2d968ee30118f5 Mon Sep 17 00:00:00 2001 From: Jude Werth <143547832+judewerth@users.noreply.github.com> Date: Wed, 25 Mar 2026 10:39:49 -0600 Subject: [PATCH 2/2] added get_probe_type function for clarity --- element_array_ephys/ephys_no_curation.py | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/element_array_ephys/ephys_no_curation.py b/element_array_ephys/ephys_no_curation.py index 38845b65..a7a046bb 100644 --- a/element_array_ephys/ephys_no_curation.py +++ b/element_array_ephys/ephys_no_curation.py @@ -148,6 +148,17 @@ def map_channel_to_electrode(probe_type="A1x32-6mm-100-177-H32_21mm", input_indi electrode_ids = lookup[input_indices] return electrode_ids + +def get_probe_type(ephys_key): + """ + Gets the probe type for a given ephys session key. EphysSessionProbe needs an entry along with the EphysSession for ephys_key + """ + probe_type = set((EphysSessionProbe * probe.Probe & ephys_key).fetch('probe_type')) + if len(probe_type) != 1: + raise ValueError( + f"Couldn't identify probe type for {ephys_key} - expected one, found {len(probe_type)}" + ) + return probe_type.pop() # ----------------------------- Table declarations ----------------------