nctrl-ephys

Overview

nctrl-ephys is a collection of tools for electrophysiology data analysis. It includes functionalities for processing and analyzing data from various sources, including SpikeGLX and Kilosort.

Installation

To install the package, use the following command:

pip install git+https://github.com/nctrl-lab/nctrl-ephys.git

Requirements

• Python 3.9 or higher
• SpikeGLX
• CatGT
• Kilosort4
• phy

Pipelines for SpikeGLX data

1. Run SpikeGLX to record Neuropixels data
2. Run CatGT (ephys catgt) to concatenate and denoise the data
3. Generate the probe map by reading meta file (ephys probe)
   • This will generate a PROBE_TYPE.mat file, which can be useful for the Kilosort GUI.
   • It also plots the probe map and channel numbers.
4. Run Kilosort (ephys runks)
   • Don't forget to switch the conda environment (conda activate kilosort) if the Kilosort is installed in a different environment.
   • If you would like to save metrics (L-ratio, isolation distance, waveform similarity, and overall scores), run ephye runks --metric.
   • If you would like to run Kilosort2, run ephye runks --ks2.
5. Run phy to curate the Kilosort results (conda activate phy2)
6. Save the results (ephys saveks)
   • This will also extract sync and event times data from SpikeGLX files.
   • The final output will be a .mat file that includes the spike, sync, and event times (from the NIDQ file).
   • NIDQ time will be synced with the spike data.
   • If you would like to save all the units, run ephys saveks --all.
   • If you would like to save the good and mua units, run ephys saveks --mua.
   • If you would like to save metrics, run ephys saveks --metric.
7. Extract behavioral data from the VR log file (ephys task)
   • This will extract 'vr', 'trial', 'task_info', 'task_parameter', and 'monitor_info'.
   • The extracted data will be merged with the previous .mat file.
   • The final file will represent one behavioral session.
   • For the avoidance task, run ephys task --pi.
8. Load the data using ephys.spike.Spike class.
   • If you provide electrode location, the module will calculate approximate location of cells and probe channels.
   • spike = Spike(path, coord=[2.0, 4.0, 0.5]) if you inserted the probe at 2.0 mm AP (anterior is positive), 0.5 mm ML (right is positive) from bregma, and 4.0 mm DV (ventral is positive) from pial surface.
     • spike.Spike.channel_position: channel location on the probe (left-lower is the origin)
     • spike.Spike.channel_region: region acronym on the Allen brain atlas
     • spike.Spike.unit_position: unit location on the probe
     • spike.Spike.unit_region: region acronym for each unit
     • spike.Spike.add_region(region_acronym='PL'): draws brain region on the brain atlas
     • spike.Spike.plot_brain(): renders brain atlas including channel and unit location

from ephys.spike import Spike
from ephys.utils import finder

path = finder(msg='Select the .mat file', pattern='.mat$')
spike = Spike(path)
spike

path
    C:\SGL_DATA\abc0\abc0_20240101_M1_g0_imec0\kilosort4\abc0_20240101_M1_g0_imec0_data.mat
spike
    time: (12,)
    frame: (12,)
    firing_rate: (12,)
    position: (12, 2)
    waveform: (12, 61, 14)
    waveform_idx: (12, 14)
    waveform_channel: (12, 14)
    waveform_position: (12, 14, 2)
    Vpp: (12,)
    n_unit: 12
    channel_map: (374,)
    channel_position: (374, 2)
    cluster_group: (370,)
    meta:
    n_channel: 384
    sample_rate: 29999.872727272726
    waveform_raw: (12, 61, 14)
    Vpp_raw: (12,)
sync
    time_imec: (407,)
    frame_imec: (407,)
    type_imec: (407,)
    time_nidq: (407,)
    frame_nidq: (407,)
    type_nidq: (407,)
nidq
    time: (408,)
    frame: (408,)
    chan: (408,)
    type: (408,)
    time_imec: (408,)
vr
    timeSecs: (29738,)
    frame: (29738,)
    timeSecsAfterSplash: (29738,)
    frameAfterSplash: (29738,)
    readTimestampMs: (29738,)
    speed: (29738,)
    rotation: (29738,)
    ballSpeed: (29738,)
    pitch: (29738,)
    roll: (29738,)
    yaw: (29738,)
    distance: (29738,)
    events: (29738,)
    position_x: (29738,)
    position_y: (29738,)
    position_z: (29738,)
trial
    timeSecs: (307,)
    frame: (307,)
    timeSecsAfterSplash: (307,)
    frameAfterSplash: (307,)
    iState: (307,)
    iTrial: (307,)
    iTrial1: (307,)
    iTrial2: (307,)
    iCorrect: (307,)
    iCorrect1: (307,)
    iCorrect2: (307,)
    iChoice: (307,)
    cChoice: (307,)
    iReward: (307,)
    delayDuration: (307,)
    rewardLatency: (307,)
    punishmentLatency: (307,)
    note: (307,)
