Other_analysis_methods

Tutorial

The goal of this tutorial is to introduce basic aspects of working with neuropixels data and to point to some useful existing code, to help you get started quickly.

This brief tutorial specifically covers Neuropixels Phase3A data in Matlab, and in some cases specifically for files recorded with SpikeGLX and/or processed with Kilosort. But it will hopefully be useful also as a guide for other situations.

The code here is also contained in the exampleScript.m file of the spikes repository. You will also need the npy-matlab repository to load kilosort/phy files.

Loading kilosort/phy data easily

Kilosort output files can be loaded directly:

spikeTimes = readNPY('spike_times.npy');

But more convenient will be to load many of the relevant files quickly with one command:

myKsDir = 'C:\...\data\myKilosortOutputDirectory';
sp = loadKSdir(myKsDir)

sp.st are spike times in seconds, and sp.clu are cluster identities.

Spikes from clusters labeled "noise" have already been omitted.

Analyzing Drift

Because the Neuropixels probes sample densely of a long continuous span, movement of the brain relative to the electrode can be well detected and quantified. Methods to correct for it are under development.

To observe whether there was drift in your recording, a useful type of plot is the “driftmap”. Make sure to zoom in on the y-axis to look carefully. A drift of 20µm can move a neuron off your site, but the probe is 4mm long.

[spikeTimes, spikeAmps, spikeDepths, spikeSites] = ksDriftmap(myKsDir);
figure; plotDriftmap(spikeTimes, spikeAmps, spikeDepths);

Note that the channel map file used here has channel 1 at y-position=0, and since channel 1 is the site nearest the tip of the probe, this plot goes from the tip of the probe at the bottom to the most superficial part at the top.

Quantification of spiking amplitudes

To see where on the probe spikes of different amplitudes were recorded, we can plot a colormap of the distribution of spikes across depth and amplitude.

depthBins = 0:40:3840;
ampBins = 0:30:min(max(spikeAmps),800);
recordingDur = sp.st(end);

[pdfs, cdfs] = computeWFampsOverDepth(spikeAmps, spikeDepths, ampBins, depthBins, recordingDur);
plotWFampCDFs(pdfs, cdfs, ampBins, depthBins);

Basic LFP characterization

Here we compute some power spectra from the recorded local field potentials and plot them across depth. The method that computes them ("lfpBandPower") can be a useful starting point for seeing how to load LFP data.

lfpD = dir(fullfile(myKsDir, '*.lf.bin')); % LFP file from spikeGLX specifically
lfpFilename = fullfile(myKsDir, lfpD(1).name);

lfpFs = 2500;  % neuropixels phase3a
nChansInFile = 385;  % neuropixels phase3a, from spikeGLX

[lfpByChannel, allPowerEst, F, allPowerVar] = ...
    lfpBandPower(lfpFilename, lfpFs, nChansInFile, []);

chanMap = readNPY(fullfile(myKsDir, 'channel_map.npy'));
nC = length(chanMap);

allPowerEst = allPowerEst(:,chanMap+1)'; % now nChans x nFreq

% plot LFP power
dispRange = [0 100]; % Hz
marginalChans = [10:50:nC];
freqBands = {[1.5 4], [4 10], [10 30], [30 80], [80 200]};

plotLFPpower(F, allPowerEst, dispRange, marginalChans, freqBands);

Computing some useful properties of the spikes and templates

One function is particularly useful for some standard computations you might want to do, like the position on the probe of each template, and the waveform duration of each template.

[spikeAmps, spikeDepths, templateYpos, tempAmps, tempsUnW, tempDur, tempPeakWF] = ...
    templatePositionsAmplitudes(sp.temps, sp.winv, sp.ycoords, sp.spikeTemplates, sp.tempScalingAmps);

Loading synchronization data

Neuropixels Phase3a probes have 16 digital inputs recorded alongside the data from the probe. In spikeGLX these are encoded as the 16 bits in the last (385th) row of the data file. This is true of both the LFP and AP band files, but it's much quicker to load the information from the LFP file (because it is smaller) and the temporal resolution should be sufficient for most applications (2.5kHz). So to load and parse them:

syncChanIndex = 385;
syncDat = extractSyncChannel(myKsDir, nChansInFile, syncChanIndex);

eventTimes = spikeGLXdigitalParse(syncDat, lfpFs);

% - eventTimes{1} contains the sync events from digital channel 1, as three cells: 
% - eventTimes{1}{1} is the times of all events
% - eventTimes{1}{2} is the times the digital bit went from off to on
% - eventTimes{1}{2} is the times the digital bit went from on to off

To synchronize the event times recorded here with those recorded on another system, like a NI board or the FPGA of another probe, it's important to not just find an offset but also a slope. I find that the clocks on the probes drift by up to 50msec/hour relative to each other. A linear regression between the times in the two systems will do this, and to make it simple you can use these functions:

[~,b] = makeCorrection(syncTimesProbe1, syncTimesProbe2, false);
correctedSpikeTimes = applyCorrection(spikeTimesProbe2, b);

Looking at PSTHs aligned to some event

Loading raw waveforms (or spike-triggered LFP, e.g.)