Tutorial: from a file to a raster

One real session, start to finish, with the intermediate output shown at each step. Runtimes are from a four-core laptop with the file on an internal SSD; they are there so you know whether a step is slow or stuck.

The data#

A ten-minute recording, 32 channels of int16 interleaved by sample, 30 kHz, 1.15 GB on disk as session12.dat. Beside it sits a sidecar describing the layout, which is what spares you from passing six keyword arguments every time:

{
  "fs": 30000,
  "n_channels": 32,
  "dtype": "int16",
  "gain_uv_per_bit": 0.195,
  "channel_groups": [[0,1,2,3], [4,5,6,7]],
  "bad_channels": [11, 26]
}

With the sidecar in place the reader takes the path alone:

import axonforge as af

rec = af.io.read("session12.dat")
print(rec.summary())

Recording  session12.dat
  32 channels, 2 marked bad (11, 26)
  18 000 000 samples at 30 000 Hz  =  600.0 s
  int16, 0.195 uV/bit, memory-mapped
  8 channel groups of 4

Referencing and filtering#

Common median referencing subtracts, sample by sample, the median across the channels of each group. It removes movement and chewing artefacts that appear on all channels at once, and unlike a common average it is not dragged around by one loud channel.

rec = af.filt.common_median(rec)              # 11 s
rec = af.filt.bandpass(rec, low=300, high=6000, order=4)
sigma = af.filt.noise_sigma(rec)              # robust estimate, per channel
print(sigma.round(1)[:8])                     # [ 7.4  7.1  7.9  7.2  6.8 14.6  7.5  7.3]

Channel 5 stands out at 14.6 µV against a 7.4 µV median. That is the kind of thing worth looking at before detection rather than after sorting: a doubled noise floor halves the effective threshold in units of amplitude, and the cluster it produces will look like a real unit with a suspiciously wide waveform.

Detection#

spikes = af.detect.threshold(rec, k=4.5, polarity="neg", dead_time_ms=1.0)
print(spikes)                                 # 38 s

SpikeTable  41 882 events on 30 channels (2 excluded)
  rate per channel: median 2.1 /s, min 0.3, max 9.8
  threshold: 4.5 x sigma, negative-going
  dead time: 1.0 ms enforced after each crossing

Nine point eight events per second on one channel against a median of two is not necessarily wrong — a well-isolated, fast-spiking cell does that — but it is worth checking that it is not the channel with the raised noise floor. Here it was channel 19, not 5.

Waveforms and alignment#

wf = af.waveform.extract(rec, spikes, pre_ms=0.5, post_ms=1.5)
print(wf.shape)                               # (41882, 60)  = 2.0 ms at 30 kHz

wf = af.waveform.align(wf, mode="peak", interp=True)

Alignment on the negative peak with parabolic interpolation across the three samples around it gives a sub-sample estimate, which matters more than it sounds: at 30 kHz one sample is 33 µs, and a one-sample jitter spreads a cluster in feature space about as much as a genuine amplitude difference of ten per cent.

Clustering#

feat  = af.sort.features(wf, n_components=3)
print(feat.explained_variance_ratio_.round(3))   # [0.612 0.183 0.074]

units = af.sort.cluster(feat, max_k=8, restarts=20, criterion="bic")
print(units.n_units)                             # 7        (24 s)

Three components carry 86.9 % of the variance here, which is typical for a single tetrode group and low for a dense probe. If the first component alone carries more than about 80 %, you are usually looking at amplitude and nothing else, and clustering will simply sort by distance from the electrode.

Curation by numbers#

axonforge does not decide which clusters are units. It computes the four metrics people argue about and leaves the decision in your script, where a reviewer can see it.

Cluster u07 fails on every axis at once: one spike in twelve falls inside the refractory period, the L-ratio is five times the next worst, and it is present in only 62 % of one-minute bins. That combination is the signature of a cluster assembled out of threshold crossings from several sources, not of a neuron.

units = units.drop("u07")
units.save("session12.afz")                   # 18 MB, spike times + waveforms + parameters

The raster#

A two-second window around a stimulus, drawn from the curated table:

Three units increase their rate inside the stimulus window, one of them (u03) with a visible latency of about 50 ms. u06 does the opposite and is suppressed. u04 fires in bursts of four to six spikes that ignore the stimulus entirely — and its 0.72 % interval violations are exactly what a genuine burster produces, which is why the refractory metric is read alongside the others and not alone.

Total time#

• Reading and referencing: 11 s
• Filtering: 19 s
• Detection: 38 s
• Waveform extraction and alignment: 14 s
• Features and clustering: 24 s
• Total: about 1 min 46 s for a ten-minute, 32-channel session

Roughly 60 % of that is spent reading the file, which is why chunked reading is the default and why a faster disk helps more than a faster processor at this channel count.