Issue with Shape Mismatch Error in calculate_drift_metrics function in metrics.py

[OlyOlympic]

Hello, @jsiegle

I am encountering an issue when using the calculate_drift_metrics function in the ecephys_spike_sorting\modules\quality_metrics\metrics.py file. The error occurs due to a shape mismatch between variables, and the shapes always differ by 1. Below is the relevant function code:
def calculate_drift_metrics(spike_times,
spike_clusters,
spike_templates,
unit_template_ids,
total_units,
pc_features,
pc_feature_ind,
channel_pos,
interval_length,
min_spikes_per_interval)

I am receiving a ValueError related to mismatched shapes when I try to run this function. The dimensions of the variables seem to differ by 1 each time, which causes the error.

_File d:\bo\ecephys\ecephys_spike_sorting-master\ecephys_spike_sorting\modules\quality_metrics\metrics.py:142, in calculate_metrics(spike_times, spike_clusters, spike_templates, amplitudes, channel_map, channel_pos, templates, pc_features, pc_feature_ind, params, epochs)
133 the_silhouette_score = calculate_silhouette_score(spike_clusters[in_epoch],
134 spike_templates[in_epoch],
135 total_units,
136 pc_features[in_epoch,:,:],
137 pc_feature_ind,
...
556 # majority template for htat cluster. These operations make copies of the
557 # arrays.
558 m_spike_clusters = spike_clusters[match_maj]

ValueError: operands could not be broadcast together with shapes (2382853,) (2382852,)_

Could you help me understand why this happens and how I can resolve the issue?
I would really appreciate any insights or suggestions for solving this problem!
Thanks in advance!

Best regards,
Oly

[jsiegle]

Looping in @jenniferColonell -- I believe you are using her forked version of this repo

[jenniferColonell]

Hi Oly -- Sorry you're having this trouble. Do you see it on all runs, or just some? What other processing of the phy data has happened (e.g. have you gone through removing duplicates?).

[jenniferColonell]

Hi again -- This issue stems from on old bug in the metric calculation that would exclude spikes with time=0, and then not properly re-dimension all of the arrays. If you grab the latest version, it will be fixed. Note that the easiest way to do that is generally to get a new clone and build a new virtual environment. create_input_json and the driver scripts have some changes, so edit the copies in the new clone rather than copying over old versions.
If you have been modifying your copy of the pipeline and prefer not to update, I can tell you what to edit to fix just this issue.

[OlyOlympic]

@jenniferColonell
Hi,

Thanks for the explanation. Yes, this issue occurs on all runs, not just some. Since my current environment involves multiple processing steps and this is just one part of the pipeline, I'd prefer not to do a full update if it can be avoided.

If it’s not too much trouble, could you please show me how to manually fix this specific issue in my current setup? That would be really helpful.

Best,
Oly

[jenniferColonell]

Well, if you are seeing this on every run, I'm not sure this is the correct answer. Are you using defined epochs in your analysis?
If you are not using epochs, these errors should only occur if there are spikes at time = 0.

The changes to the code are in metrics.py. Keep in mind that these edits will only take effect if you have an editable version of the pipeline. If you've moved the code to your own pipeline, of course that works too.

In metrics.py, line 88, make sure the that the spike in the epoch are defined with >= (rather than just >)

    in_epoch = (spike_times >= epoch.start_time) * (spike_times <= epoch.end_time)

In the call to calculate drift metrics, make sure spike_templates is just the spikes in the current epoch:

max_drift, cumulative_drift = calculate_drift_metrics(spike_times[in_epoch],
spike_clusters[in_epoch],
spike_templates[in_epoch],
template_ids,
total_units,
pc_features[in_epoch,:,:],
pc_feature_ind,
channel_pos,
params['drift_metrics_interval_s'],
params['drift_metrics_min_spikes_per_interval'])

[OlyOlympic]

Hi, @jenniferColonell

Thanks for your reply. I tested several different sessions and all of them showed the same error, though I’m not completely sure whether ‘every’ session would produce it.

I was using calculate_metrics with epochs=None, so I wasn’t using defined epochs. I modified line 88 in metrics.py as you suggested, but after rerunning the code, I still encountered the same error:

File "d:\ecephys\ecephys_spike_sorting-master\ecephys_spike_sorting\modules\quality_metrics\metrics.py", line 541, in calculate_drift_metrics
maj_tid = unit_template_ids[spike_clusters]
match_maj = spike_templates == maj_tid
ValueError: operands could not be broadcast together with shapes (2382853,) (2382852,)

Below is the code I’m using to call metrics:

import numpy as np
import os
from ecephys_spike_sorting.modules.quality_metrics import metrics
import pandas as pd
from collections import OrderedDict

folder = r'D:\EphysData\real\sound2\kilosort4'
print(folder)
print("loading data")

def load_data(folder):
spike_times = np.load(os.path.join(folder,'spike_times.npy'))
spike_clusters = np.load(os.path.join(folder,'spike_clusters.npy'))
spike_templates = np.load(os.path.join(folder, 'spike_templates.npy'))
amplitudes = np.load(os.path.join(folder,'amplitudes.npy'))
channel_map = np.load(os.path.join(folder, 'channel_map.npy'))
channel_pos = np.load(os.path.join(folder, 'channel_positions.npy'))
templates = np.load(os.path.join(folder, 'templates.npy'))
pc_features = np.load(os.path.join(folder,'pc_features.npy'))
pc_feature_ind = np.load(os.path.join(folder,'pc_feature_ind.npy'))

return spike_times, spike_clusters, spike_templates, amplitudes, channel_map, channel_pos, templates, pc_features, pc_feature_ind
spike_times, spike_clusters, spike_templates, amplitudes, channel_map, channel_pos, templates, pc_features, pc_feature_ind = load_data(folder)

params = {}
params['include_pcs'] = True
params['isi_threshold'] = 0.0015
params['min_isi'] = 0.000166
params['tbin_sec'] = 0.001
params['max_radius_um'] = 250
params['max_spikes_for_unit'] = 500
params['max_spikes_for_nn'] = 10000
params['n_neighbors'] = 4
params['n_silhouette'] = 10000
params['drift_metrics_interval_s'] = 51
params['drift_metrics_min_spikes_per_interval'] = 10

metrics.calculate_metrics(
spike_times=spike_times,
spike_clusters=spike_clusters,
spike_templates=spike_templates,
amplitudes=amplitudes,
channel_map=channel_map,
channel_pos=channel_pos,
templates=templates,
pc_features=pc_features,
pc_feature_ind=pc_feature_ind,
params=params,
epochs=None
)

Do you see any issues with how I’m using metrics? Or is there something incorrect in my code that might explain why the fix to line 88 didn’t resolve the error?

I really appreciate your help, thank you so much!

Best,
Oly

[jenniferColonell]

Hi Oly,

It's not causing your problem, but the drift metrics do assume that the spike times are in seconds. You can see how the pipeline loads data here:

https://github.com/jenniferColonell/ecephys_spike_sorting/blob/master/ecephys_spike_sorting/common/utils.py#L246

Obviously I've run this code quite a bit, and it usually works. So, there's something unique about your input data. Specifically, it looks like there is an extra entry in spike_templates. Can you try loading your spike_clusters.npy and spike_templates.npy files and see if they match? If the input arrays match, you can try just running the calculation that gets the majority template for each cluster, i.e.:

total_units = np.max(spike_clusters) + 1
template_ids = np.zeros((total_units,), dtype='uint16')
spike_templates = np.squeeze(spike_templates)
cluster_ids = np.unique(spike_clusters)

for cid in cluster_ids:
cluster_templates = spike_templates[np.where(spike_clusters==cid)]
template_ids[cid] = np.argmax(np.bincount(cluster_templates))

maj_tid = template_ids[spike_clusters]
match_maj = spike_templates==maj_tid

What's happening here is that we're getting the index of the template that was used to extract the majority of the spikes in a cluster -- a cluster can contain spikes extracted with different templates if there has been manual curation in phy. We then make a list of spikes that match the majority template, and that it's those that are used in the drift calculation.

If your data has some unexpected properties, there might be some odd results. What sorter are you running?

The line 88 fix is specific to using defined epochs -- so not related to your problem.

[OlyOlympic]

Hi, @jenniferColonell

I’ve checked and confirmed that the spike_clusters.npy and spike_templates.npy files have the same shape. My data was processed using Kilosort 4.

I’ve resolved the issue by extracting and using the specific functions and tools from the jsiegle (https://github.com/AllenInstitute/ecephys_spike_sorting) directly. Thanks a lot for your help!

Best,
Oly

[jenniferColonell]

HI @OlyOlympic Sorry you've had all this trouble -- I do run this code on KS4 data regularly, so I think there must have been an issue with how you were reading in your arrays. Glad Josh's earlier code works for your purpose, I agree it should be fine if you don't need to account for manual curation.