Kilosort output for sorting analyzer

[Djoels]

As mentioned in the KS (2.5/4) documentation, part of the sorting output prepared for curation is a number of npy files containing templates, similar templates, PC's, ... which are also calculated using a sorting analyzer for a number of metrics.

Could these features somehow be reused, or are they incompatible with sorting analyzer way of working?

Furthermore, what would be the advice when analyzing a recording processed by Kilosort: use the original recording or use the temp_wh.dat file produced by KS which contains already preprocessed (car, whitening, highpass filtered) data?

[Djoels]

Hmm, I already stumbled upon an indication of whether the temp_wh.dat file should be used: I guess not then, because it needs to be "unwhitened" before applying sorting analyzer?

[chrishalcrow]

Hey @Djoels I'm actively working on this! There's already some functionality to import kilosort output as an analyzer (https://spikeinterface.readthedocs.io/en/latest/how_to/import_kilosort_data.html) You'll need use the main branch of spikeinterface (locally clone and install rather than just "pip install spikeinterface"). At the moment it just loads the "easy" stuff.

But I have a branch here which adds the "hard" stuff: main...chrishalcrow:spikeinterface:add-amps-to-kilosort-convert but it's not quite ready for a proper PR yet. You've already figured out the hard stuff is: some output is whitened, some is in PC space, some is in physical space, and we can only approximately go between the spaces. But think I've figured it all out for KS4 output now.

For an analyzer, I take a raw recording and lightly preprocess it (bandpass + CMR) then combine with the kilosort sorting output - this seems to work pretty well for my data (chronic NP2).

[Djoels]

That is great news!
At our labs we're mostly still using KS2.5 so I'm keen to find out whether this would also work.
Our recordings are mostly on custom platform or Maxwell recordings (so in vitro HD-MEA's).
Will give it a try shortly, thanks!

[pas-calc]

What happens if templates were split or merged after kilosort as a manual curation? I think in this case it is crucial to only read the spike_templates.npy and not the spike_clusters.npy file. The BasePhyKilosortSortingExtractor tries to load spike_clusters.npy if present. But here it is important that metrics are only calculated for templates and are not updated after manual curation.
https://phy.readthedocs.io/en/latest/terminology/#cluster-vs-templates
https://kilosort.readthedocs.io/en/latest/export_files.html

Alternatively we can load both files and then decompose a cluster into the template like this:

# load
st = np.load("spike_templates.npy")
sc = np.load("spike_clusters.npy")
clusterids = np.unique(sc)

# decompose cluster into templates
def decompose_cluster(c):
    """
    phy split/merge leads to cluster = sum of templates. Function to decompose a cluster in its templates
    Returns: templates with its weights (proportion of number of spikes)
    """
    if not c in clusterids:
        raise ValueError("invalid cluster: {} from {} clusters".format(c,len(clusterids)))
    s_ind = (sc==c) # get spikes associated to that cluster c (np.where)
	s_templates = st[s_ind] # get templates associated to these spikes
	# now you have: cluster -> spikes -> templates, and thus a mapping from the cluster c to the templates on which the corresponding spikes were detected
	unique, counts = np.unique(s_templates, return_counts=True) # get the distribution of the templates from the cluster c
	weights = counts / np.sum(counts) # normalize
	return unique, weights
            
# you could run this to decompose each cluster in its templates
for c in clusterids:
    print("cluster:",c)
    unique, weights = decompose_cluster(c)
    print(dict(zip(unique, weights)))
    print("")

We could apply these weights on the data (waveforms...) but it won't be accurate. For example if we split a template into two clusters, we don't know the waveform of each. It would assign the same

[chrishalcrow]

Hello, that's beyond the scope of what I'm trying to do. I just want to be able to load in plain, untouched kilosort output (for now!)

[pas-calc]

Hello, that's beyond the scope of what I'm trying to do. I just want to be able to load in plain, untouched kilosort output (for now!)

Yes I can understand, and I think it is better to use the actual output, not the manual refinement.
I just did not see it either here #4202 but just one should make sure to only use the spike_templates.npy file

Just want to warn about sorting = read_phy(phy_path) in read_kilosort_as_analyzer because it guesses the npy file

spikeinterface/src/spikeinterface/extractors/phykilosortextractors.py

Lines 80 to 83 in 24c4c0d

	if (phy_folder / "spike_clusters.npy").is_file():
	spike_clusters = np.load(phy_folder / "spike_clusters.npy")
	else:
	spike_clusters = np.load(phy_folder / "spike_templates.npy")