Shuffled channels when concatenating sliced recordings

[hornauerp]

I tried to concatenate recording slices (subselection of channels) and realized that my sorting (Kilosort2.5) did not find consistent units across the recordings:

My workflow is as follows: I have several recordings with some core electrodes that stay constant throughout all recordings, while I sample from the rest of the electrodes "at random". I match the core electrodes and concatenate the recordings over the matching channels, which I then spike sort to subsequently extract footprints from the randomly selected channels.

I went back and checked how the traces looked like after concatenating and filtering. Interestingly, the channels are only shuffled across recordings, when I plot (si.plot_timeseries) or retrieve (si.get_traces) several traces at a time, but not when I look at individual ones:

vs

I tried to backtrace the error but could not figure out where the difference was between retrieving one vs many traces. One thing that I found, however, was that the shuffling did not take place, when the original channel ids (before concatenating) matched.

I realize that this is a pretty niche issue, but I would appreciate any help, as I would like to keep my workflow within the spikeinterface environment.

Thank you very much in advance!

[alejoe91]

Hi Philipp!

Thanks for all the info, this is helpful. I guess something is buggy when re-sorting the channels..can you run this to double check?

# rec_concatenated  is your concatenated recording

trace_one_channel_only = rec_concatenated.get_traces(channel_ids=rec_concatenated.channel_ids[0])[:, 0]
trace_multi_channels = rec_concatenated.get_traces(channel_ids=rec_concatenated.channel_ids[:2])[:, 0]

np.testing.assert_array_equal(trace_one_channel_only, trace_multi_channels)

Can I ask you some clarification on?

One thing that I found, however, was that the shuffling did not take place, when the original channel ids (before concatenating) matched.

[hornauerp]

Hi Alessio!

When I try to run your code I get the following error:

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[13], line 1
----> 1 trace_one_channel_only = multirecording.get_traces(channel_ids=multirecording.channel_ids[0])[:, 0]
      2 trace_multi_channels = multirecording.get_traces(channel_ids=multirecording.channel_ids[:2])[:, 0]

File ~/miniconda3/envs/si_env/lib/python3.9/site-packages/spikeinterface/core/baserecording.py:158, in BaseRecording.get_traces(self, segment_index, start_frame, end_frame, channel_ids, order, return_scaled, cast_unsigned)
    126 """Returns traces from recording.
    127 
    128 Parameters
   (...)
    155     If return_scaled is True, but recording does not have scaled traces
    156 """
    157 segment_index = self._check_segment_index(segment_index)
--> 158 channel_indices = self.ids_to_indices(channel_ids, prefer_slice=True)
    159 rs = self._recording_segments[segment_index]
    160 traces = rs.get_traces(start_frame=start_frame, end_frame=end_frame, channel_indices=channel_indices)

File ~/miniconda3/envs/si_env/lib/python3.9/site-packages/spikeinterface/core/base.py:102, in BaseExtractor.ids_to_indices(self, ids, prefer_slice)
    100 else:
    101     _main_ids = self._main_ids.tolist()
--> 102     indices = np.array([_main_ids.index(id) for id in ids], dtype=int)
    103     if prefer_slice:
    104         if np.all(np.diff(indices) == 1):

TypeError: 'numpy.int64' object is not iterable

I I replace channel_ids[0] by channel_ids[0:1] it works and the assertion does not fail.

Regarding this:

One thing that I found, however, was that the shuffling did not take place, when the original channel ids (before concatenating) matched.

I checked again and this seems to have been a concidence, but instead I now found that also retrieving 2 traces does not shuffle them, but only if it is >=3:

[alejoe91]

Ok, getting closer!

So basically the same example above but with the secon line being:

trace_multi_channels = rec_concatenated.get_traces(channel_ids=rec_concatenated.channel_ids[:3])[:, 0]

fails? That is interesting!

[hornauerp]

Yes exactly, this gives this shuffled trace:

[alejoe91]

ok..what are the channel_ids of the individual recordings?

[hornauerp]

Those should be:

[alejoe91]

can you try to plot the traces for the single recordings separately for the same channels? Thanks for helping out here! Trying to understand what's going on

[hornauerp]

So what I tried before was:

traces_list = [multirecording.recording_list[x].get_traces() for x in range(3)]
conc_traces = np.concatenate(traces_list)

If I plot those it is still shuffled:

Would that answer your question?

Ah, sorry I misread your question, one second.

So I guess you are referring to something like this?

traces_list = [multirecording.recording_list[x].get_traces(channel_ids=[0]) for x in range(3)]
conc_traces = np.concatenate(traces_list)

This seems to look fine:

[alejoe91]

I would try this 2 options:

1. get all traces and slice later

traces_list = [multirecording.recording_list[x].get_traces()[:, 0] for x in range(3)]

2. get one trace at a time:

traces_list = [multirecording.recording_list[x].get_traces(channel_ids=[multirecording.recording_list[x].channel_ids[0]]) for x in range(3)]

Could you also share the code you used for selecting and sorting the channels and to concatenate them?

[hornauerp]

Your first option produces:

The second one:

To select and concatenate the channels I do the following:

rec_list = []
for r in range(rec_id): 
    rec_name = 'rec' + '%0*d' % (4, r)
    rec = si.MaxwellRecordingExtractor(full_path,stream_name=stream_name,rec_name=rec_name)
    rec_el = rec.get_property("contact_vector")["electrode"]
    
    chan_idx = [np.where(rec_el == el)[0][0] for el in common_el[2:6]] #Just for debugging
    sel_channels = rec.get_channel_ids()[chan_idx]
    rec_filtered = si.filter(rec, band=[150,4999])
    rec_list.append(rec_filtered.channel_slice(sel_channels, renamed_channel_ids=list(range(len(chan_idx)))))

multirecording = si.concatenate_recordings(rec_list)

I just do the filtering here to get rid of the offset and it is actually faster than si.center().

[alejoe91]

And the second option seems to work well for all 4 channels, right?

[hornauerp]

Yes, it works for all 4 individually.

[alejoe91]

Thanks. I think I need to ask you to share the data with me! Otherwise we can try to have a debugging by the end of next week. Let me know what you prefer!