SpikeInterface / spikeinterface

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

Missing channels in phy #566

Closed mfvd closed 2 years ago

mfvd commented 2 years ago

Hi

When doing the manual curation in phy, the WaveformView is missing a number of channels (pic below). Curiously, the missing channels are not consistent between units and there doesn't seem to be a discernible pattern.

Also, when checking the units in phy, it seems like the channel position might have changed. I suspect this might have happened only because phy displays high amplitude waveforms of a given unit in fairly distant channels (e.g. channels 30 and 22 in the pic).

I used the following code to define my probe to device map:

manufacturer = 'cambridgeneurotech'
probe_name = 'ASSY-37-Fb'
probe = pi.get_probe(manufacturer, probe_name)

# Match to device (RHD headstage) map. In this case it was manually corrected because probeinterface might be wrong.
device_channel_indices = [
    29, 19, 10, 18, 28, 30, 20, 17, 21, 31, 22, 16, 23, 27, 26, 25, 24, 7, 6, 5,
    4, 8, 9, 3, 11, 2, 12, 1, 13, 0, 14, 15]
# set channel ids based on manual probe map from cambridge site and RHD map.
probe.set_device_channel_indices(device_channel_indices)

# set recording channels to new probe>RHD map.
rec_probe = rec.set_probe(probe, 
                          group_mode='by_shank')

and the following code to export to phy:

waveforms_ms4 = si.WaveformExtractor.create(rec_cashed, 
                                            sorted_ms4,
                                            'waveforms', 
                                            remove_if_exists = True)

# set params
waveforms_ms4.set_params(ms_before = 2, 
                         ms_after = 2, 
                         max_spikes_per_unit = 1000,
                         return_scaled=True)

# run  extractor to get waveforms
waveforms_ms4.run_extract_waveforms(n_jobs = -1, 
                                    chunk_size = 30000)

# export spikes to phy
export_to_phy(waveforms_ms4, 
              './phy_ms4',
              compute_pc_features = True, 
              compute_amplitudes = True,
              progress_bar = True)

Any idea of what might have gone wrong? Cheers

waveforms

alejoe91 commented 2 years ago

Hi @mfvd

When are you saving the rec_cashed object? Can you share the entire code?

samuelgarcia commented 2 years ago

Hi. Note that when we export to phy only best X best channels per units are exported. You can control this by max_channels_per_template which is 16 by default. So the sparsity on waveform is something expected. Loking at you display, there is maybe a bug in this sparsity and/or the channel position. I don't know if it is in phy or not.

Mor importantly : you set the channel mapping (wiring) manually because you say that the RDH mapping in probeinterface is wrong. If we are wrong it is super important to make a feedback on this!!!!!

Could you do this with your mapping and the probeinterface one so we could double check this:

dataframe_papping = probe.to_dataframe(complete=True).loc[:, ["contact_ids", "shank_ids", "device_channel_indices"]]
print(dataframe_papping)

If the mapping is wring from the beginning then it could explain that channels position are super strange. This mapping is critical!

mfvd commented 2 years ago

Hi again,

@alejoe91 sorry, I assumed the error might be in the code above. Here's the rest of the code:

rec_f = st.bandpass_filter(rec_probe,
                           freq_min=300,
                           freq_max=6000,
                           ftype='butter')

rec_cmr = st.common_reference(rec_f,
                              reference='global',
                              operator='median')

# save data
rec_cashed = rec_cmr.save(format='binary',
                          folder="preprocessed",
                          n_jobs=8,
                          total_memory="3G")

# get dafault params for ms4
default_ms4_params = ss.get_default_params('mountainsort4')
# set params
default_ms4_params['detect_threshold'] = 4
default_ms4_params['detect_sign'] = -1
default_ms4_params['filter'] = False
default_ms4_params['whiten'] = True

# run sorter
sorted_ms4 = ss.run_mountainsort4(recording = rec_cashed,
                                  output_folder = 'ms4_output', 
                                  verbose = False,
                                  **default_ms4_params)

@samuelgarcia thanks, that explains the missing channels and matches the 16 waveforms I see for each unit.

I reported a possible mismatch in the channel map in https://github.com/SpikeInterface/probeinterface/issues/92

In the meanwhile I did the mapping manually in such a way that the device_channel_indices match the Assy-37-Fb to RHD wiring (they are in the correct locations when plotting the probe with device channels and contact_ids). For this I used the wiring diagrams provided by cambridge neurotech and Intan. See pics below.

probe_map

shank0

shank1_2

samuelgarcia commented 2 years ago

OK sorry I did not make the link between this issue and the one you did in probeinterface sorry. So if I understand correctly, you do this mapping to compensate the error in the probe layout we implemented.

mfvd commented 2 years ago

@samuelgarcia no problem! I should have linked to that issue in this one.

Yes. I get the probe map from probeinterface and do the wiring to the headstage by manually setting the device_channel_indices. That is the only thing I change in the probe object. The assumption is that, both the sorter and phy, only care about the device channel ids.

alejoe91 commented 2 years ago

@mfvd do you get this for all units? One possible explanation could be that the unit is overmerged, i.e., it actually has spikes from more than one neuron!

mfvd commented 2 years ago

@alejoe91 unfortunately I do and in all units I have spikes from contacts across different shanks. I'll try to run the sorter by shank and see what changes.

mfvd commented 2 years ago

When comparing the group/shank_ids in the probe object with the cashed recording I found that they are different.

I'm not sure but it seems that the channels are somehow getting misplaced when assigning the probe information to the recording object. Channel groups do not match the shank_ids in the probe (see pic below).

I'm now wondering if I'm using the correct method to manually redo the wiring (probe to headstage).

Code to assign probe information to my recording

manufacturer = 'cambridgeneurotech'
probe_name = 'ASSY-37-Fb'
probe = pi.get_probe(manufacturer, probe_name)

# Match to device (RHD headstage) map. In this case it was manually corrected because probeinterface might be wrong.
device_channel_indices = [
    29, 19, 10, 18, 28, 30, 20, 17, 21, 31, 22, 16, 23, 27, 26, 25, 24, 7, 6, 5,
    4, 8, 9, 3, 11, 2, 12, 1, 13, 0, 14, 15]
# set channel ids based on probe map from cambridge site and RHD map.
probe.set_device_channel_indices(device_channel_indices)

# set recording channels to new probe>RHD map.
# by shank to perform sorting by shank
rec_probe = rec.set_probe(probe, 
                          group_mode='by_shank')

All recording objects through the preprocessing (filtering, referencing, etc.) have the same channel group information.

# rec objects
print(rec_probe.get_channel_groups())
print(rec_f.get_channel_groups())
print(rec_cmr.get_channel_groups())
print(rec_cashed.get_channel_groups())
# probe info
probe.to_dataframe()

shank_ids

alejoe91 commented 2 years ago

So you need to sort the probe dataframe by device_channel_index. Can you also show that column in the screenshot? Is it from 0 to N-1?

mfvd commented 2 years ago

Sorry, stupid error on my part. Forgot the sort the dataframe... It does make sense.

shank_ids_sorted