Amplitude Cutoff not working

[karmenrai2]

I am trying to calculate the amplitude cutoff for a set of data. Whenever I try to compute it I get an error.

import spikeinterface as si
import spikeinterface.extractors as se
import spikeinterface.widgets as sw
import spikeinterface.qualitymetrics as sqm
from spikeinterface.extractors import read_phy
from probeinterface import get_probe
from pathlib import Path
import numpy as np
import csv

# Cluster iterator
i = 0

# File paths 
rec_file_path = Path(r"C:\Users\cognp\Documents\Kilosort4\kilosort4data_day4\continuous1.dat") # Recording file
sort_file_path = Path(r"C:\Users\cognp\Documents\Kilosort4\kilosort4data_day4") # Folder with the kilosort data
cluster_label = Path(r"C:\Users\cognp\Documents\Kilosort4\kilosort4data_day4\spike_clusters.npy") # spike_clusters npy file
pcs = Path(r"C:\Users\cognp\Documents\Kilosort4\kilosort4data_day4\spike_positions.npy") # spike_positions npy file

# Define recording parameters
sampling_frequency = 30_000.0  
num_channels = 384  
dtype = "float64"  

# Load data using SpikeInterface
recording = si.read_binary(file_paths=rec_file_path, sampling_frequency=sampling_frequency,
                           num_channels=num_channels, dtype=dtype)
sorting = read_phy(sort_file_path)

# Load data using numpy
all_labels = np.load(cluster_label)
all_pcs = np.load(pcs)

# Set the probe
other_probe = get_probe('cambridgeneurotech', 'ASSY-37-E-1')
other_probe.set_device_channel_indices(np.arange(32))
recording = recording.set_probe(other_probe, group_mode='by_shank')

# Create sorting analyzer
analyzer = si.create_sorting_analyzer(sorting=sorting, recording=recording, format="memory")
analyzer.compute(['random_spikes', 'templates'])
amplitudes = analyzer.compute(input="spike_amplitudes", peak_sign="neg")

# Quality Metrics
presence_ratio = sqm.compute_presence_ratios(sorting_analyzer=analyzer)
amplitude_cutoff = sqm.compute_amplitude_cutoffs(sorting_analyzer=analyzer, peak_sign="neg")
isi_violations_ratio, isi_violations_count = sqm.compute_isi_violations(sorting_analyzer=analyzer, isi_threshold_ms=1.0)

# CSV file 
with open('spikesorting.csv', 'w') as file:
    writer = csv.writer(file)
    writer.writerow(["Cluster", "Presence Ratio", "Amplitude Cutoff", "Isolation Distance", "ISI Violations ratio", "ISI Violations count"])
    for cluster in presence_ratio.keys(): 
        isolation_distance, _ = sqm.mahalanobis_metrics(all_pcs=all_pcs, all_labels=all_labels, this_unit_id=cluster)
        i+=1
        print(str(i) + "th cluster done")
        writer.writerow([cluster, presence_ratio[cluster], amplitude_cutoff[cluster], isolation_distance, isi_violations_ratio[cluster], isi_violations_count[cluster]])

# Print data for verification
print(recording)
print(sorting)
print(analyzer)
print(presence_ratio)
print(type(amplitudes))
print(len(amplitudes))
print(amplitudes)
exit()