......

spike.plot() # this will generate an interactive raster and PSTH figure to browse the data.

Pipeline for BMI data

1. Record BMI data
2. Run ephys bmi to merge the BMI data into a binary file to run Kilosort.
3. Run Kilosort (ephys runks)
4. Run phy to curate the Kilosort results and save the result (ephys saveks --bmi)
5. Load the data using ephys.spike.Spike class.

Pipeline for Bruker 2p data

1. Record two-photon imaging with PrairieView. One session folder holds:
   • SESSION.xml: the per-frame table (index, absolute/relative time, file, page, channel) and the microscope state.
   • SESSION.env: the acquisition environment. This is what the module looks for to locate a session.
   • SESSION.companion.ome: OME metadata with the per-frame timestamp and stage position, which the PrairieView XML does not carry.
   • SESSION_CycleNNNNN_ChN_NNNNNN.ome.tif: the images.
   • SESSION_CycleNNNNN_VoltageRecording_NNN.xml and its binary: the analog/TTL channels.
2. Inspect the imaging session (ephys ophys)
   • This prints the microscope settings, the frame table, and the measured frame rate.
   • It counts frames three independent ways (PrairieView XML, companion OME, TIFF pages) and marks a disagreement as MISMATCH, which is the quickest way to catch a truncated or still-writing session.
   • Nothing is written; this step is read-only.
   • Counting TIFF pages directly (Ophys.count_tiff(exact=True)) needs tifffile, which is not installed by default. Without it, the count comes from the companion OME.
3. Extract the voltage recording and the behavior (ephys vrec)
   • This turns every channel into on/off pulse times, then parses the trial structure of the VR task (--task vr, the default).
   • All times are realigned so that the first imaging frame onset (AI 1) is t = 0, which puts them on the same clock as the PrairieView frame table.
   • The result is saved to FOLDER_data.mat under a single vrec struct, so nothing collides with the trial that ephys task writes into the same file. Re-running replaces vrec after asking to overwrite.
     • vrec.channels: the channel names, as a cell array.
     • vrec.data: a cell array in the same order, each cell an (n_pulse, 2) array of on/off times in seconds. A channel that never went high gets an empty array, so the two stay index-aligned and vrec.data{c} always belongs to vrec.channels{c}.
     • vrec.trial: the parsed trials, absent when the task is not parsed.
   • The TTL channel map the VR task expects:
     • AI 1: 2p imaging frame, high once per frame
     • AI 2: task start (on) and end (off)
     • AI 3: delay start (on), cue start (off)
     • AI 4: choice/ITI start (on), next trial delay start (off)
     • AI 5: direction, read during the cue and the choice window (1: left, 2: right)
     • AI 6: water reward
   • Crosstalk puts sub-millisecond pulses on otherwise idle lines, and each one would be counted as a phantom trial or reward. Pulses shorter than jitter (2 ms by default) are dropped, and the dropped count is reported per channel. Real TTLs here are 31 ms and longer, so the margin is wide; retune per rig with VRec.parse_data(jitter=...) and then re-run parse_task().
4. Run suite2p on the images, then collect its output (ephys s2p)
   • This reads F, Fneu, spks, stat, iscell and ops from a suite2p/planeN folder, and adds the neuropil-corrected trace (F - neucoeff * Fneu, --neucoeff 0.7 by default) and dF/F over suite2p's own maximin baseline.
   • Everything is saved to the same FOLDER_data.mat under a suite2p struct, next to vrec. ops keeps the settings, the summary images and the rigid offsets; the registration QC (regPC alone is ~57 MB) is left out.
   • Check the frame rate suite2p ran with. It is whatever was typed into the GUI, and it sets the deconvolution kernel and the baseline window. Compare it with the rate ephys ophys reports and pass --fs to override it; the override is printed in red so a wrong rate does not pass unnoticed. F and Fneu do not depend on it, but spks does.
   • dF/F is left as nan for an ROI whose baseline is not positive, which happens when the ROI is no brighter than the neuropil around it. Read corrected for those instead of trusting a ratio to a near-zero baseline.
5. Load the data using the ephys.ophys.Ophys, ephys.ophys.VRec and ephys.ophys.Suite2p classes.
   • VRec.data is indexed like VRec.channels; VRec.pulse('AI 3') looks a line up by name instead.
   • VRec.frame2time(idx) gives the time of an imaging frame, and VRec.time2frame(t) gives the frame a task event falls in. Use these to put imaging and behavior on one index.
   • Suite2p.cell is the boolean mask of the ROIs the classifier accepted, so s2p.dff[s2p.cell] is the traces you usually want.

from ephys.ophys import Ophys, VRec, Suite2p

path = '/path/to/20260711_DIT09_4X-136'
ophys = Ophys(path)
print(ophys)

