SpikeInterface / spikeinterface

A Python-based module for creating flexible and robust spike sorting pipelines.
https://spikeinterface.readthedocs.io
MIT License
527 stars 187 forks source link

Incorrect multi-comparison spike trains #1725

Open chimkasinma opened 1 year ago

chimkasinma commented 1 year ago

I am attempting to sort multiple recordings using kilosort2, ironclust, and mountainsort4. However, the agreement spike train (intersection) has spikes even at times when there is no spike from mountainsort4.

image

code:

sorting_comparison = sc.compare_multiple_sorters(
    sorting_list=[KS2_SORTING, MS4_SORTING, IC_SORTING],
    name_list=['KS2', 'MS4', 'IC'],
    spiketrain_mode='intersection',
    verbose=True
)

MIN_AGREE = 3
agreement_units = sorting_comparison.get_agreement_sorting(minimum_agreement_count=MIN_AGREE)
print('Units in agreement for at least {0} sorters: '.format(MIN_AGREE), agreement_units.get_unit_ids())

Thanks for your help!

alejoe91 commented 1 year ago

Hi @chimkasinma

The default mode for spike trains is the union! You can change it with spiketrain_mode='intersection' when calling the compare_* function

chimkasinma commented 1 year ago

Hi @alejoe91! Thank you for your quick response.

I actually did use spiketrain_mode='intersection' in the code, as seen above. I looked at the other plots that were ran and I figured out that their agreements were not unions either. image

ghost commented 1 year ago

I also have a similar issue. I'm comparing the outputs of three spike sorters and want to get the spike trains for any matching units, however, the AgreementSortingExtractor object seems to have incorrect spike trains.

Here's some example code:

pykilo_sort = si.load_extractor('pykilosort/in_container_sorting')
klusta_sort = si.load_extractor('klusta/in_container_sorting')
ironclust_sort = si.load_extractor('ironclust/in_container_sorting')

comp = sc.compare_multiple_sorters(
    sorting_list = [pykilo_sort,klusta_sort,ironclust_sort],
    name_list = ['pykilosort','klusta','ironsort'],
    spiketrain_mode='intersection',
    delta_time=1
)

agr = comp.get_agreement_sorting(minimum_agreement_count=3)

In my case I only have one unit in common between all three sorters. Plotting a section of the sorter spike trains alongside the AgreementSortingExtractor's spike train gives the following:

spike_comparison

I've tried to change the delta_time parameter but this doesn't seem to change the output. I believe it should have far more spikes in the intersection based on my plot. Any ideas why this could be? Thanks in advance for any help!

alejoe91 commented 1 year ago

Hi guys!

I think this was recently fixed by @cristofer-holobetz here https://github.com/SpikeInterface/spikeinterface/pull/1560

Can you try to install from the current main branch? We will make a release in the first week of July ;)

ghost commented 1 year ago

spike_comparison_2

Brilliant! Thanks for the quick respose @alejoe91, and thanks for the fix @cristofer-holobetz!

cristofer-holobetz commented 1 year ago

Hi @alejoe91!

I'm actually working with @chimkasinma, so her version of spikeinterface includes the fix I made. The problem she's experiencing is something else. Thank you for any insight you might have about the issue!

I also wonder if this bug only occurs at the level of creating intersection spiketrains or if it may also be changing the number of shared units identified across sorters. Of course we can't answer this question until we know the source of the bug, but it's something I'd be interested to know.

@rbedfordwork I'm glad it fixed your problem!

alejoe91 commented 1 year ago

Ahhh now I recall! The spike times are taken by the intersection between the two sorters with the highest agreement on a unit-by-unit basis: see https://github.com/SpikeInterface/spikeinterface/blob/main/src/spikeinterface/comparison/multicomparisons.py#L116-L135

cristofer-holobetz commented 1 year ago

Ah thanks for finding that!

Would you all be willing to implement a feature allowing one to choose how conservative an intersection users would like?

Just as one can choose the minimum_agreement_count when finding units in an AgreementSorting, it might be nice to be able to specify that we want the intersection of spikes from the k highest agreement sorters for a particular unit. If that level of granularity is non-trivial to implement, could there be just two modes? One that follows the current behavior and another to extract the "true" intersection across all sorters?

Thanks :)

alejoe91 commented 1 year ago

@cristofer-holobetz yes that sounds like a useful addition. Wanna give it a try?