From histology to probemap.

[Gfernandezv]

Hello, I'm relatively new to sorting data, and I'm currently working on creating probes for analyzing my data. The records consist of four groups of four single electrodes (not tetrodes). Their positions are obtained through histological analysis, which provides the Anteroposterior (AP), Dorsoventral (DV), and Temporolateral (TL) stereotactic coordinates in millimeters. How can I convert these positions into probes?
As a preliminary step, I attempted to map AP to x, ML to y, and DV to z, but I notice that the probe interface uses Y as the depth (DV axis)

[zm711]

@Gfernandezv, A lot of the sorters don't play well with 3D probes in general. So typically the field uses planar probes where you have laterality on the x and depth on the y. Then you could add a 3 dimension if you add shanks in 3d, but then you often need to use other spikeinterface/probeinterface machinery to only do 2 dims at a time. Does that make sense for the broad issue?

That being said we can usually help more if you post your script that you've tried so far so we can point out specific places where you might need to adjust. You can definitely look over the general docs, but if we see your script we can help more.

[Gfernandezv]

Thank you for your response. After reading some papers, I agree with you. As a strategy, I have created a shank with four probes, each placed far enough apart to avoid interference during the analysis. Could this be a good strategy? This time, I attach the code.
Thanks in advance.
greetings.

%matplotlib widget
import numpy as np
import matplotlib.pyplot as plt

from probeinterface import Probe, ProbeGroup
from probeinterface.plotting import plot_probe
from probeinterface import generate_multi_shank, combine_probes

def create_probe(positions, ndim=2, si_units='um', shapes='circle', shape_params={'radius': 10}):
    probe = Probe(ndim=ndim, si_units=si_units)
    probe.set_contacts(positions=positions, shapes=shapes, shape_params=shape_params)
    probe.create_auto_shape(probe_type='tip')
    return probe
   
positions1 = np.array([
    [120.34, -708.82],
    [-2.80, -92.74],
    [26.19,	-76.18],
    [0, 0]
])

positions2 = np.array([
    [-137.03, -463.67],
    [-107.65, -375.10],
    [-71.06, -291.74],
    [0, 0]
])

positions3 = np.array([
    [0,	0],
    [22,0],
    [44,0],
    [66,0]
])

positions4 = np.array([
    [-29.13, -464.47],
    [42.42, -356.22],
    [-66.12, -268.36],
    [0.00, 0.00]
])

probe1 = create_probe(positions1)
probe2 = create_probe(positions2)
probe3 = create_probe(positions3)
probe4 = create_probe(positions4)

probe2.move([1000, 0])
probe3.move([2000, 0])
probe4.move([3000, 0])

multi_shank = combine_probes([probe1, probe2, probe3, probe4])

[zm711]

Yeah that could work. The other option would be to just use the run_sorter_by_property function which allows you to sort each "probe" that you've generate separately. But your strategy was what the kilosort repo (pre-KS4) would recommend for multishanks (ie for KS2, 2.5, and 3).

[Gfernandezv]

Your proposal seems better, will it be compatible with KS4?

[zm711]

Yes it is compatible. My only warning is that Kilosort in general hasn't always performed the best with low channel count probes their biorxiv so for low channel count they were recommending KS2, but you could definitely try KS4, but you will absolutely need to shut off drift correction since you don't have enough vertical spread of electric contacts.

[Gfernandezv]

Hi, it's me...again.
I been trying to implement the analysis by group, but I'm having some troubles of the type:

"ValueError: could not broadcast input array from shape (0,4,6) into shape (0,10,6)"

this is the code that generate the problem:

probe code:

%matplotlib widget
import numpy as np
import matplotlib.pyplot as plt
from probeinterface import Probe, ProbeGroup, write_probeinterface
from probeinterface.plotting import plot_probe, plot_probe_group

def create_probe(positions, ndim=2, si_units='um', shapes='circle', shape_params={'radius': 10}):
    probe = Probe(ndim=ndim, si_units=si_units)
    probe.set_contacts(positions=positions, shapes=shapes, shape_params=shape_params)
    return probe

n = 4
all_positions = []

# probes position:
for j in range(4):
    positions = np.zeros((n, 2))
    for i in range(n):
        x = j * 1000  # Avanzar x en saltos de 100
        y = i * 40   # Avanzar y en saltos de 40
        positions[i] = [x, y]
    all_positions.append(positions)

# setup probes
Anillo_probegroup = ProbeGroup()
for k, positions in enumerate(all_positions):
    probe = create_probe(positions)
    channel_indices = np.arange(k * n, (k + 1) * n)
    probe.set_contact_ids(channel_indices)
    probe.set_device_channel_indices(channel_indices)
    probe.set_shank_ids(np.zeros(4))
    Anillo_probegroup.add_probe(probe)

# checking.
Anillo_probegroup.to_dataframe()
print(Anillo_probegroup.get_global_device_channel_indices())
fig, ax = plt.subplots()
plot_probe_group(Anillo_probegroup, with_device_index=True, with_contact_id=True, same_axes=True, ax=ax)
plt.show()

#saving.
write_probeinterface('Anillo_probegroup.json', Anillo_probegroup)

now, sorting (for the example, I use KS4)...
a code for concatenation is included, because I create more than one file sometimes.

#importar librerias
%matplotlib widget
import os, shutil
import glob
import numpy as np
from pathlib import Path
import matplotlib.pyplot as plt
import spikeinterface.full as si
import probeinterface as pi
import spikeinterface.sorters as ss
import spikeinterface.extractors as se
import spikeinterface.preprocessing as prep

#parameters
n_cpus = os.cpu_count()
n_jobs = n_cpus
global_job_kwargs = dict(n_jobs=n_jobs, progress_bar=True, chunk_duration="1s")
si.set_global_job_kwargs(**global_job_kwargs)

# probe
probegroup = pi.read_probeinterface("Anillo_probegroup.json")

# rhd file list
rhd_files = [
    str(r'C:\Users\LabPF\OneDrive - Universidad Católica de Chile\Anillo\Records\NE01\Sleep\Dia01\NE01_Dia01_sleep_230515_114413\NE01_Dia01_sleep_230515_114413.rhd')
]

# detected files:
for rhd_file in rhd_files:
    print(f"Archivo detectado: {rhd_file}")

# read and concatenate.
raw_rec = [se.read_intan(file, stream_name='RHD2000 amplifier channel') for file in rhd_files]
raw_rec = si.concatenate_recordings(raw_rec)
raw_rec = raw_rec.set_probegroup(probegroup, group_mode='by_probe')
recording_f = prep.bandpass_filter(raw_rec, freq_min=300, freq_max=9000)

#sorting.
sorting_ks = ss.run_sorter_by_property(sorter_name='kilosort4',
                                       do_correction='False',
                                       recording=recording_f,
                                      grouping_property="group",
                                      working_folder="sort_by_group")

the error does not appear when I don't use groups:

sorting_ks = ss.run_sorter(sorter_name='kilosort4',
                                       do_correction='False',
                                       recording=recording_f)