vrec = VRec(path, task='vr')
print(vrec)
vrec.save()

s2p = Suite2p(path, fs=ophys.frame_rate) # suite2p ran at the wrong rate here
print(s2p)
s2p.save()

ophys
    Ophys: 20260711_DIT09_4X-136.xml
        laserPower:
            Uncaging Laser Power Modulator: 0
            Imaging Laser Power Modulator: 250
        pmtGain:
            Detector 1: 748.192749
            Detector 2: 900
        framePeriod: 0.033182
        rastersPerFrame: 2
        linesPerFrame: 512
        pixelsPerLine: 512
        opticalZoom: 4
    frames: (35004, 7) ['cycle', 'index', 'absoluteTime', 'relativeTime', 'filename', 'page', 'channel']
        1 cycle(s), 1 channel(s), 15.07 Hz
    ome: (35004, 9) ['t', 'z', 'c', 'deltaTime', 'positionX', 'positionY', 'positionZ', 'filename', 'ifd']
    n_frame: xml 35004, ome 35004, tiff 35004
vrec
    VRec: 20260711_DIT09_4X-136_Cycle00001_VoltageRecording_001.xml
    23230464 samples @ 10000 Hz = 2323.0 s
        AI 1: 35006 pulses, duration: 0.066 s
        AI 2: 1 pulses, duration: 2315.481 s
        AI 3: 126 pulses, duration: 4.603 s
        AI 4: 125 pulses, duration: 8.164 s
        AI 5: 83 pulses, duration: 12.063 s
        AI 6: 89 pulses, duration: 0.065 s
        AI 7: 0 pulses, duration: 0.000 s
        AI 8: 0 pulses, duration: 0.000 s
    trial: 125 trials, 89 rewarded
        ['nTrial', 'timeStart', 'timeCue', 'timeChoice', 'timeEnd', 'cue', 'choice', 'result']
s2p
    Suite2p: plane0
    102 ROIs (45 cells), 35004 frames @ 15.0685 Hz = 2323.0 s
        F: (102, 35004)
        Fneu: (102, 35004)
        spks: (102, 35004)
        corrected: (102, 35004)
        dff: (102, 35004) (undefined for 66 ROIs)
        neucoeff: 0.7, tau: 1
        fs overridden: 10 -> 15.0685 Hz

• AI 1 usually carries a few more pulses than there are saved frames: PrairieView keeps the frame clock running for a moment past the last frame it writes. The extra pulses sit at the end, so frame indices still line up.
• The binary voltage file can be longer than SamplesAcquired in its XML. The trailing rows are a stale acquisition buffer, not data, and are dropped on load; the sample count in the XML is authoritative.

Usage

Command Line Interface

The package provides a command-line interface (CLI) for various operations. Below are some examples:

Running CatGT

ephys catgt --path /path/to/data

• You can omit the --path option and the command will ask you the path to the data.

Running Kilosort

ephys runks --path /path/to/data

• You can omit the --path option and the command will ask you the path to the data.

Saving Kilosort Results

• This command will generate a '.mat' file containing only the good units that were curated by phy.
• This command also saves the waveform data by reading the raw '.bin' files (by default, it will read the first 60 seconds of data).

ephys saveks --path /path/to/data

Python API

You can also use the functionalities provided by nctrl-ephys directly in your Python scripts. Check the example.ipynb for more details.

Reading and plotting SpikeGLX Data

from ephys.spikeglx import read_meta, read_analog, read_digital
from ephys.utils import finder

# Finding the data file
fn = finder("C:\\SGL_DATA")

# Loading the meta data
meta = read_meta(fn)

# Loading the Neuropixels data
data = read_analog(fn, sample_range=(0, 3000))

# Plotting the Neuropixels data
plt.imshow(data.T, vmin=-200, vmax=200, cmap='bwr', aspect='auto', interpolation='none')
plt.colorbar()
plt.title('Raw Data')
plt.xlabel('Time (samples)')
plt.ylabel('Channel')
plt.show()

# Loading digital data to get the sync pulse times
data_event = read_digital(fn)
time_sync = data_event.query('chan == 6').times.values

Running Kilosort

from ephys.ks import run_ks4
run_ks4(path='/path/to/data')

Loading and Plotting MUA Data

from ephys.bmi import BMI

bmi = BMI(path='/path/to/data')
bmi.load_mua()
bmi.plot_mua()

Ploting raster and PSTH

# load mat file
path = finder(path="C:\SGL_DATA", msg='Select a session file', pattern=r'.mat$')
spike = Spike(path)

# plot raster and psth
time_spike = spike.spike['time'][0]
time_event = spike.nidq.query('chan == 2 and type == 1')['time_imec'].values
plot_raster_psth(time_spike, time_event)
plt.show()

License

This project is licensed under the MIT License. See the LICENSE file for details.