SpikeInterface / spikeinterface

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

Exporting to phy changes cluster channels #3394

Open RikkelBob opened 1 month ago

RikkelBob commented 1 month ago

After manual curation, I'm recalculating cluster waveforms and some quality metrics using the script below. However, when I export to phy, the channel numbers are wrong. When I pause/debug the script before export_to_phy, channel numbers are correct, but inspecting cluster_info.tsv after exporting, they're different. I would fix this post-hoc, but because the channel numbers are wrong, the calculated waveforms and subsequent waveform based quality metrics are also incorrect. I'd very much appreciate any help!

                  rec_folder = Path(probe_folder, probe)
                        probe_number = find_numbers_in_string(probe)
                        phy_dir = rec_folder
                        ac_dir = Path(phy_dir, 'postprocessed')
                        ac_dir_wave = Path(ac_dir, 'waveforms')
                        ac_dir_phy = Path(ac_dir, 'phy')

                        if not os.path.isdir(ac_dir_phy):

                            waveform_folder = Path(phy_dir, "analyzer")
                            report_folder = Path(phy_dir, 'report')

                            # Load your sorted units from the Phy folder
                            sorting = read_phy(phy_dir)

                            # Get channel locations from .npy file
                            path_channel_location = Path(phy_dir, 'channel_positions.npy')
                            channel_locations = np.load(path_channel_location)
                            num_channels = channel_locations.shape[0]

                            # Load the original recording
                            recording = si.read_binary(Path(ses_folder, "ksData_probe" + str(probe_number[0]) + ".dat"),
                                                       sampling_frequency=fs,
                                                       num_channels=num_channels,
                                                       dtype="int16")

                            probe_64, locations = utils.gen_probe()
                            recording.set_probe(probe_64, in_place=True)

                            #recording.set_dummy_probe_from_locations(channel_locations)
                            recording.annotate(is_filtered=True)

                            sorting = remove_excess_spikes(sorting, recording)
                            sorting = remove_duplicated_spikes(sorting=sorting, censored_period_ms=0.5, method="keep_first")
                            #sorting = sorting.remove_empty_units()

                            print("building analyzer")
                            analyzer = si.create_sorting_analyzer(sorting=sorting, recording=recording, return_scaled=False)

                            analyzer.compute("random_spikes")
                            analyzer.compute("waveforms", save=True)
                            templates = analyzer.compute("templates", save=True)
                            analyzer.compute(
                                "spike_amplitudes",
                                peak_sign="pos"
                            )
                            analyzer.compute("noise_levels")

                            amp_cutoff = analyzer.compute(
                                "quality_metrics",
                                qm_params=dqm_params
                            )

                            print("Saving analyizer")
                            analyzer.save_as(format="binary_folder", folder=ac_dir_wave)
                            print("Exproting to phy")
                            export_to_phy(analyzer, output_folder=ac_dir_phy, copy_binary=False, additional_properties=['quality'], remove_if_exists=True)
                            print("Finished" + str(ac_dir_phy))
zm711 commented 1 month ago

Maybe related #3392? I'll have to read a bit more later. Are you using a probe group or probe?

RikkelBob commented 1 month ago

@zm711 I don't use probe groups, just a probe - recording.set_probe(probe_64, in_place=True)

zm711 commented 1 month ago

Internally we are representing probes as a probegroup (see here https://github.com/SpikeInterface/spikeinterface/blob/73f4d58d2d542e607155d49d30716cc28215e12f/src/spikeinterface/core/sortinganalyzer.py#L1101-L1106) so if there is some sort of bug happening when we get the probegroup (to get your probe) it could affect this as well. Do you want to try exporting from that PR and see if it fixes things or if it is still broken.

RikkelBob commented 1 month ago

I've tried it now, but the issue persists

zm711 commented 1 month ago

Okay thanks. We will read this more closely soon and see if we can find the place soon.

RikkelBob commented 3 weeks ago

Would love to hear it if anything has been found related to this :)

samuelgarcia commented 3 weeks ago

@alejoe91 could you check this?