The error I'm getting is:
`Traceback (most recent call last):
  File "c:\Users\cognp\Documents\npyfiles\sort.py", line 46, in <module>
    amplitude_cutoff = sqm.compute_amplitude_cutoffs(sorting_analyzer=analyzer, peak_sign="neg")
                       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\cognp\AppData\Local\Packages\PythonSoftwareFoundation.Python.3.12_qbz5n2kfra8p0\LocalCache\local-packages\Python312\site-packages\spikeinterface\qualitymetrics\misc_metrics.py", line 896, in compute_amplitude_cutoffs
    amplitudes_by_units = _get_amplitudes_by_units(sorting_analyzer, unit_ids, peak_sign)
                          ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\cognp\AppData\Local\Packages\PythonSoftwareFoundation.Python.3.12_qbz5n2kfra8p0\LocalCache\local-packages\Python312\site-packages\spikeinterface\qualitymetrics\misc_metrics.py", line 814, in _get_amplitudes_by_units
    amplitudes_by_units[unit_id] = all_amplitudes[spike_mask]
                                   ~~~~~~~~~~~~~~^^^^^^^^^^^^
IndexError: boolean index did not match indexed array along dimension 0; dimension is 553228 but corresponding boolean dimension is 2363657

I'm assuming there is some type of discrepancy between the size of the computed spike amplitudes and the way the amplitude cutoff is calculated.

[zm711]

Kilosort often has a problem of generating spikes after the recording (maybe due to padding?). We suggest running remove_excess_spikes so

import spikeinterface.curation as scur
clean_sorting = scur.remove_excess_spikes(sorting, recording)

Then analyze your data with this :)

[karmenrai2]

I've added these lines into my code like this:

import spikeinterface.extractors as se
import spikeinterface.widgets as sw
import spikeinterface.qualitymetrics as sqm
import spikeinterface.curation as scur
from spikeinterface.extractors import read_phy
from probeinterface import get_probe
from pathlib import Path
import numpy as np
import csv

# Cluster iterator
i = 0

# File paths 
rec_file_path = Path(r"C:\Users\cognp\Documents\Kilosort4\kilosort4data_day4\continuous1.dat") # Recording file
sort_file_path = Path(r"C:\Users\cognp\Documents\Kilosort4\kilosort4data_day4") # Folder with the kilosort data
cluster_label = Path(r"C:\Users\cognp\Documents\Kilosort4\kilosort4data_day4\spike_clusters.npy") # spike_clusters npy file
pcs = Path(r"C:\Users\cognp\Documents\Kilosort4\kilosort4data_day4\spike_positions.npy") # spike_positions npy file

# Define recording parameters
sampling_frequency = 30_000.0  
num_channels = 384  
dtype = "float64"  

# Load data using SpikeInterface
recording = si.read_binary(file_paths=rec_file_path, sampling_frequency=sampling_frequency,
                           num_channels=num_channels, dtype=dtype)
sorting = read_phy(sort_file_path)

# Load data using numpy
all_labels = np.load(cluster_label)
all_pcs = np.load(pcs)

# Set the probe
other_probe = get_probe('cambridgeneurotech', 'ASSY-37-E-1')
other_probe.set_device_channel_indices(np.arange(32))
recording = recording.set_probe(other_probe, group_mode='by_shank')

# Remove the excess spikes from the kilosort data 
clean_sorting = scur.remove_excess_spikes(sorting, recording)

# Create sorting analyzer
analyzer = si.create_sorting_analyzer(sorting=sorting, recording=recording, format="memory")
analyzer.compute(['random_spikes', 'templates'])
amplitudes = analyzer.compute(input="spike_amplitudes", peak_sign="neg")

# Quality Metrics
presence_ratio = sqm.compute_presence_ratios(sorting_analyzer=analyzer)
amplitude_cutoff = sqm.compute_amplitude_cutoffs(sorting_analyzer=analyzer, peak_sign="neg")
isi_violations_ratio, isi_violations_count = sqm.compute_isi_violations(sorting_analyzer=analyzer, isi_threshold_ms=1.0)

# CSV file 
with open('spikesorting.csv', 'w') as file:
    writer = csv.writer(file)
    writer.writerow(["Cluster", "Presence Ratio", "Amplitude Cutoff", "Isolation Distance", "ISI Violations ratio", "ISI Violations count"])
    for cluster in presence_ratio.keys(): 
        isolation_distance, _ = sqm.mahalanobis_metrics(all_pcs=all_pcs, all_labels=all_labels, this_unit_id=cluster)
        i+=1
        print(str(i) + "th cluster done")
        writer.writerow([cluster, presence_ratio[cluster], amplitude_cutoff[cluster], isolation_distance, isi_violations_ratio[cluster], isi_violations_count[cluster]])

# Print data for verification
print(recording)
print(sorting)
print(analyzer)
print(presence_ratio)
print(type(amplitudes))
print(len(amplitudes))
print(amplitudes)
exit()```

however I'm still getting the same error. 

[zm711]

analyzer = si.create_sorting_analyzer(sorting=sorting, recording=recording, format="memory")

This is wrong. You need to put in your clean_sorting :)

[karmenrai2]

