Add trial number to unit

[ahabdallah27]

When I run my experiments I run a given trial for some time, I move on to another trial. In a given series of tests on any given day, I may have 1000's trials. I am also using a Utah Array so I have 96 channels that I am recording data from. Therefore for each trial I have 96 spike times. Is it possible for me to add more than one channel for a row in the unit table? By doing it this way, the unit table would have only 1 row per trial (so 1000 rows, for example). Or am I required to make a unique row for each row and each trial resulting in 96*1000 rows (96000 rows in this case). The code below is an example with 3 channels, and I get the error shown in the screenshot when I try to append the spike times together.

If neither approach is recommended, please do let me know the best way to set this up.

from datetime import datetime
from dateutil.tz import tzlocal
from uuid import uuid4
import numpy as np
import pandas
from pynwb import NWBHDF5IO, NWBFile
from pynwb.file import Subject
from pynwb.device import Device
from pynwb.ecephys import LFP, ElectricalSeries, TimeSeries, ElectrodeGroup


subject = Subject(description='Primate',
                   subject_id='P')

device = Device(name="Utah Array",
                description='96 channel Utah Array',
                manufacturer="Blackrock Neurotech")

electrodegroup = ElectrodeGroup(name='BMI',
                                description='Electrodes',
                                device=device,
                                location='Brain')

nwbfile = NWBFile(session_description='Lab',
                  identifier='identifier',
                  session_start_time=datetime.now(tzlocal()), #these three are required
                  file_create_date=datetime.now(tzlocal()),
                  experimenter=['Person 1'],
                  lab='My Lab',
                  subject=subject,
                  devices=[device],
                  electrode_groups=[electrodegroup])

nwbfile.add_electrode_column(name="label", description="label of electrode")

channels = 96
for i in range(channels):
    nwbfile.add_electrode(
            group=electrodegroup,
            label="Channel {}".format(i),
            location="Brain",
        )
    
all_table_region = nwbfile.create_electrode_table_region(
    region=list(range(channels)),  # reference row indices 0 to N-1
    description="all electrodes",
)

nwbfile.add_trial_column(name='Mask',
                         description='Mask applied to data')
nwbfile.add_trial_column(name='TrialNumber',
                         description='Trial')

nwbfile.add_unit_column(name='TrialNumber',
                         description='Trial')

#                        trial 1                    trial 2                   trial 3                      
trials = {'time':        [list(range(101, 201, 1)), list(range(225, 265, 1)), list(range(450, 812, 1))], 
          'TrialNumber': [5,6,7],
        'Mask':          [[0,0,1,0,1],              [0,0,1,0,1],              [0,0,1,0,1]],
#                     Using only two channels as an example for spike times
        'SpikeTimes':    [[[],[]],             [[230, 238, 260], [244]], [[],[557, 697]]]
        }
for idx in range(len(trials['time'])):
    nwbfile.add_trial(start_time=float(trials['time'][idx][0]),
                    stop_time=float(trials['time'][idx][-1]),
                    timeseries=TimeSeries(
                                        name='Trial ' + str(trials['TrialNumber'][idx]),
                                        data=trials['time'][idx],
                                        timestamps=trials['time'][idx],
                                        unit="s",
                                        conversion=.001
                                        ),
                    Mask=trials['Mask'][idx],
                    TrialNumber=trials['TrialNumber'][idx])
    spikeTimes = []
    for elems in trials['SpikeTimes'][idx]:
        spikeTimes.append(elems)
    nwbfile.add_unit(spike_times=spikeTimes,
                     electrode_group=nwbfile.get_electrode_group(name='BMI'),
                     TrialNumber=trials['TrialNumber'][idx])
    

[CodyCBakerPhD]

Hi @ahabdallah27

What I would strongly suggest from a data modeling point of view, is to decouple the trials and units tables to be separate

A trials table is native to the NWBFile object, check out this tutorial for an easy way to add trial times (relative to the global session_start_time)

Likewise, the units table is native to the NWBFile object, but all it does is specify the times when spikes occurred (again, with respect to the global session_start_time) and associate a single units spike train with additional info, including the electrodes responsible for detecting the unit

And yes, this can accept multiple electrode indices! For example,

nwbfile.add_unit(spike_times=[1,2,3], electrodes=[5,6,7])

(notice these indices don't even need to be explicit DynamicTableRegions - the helper function .add_unit will take care of that for you)

Also, for the location column of the electrodes table, that is intended to be the specific subregion of the brain you are recording from, preferable using a common abbreviation framework (an example is Allen CCF for mice); for example, "VISrll" for "Rostrolateral lateral visual area". The full name, if an abbreviation is not standard or typical, is also fine. If you don't know the specific region, then setting it to "unknown" is preferred

[ahabdallah27]

With that structure, how do I easily get all the units from a single or select number of trials? My understanding is that I would have to iterate through every row the individual electrode and then search for the timestep related to the trial and that append each one of those together to see trends across trials.

[CodyCBakerPhD]

how do I easily get all the units from a single or select number of trials?

If you mean 'spike times' of a single unit (a single row of units table) within a given trial, yes you would filter the train based on the trial times

An easy way to do this with numpy is to use searchsorted; here's a basic demo

import numpy

spike_times_of_unit = numpy.array([1.1,1.2,1.3,1.4,1.5])
single_trial_interval = [1.15, 1.45]

spiking_within_trial_start, spiking_within_trial_stop = numpy.searchsorted(
    spike_times_of_a_unit, a_single_trial_interval
)
> (1, 4)  # indices of the spikes that occur during the trial

spike_times_of_unit_within_trial = spike_times_of_a_unit[spiking_within_trial_start:spiking_within_trial_stop]
> array([1.2, 1.3, 1.4])  # Actual subset of spikes during the trial

A huge advantage of storing data this way is how it generalizes to multiple data streams within the file; by not trializing the source storage the only task to achieve for each individual stream is to ensure it is synchronized to the common session_start_time; then we can easily generate whatever alignment types we need for analysis

It also makes it easier to add hyperparameter shifts around those alignment points, such as a similar data view for events occurring some period of time before and after each trial

Example of existing visualizers for aligned spiking activity: https://flatironinstitute.github.io/neurosift/?p=/nwb&url=https://api.dandiarchive.org/api/assets/2abd0c2b-1190-4e79-a47c-9636e8ec4160/download/&dandisetId=000409&dandisetVersion=draft&tab=view:PSTH|/intervals/trials

And note in that file, there are many additional data streams for behavior too, so you can make joint plots such as this one: https://flatironinstitute.github.io/neurosift/?p=/nwb&url=https://api.dandiarchive.org/api/assets/2abd0c2b-1190-4e79-a47c-9636e8ec4160/download/&dandisetId=000409&dandisetVersion=draft&tab=neurodata-items:neurodata-item:/intervals/trials|TimeIntervals@neurodata-item:/processing/behavior/WheelVelocity|TimeSeries@neurodata-item:/processing/behavior/WheelAcceleration|TimeSeries&tab-time=1467.8397655900355,1476.8397655900355,1470.336463080128 (trial regions annotated on top)

[ahabdallah27]

I was expecting a lot of helper functions/tutorials for common things that are done with nwb files. For example, binning the data into longer time steps or removing data across the timeseries/trials based on certain criteria. I don't see things like that in the tutorials on this page. Do you know where I'd find them or do I have to make my own helper functions? I recognize things can take a couple lines like you showed above, but I'd like to turn those and longer functions like the two I mentioned into helper functions if that does not already exist.

[ahabdallah27]

Understood, thank you, I appreciate it!