After using this, much of the data I get is NaN. At first it wouldn't let me compute the amplitudes at all and gave me this error:

    amplitude_cutoff = sqm.compute_amplitude_cutoffs(sorting_analyzer=analyzer, peak_sign="neg")
                       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\cognp\AppData\Local\Packages\PythonSoftwareFoundation.Python.3.12_qbz5n2kfra8p0\LocalCache\local-packages\Python312\site-packages\spikeinterface\qualitymetrics\misc_metrics.py", line 903, in compute_amplitude_cutoffs
    all_fraction_missing[unit_id] = amplitude_cutoff(
                                    ^^^^^^^^^^^^^^^^^
  File "C:\Users\cognp\AppData\Local\Packages\PythonSoftwareFoundation.Python.3.12_qbz5n2kfra8p0\LocalCache\local-packages\Python312\site-packages\spikeinterface\qualitymetrics\misc_metrics.py", line 1266, in amplitude_cutoff
    h, b = np.histogram(amplitudes, num_histogram_bins, density=True)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\cognp\AppData\Local\Packages\PythonSoftwareFoundation.Python.3.12_qbz5n2kfra8p0\LocalCache\local-packages\Python312\site-packages\numpy\lib\histograms.py", line 780, in histogram
    bin_edges, uniform_bins = _get_bin_edges(a, bins, range, weights)
                              ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\cognp\AppData\Local\Packages\PythonSoftwareFoundation.Python.3.12_qbz5n2kfra8p0\LocalCache\local-packages\Python312\site-packages\numpy\lib\histograms.py", line 426, in _get_bin_edges
    first_edge, last_edge = _get_outer_edges(a, range)
                            ^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\cognp\AppData\Local\Packages\PythonSoftwareFoundation.Python.3.12_qbz5n2kfra8p0\LocalCache\local-packages\Python312\site-packages\numpy\lib\histograms.py", line 323, in _get_outer_edges
    raise ValueError(
ValueError: autodetected range of [nan, nan] is not finite 

In order to get some output I changed line 1266 in "misc_metrics.py" to this:

h, b = np.histogram(amplitudes, num_histogram_bins, (np.nanmin(amplitudes), np.nanmax(amplitudes)), density=True) 

Now my output for amplitude cutoff looks something like this:

Amplitude Cutoff

 
nan
 
0.227006181187956
 
nan
 
nan
 
nan
 
nan
 
nan
 
nan
 
nan
 
nan
 
nan
 
nan
 
nan
 
nan
 
nan
 
0.2318217716589
 
nan
 
nan
 
nan
 
nan
 
nan
 
nan
 
nan

[zm711]

@karmenrai2,

amplitude cutoff is a pretty strict metric (I often find that my neurons don't fit with the requirements). How many spikes are you getting/neuron. Looking at your script more I'm still seeing a lot mixed analysis inside of spikeinterface and outside of spikeinterface (for example I don't think your "pcs" that you are loading make sense.

I would suggest cleaning up the script a bit so we can make sure things are working. I would suggest starting out by looking at something like snr. That is an easy quality metric that doesn't have any assumptions. It just looks at your spikes and the noise levels. More "stat" metrics require your data to fit assumptions and without seeing your data itself it's hard to know why it isn't working.

[karmenrai2]

I was trying now to get the raw amplitudes sorted by unit. In the documentation it stated:

outputs“numpy” | “by_unit”, default: “numpy” The output format, either concatenated as numpy array or separated on a per unit basis

I tried doing this to get the raw amplitudes separated per unit:

amplitudes = si.compute_spike_amplitudes(sorting_analyzer=analyzer, outputs="by_unit")

However I'm getting this error:

Traceback (most recent call last): File "c:\Users\cognp\Documents\npyfiles\sort.py", line 47, in <module> amplitudes = si.compute_spike_amplitudes(sorting_analyzer=analyzer, outputs="by_unit") ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "C:\Users\cognp\AppData\Local\Packages\PythonSoftwareFoundation.Python.3.12_qbz5n2kfra8p0\LocalCache\local-packages\Python312\site-packages\spikeinterface\core\sortinganalyzer.py", line 1971, in __call__ ext = sorting_analyzer.compute(cls.extension_name, *args, **kwargs) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "C:\Users\cognp\AppData\Local\Packages\PythonSoftwareFoundation.Python.3.12_qbz5n2kfra8p0\LocalCache\local-packages\Python312\site-packages\spikeinterface\core\sortinganalyzer.py", line 1324, in compute return self.compute_one_extension(extension_name=input, save=save, verbose=verbose, **kwargs) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ File "C:\Users\cognp\AppData\Local\Packages\PythonSoftwareFoundation.Python.3.12_qbz5n2kfra8p0\LocalCache\local-packages\Python312\site-packages\spikeinterface\core\sortinganalyzer.py", line 1401, in compute_one_extension extension_instance.set_params(save=save, **params) File "C:\Users\cognp\AppData\Local\Packages\PythonSoftwareFoundation.Python.3.12_qbz5n2kfra8p0\LocalCache\local-packages\Python312\site-packages\spikeinterface\core\sortinganalyzer.py", line 2331, in set_params params = self._set_params(**params) ^^^^^^^^^^^^^^^^^^^^^^^^^^ TypeError: ComputeSpikeAmplitudes._set_params() got an unexpected keyword argument 'outputs'