Merge pull request #55 from int-brain-lab/develop

Develop

Author: mayofaulkner

atlaselectrophysiology/alignment_with_easyqc.py

@@ -1,7 +1,7 @@
 from easyqc.gui import viewseis
 from ibllib.dsp import voltage
 from ibllib.ephys import neuropixel
-from viewspikes.data import stream, get_ks2
+from viewspikes.data import stream
 from viewspikes.plots import overlay_spikes
 import scipy
 from PyQt5 import QtCore, QtGui
@@ -35,6 +35,10 @@ def __init__(self, probe_id=None, one=None, histology=False, spike_collection=No
         self.lines_tracks = []
         self.points = []
 
+        self.plotdata.channels = Bunch()
+        self.plotdata.channels['localCoordinates'] = self.plotdata.chn_coords_all
+        self.plotdata.channels['rawInd'] = self.plotdata.chn_ind_all
+
     def on_mouse_double_clicked(self, event):
         if not self.offline:
             if event.double() and event.modifiers() and QtCore.Qt.ShiftModifier:
@@ -136,12 +140,11 @@ def stream_ap(self, t):
         sos = scipy.signal.butter(3, 300 / self.ap.fs / 2, btype='highpass', output='sos')
         butt = scipy.signal.sosfiltfilt(sos, raw)
         destripe = voltage.destripe(raw, fs=self.ap.fs)
-        ks2 = get_ks2(raw, dsets, self.loaddata.one)
+
         self.eqc['butterworth'] = viewseis(butt.T, si=1 / self.ap.fs, h=h, t0=t0, title='butt',
                                            taxis=0)
         self.eqc['destripe'] = viewseis(destripe.T, si=1 / self.ap.fs, h=h, t0=t0, title='destr',
                                         taxis=0)
-        self.eqc['ks2'] = viewseis(ks2.T, si=1 / self.ap.fs, h=h, t0=t0, title='ks2', taxis=0)
 
         overlay_spikes(self.eqc['butterworth'], self.plotdata.spikes, self.plotdata.clusters,
                        self.plotdata.channels)
@@ -227,10 +230,19 @@ def add_clust_scatter(self):
 def load_extra_data(probe_id, one=None, spike_collection=None):
     one = one or ONE()
     eid, probe = one.pid2eid(probe_id)
-    if spike_collection:
+    if spike_collection == '':
+        collection = f'alf/{probe}'
+    elif spike_collection:
         collection = f'alf/{probe}/{spike_collection}'
     else:
-        collection = f'alf/{probe}'
+        # Pykilosort is default, if not present look for normal kilosort
+        # Find all collections
+        all_collections = one.list_collections(eid)
+
+        if f'alf/{probe}/pykilosort' in all_collections:
+            collection = f'alf/{probe}/pykilosort'
+        else:
+            collection = f'alf/{probe}'
 
     _ = one.load_object(eid, obj='spikes', collection=collection,
                         attribute='samples')

atlaselectrophysiology/ephys_atlas_gui.py

@@ -1,3 +1,10 @@
+import os
+import platform
+
+if platform.system() == 'Darwin':
+    if platform.release().split('.')[0] == '20':
+        os.environ["QT_MAC_WANTS_LAYER"] = "1"
+
 from PyQt5 import QtWidgets, QtCore, QtGui
 import pyqtgraph as pg
 import pyqtgraph.exporters
@@ -11,7 +18,6 @@
 import atlaselectrophysiology.ephys_gui_setup as ephys_gui
 from atlaselectrophysiology.create_overview_plots import make_overview_plot
 from pathlib import Path
-import os
 import qt
 import matplotlib.pyplot as mpl  # noqa  # This is needed to make qt show properly :/
 
@@ -33,19 +39,20 @@ def _get_or_create(title='Electrophysiology Atlas', **kwargs):
         return av
 
     def __init__(self, offline=False, probe_id=None, one=None, histology=True,
-                 spike_collection=None):
+                 spike_collection=None, remote=False):
         super(MainWindow, self).__init__()
 
         self.init_variables()
         self.init_layout(self, offline=offline)
         self.configure = True
+        one_mode = 'remote' if remote else 'auto'
         if not offline and probe_id is None:
-            self.loaddata = LoadData()
+            self.loaddata = LoadData(mode=one_mode)
             self.populate_lists(self.loaddata.get_subjects(), self.subj_list, self.subj_combobox)
             self.offline = False
         elif not offline and probe_id is not None:
             self.loaddata = LoadData(probe_id=probe_id, one=one, load_histology=histology,
-                                     spike_collection=spike_collection)
+                                     spike_collection=spike_collection, mode=one_mode)
             self.current_shank_idx = 0
             _, self.histology_exists = self.loaddata.get_info(0)
             self.feature_prev, self.track_prev = self.loaddata.get_starting_alignment(0)
@@ -1266,6 +1273,11 @@ def data_button_pressed(self):
             self.probe_rfmap, self.rfmap_boundaries = self.plotdata.get_rfmap_data()
             self.img_stim_data = self.plotdata.get_passive_events()
 
+            if not self.offline:
+                self.img_raw_data = self.plotdata.get_raw_data_image(self.loaddata.probe_id, one=self.loaddata.one)
+            else:
+                self.img_raw_data = {}
+
             if self.histology_exists:
                 self.slice_data = self.loaddata.get_slice_images(self.ephysalign.xyz_samples)
             else:
@@ -2106,11 +2118,12 @@ def viewer(probe_id, one=None, histology=False, spike_collection=None, title=Non
 
     parser = argparse.ArgumentParser(description='Offline vs online mode')
     parser.add_argument('-o', '--offline', default=False, required=False, help='Offline mode')
+    parser.add_argument('-r', '--remote', default=False, required=False, action='store_true', help='Remote mode')
     parser.add_argument('-i', '--insertion', default=None, required=False, help='Insertion mode')
     args = parser.parse_args()
 
     app = QtWidgets.QApplication([])
-    mainapp = MainWindow(offline=args.offline, probe_id=args.insertion)
+    mainapp = MainWindow(offline=args.offline, probe_id=args.insertion, remote=args.remote)
     # mainapp = MainWindow(offline=True)
     mainapp.show()
     app.exec_()

atlaselectrophysiology/ephys_gui_setup.py

@@ -105,6 +105,14 @@ def init_menubar(self):
         img_options.addAction(scatter_amp)
         self.img_options_group.addAction(scatter_amp)
 
+        raw_type = list(self.img_raw_data.keys())
+        for raw in raw_type:
+            img = QtWidgets.QAction(raw, self, checkable=True, checked=False)
+            img.triggered.connect(lambda checked, item=raw: self.plot_image(
+                                  self.img_raw_data[item]))
+            img_options.addAction(img)
+            self.img_options_group.addAction(img)
+
         stim_type = list(self.img_stim_data.keys())
         for stim in stim_type:
             img = QtWidgets.QAction(stim, self, checkable=True, checked=False)

atlaselectrophysiology/example_code.py

@@ -6,6 +6,3 @@
 
 
 av = viewer(probe_id, one=one)
-
-
-# To add trials to the window

atlaselectrophysiology/extract_files.py

@@ -82,20 +82,21 @@ def extract_rmsmap(fbin, out_folder=None, spectra=True):
         out_folder = Path(fbin).parent
     else:
         out_folder = Path(out_folder)
-    alf_object_time = f'_iblqc_ephysTimeRms{sglx.type.upper()}'
-    alf_object_freq = f'_iblqc_ephysSpectralDensity{sglx.type.upper()}'
+    alf_object_time = f'ephysTimeRms{sglx.type.upper()}'
+    alf_object_freq = f'ephysSpectralDensity{sglx.type.upper()}'
 
     # crunch numbers
     rms = rmsmap(fbin, spectra=spectra)
     # output ALF files, single precision with the optional label as suffix before extension
     if not out_folder.exists():
         out_folder.mkdir()
     tdict = {'rms': rms['TRMS'].astype(np.single), 'timestamps': rms['tscale'].astype(np.single)}
-    alfio.save_object_npy(out_folder, object=alf_object_time, dico=tdict)
+    alfio.save_object_npy(out_folder, object=alf_object_time, dico=tdict, namespace='iblqc')
     if spectra:
         fdict = {'power': rms['spectral_density'].astype(np.single),
                  'freqs': rms['fscale'].astype(np.single)}
-        alfio.save_object_npy(out_folder, object=alf_object_freq, dico=fdict)
+        alfio.save_object_npy(
+            out_folder, object=alf_object_freq, dico=fdict, namespace='iblqc')
 
 
 def _sample2v(ap_file):

atlaselectrophysiology/load_data.py

@@ -19,8 +19,8 @@
 
 class LoadData:
     def __init__(self, one=None, brain_atlas=None, testing=False, probe_id=None,
-                 load_histology=True, spike_collection=None):
-        self.one = one or ONE(base_url=ONE_BASE_URL)
+                 load_histology=True, spike_collection=None, mode='auto'):
+        self.one = one or ONE(base_url=ONE_BASE_URL, mode=mode)
         self.brain_atlas = brain_atlas or atlas.AllenAtlas(25)
         self.download_hist = load_histology  # whether or not to look for the histology files
         self.spike_collection = spike_collection
@@ -435,7 +435,7 @@ def upload_data(self, xyz_channels, channels=True):
             # Create new trajectory and overwrite previous one
             histology.register_aligned_track(self.probe_id, xyz_channels,
                                              chn_coords=self.chn_coords, one=self.one,
-                                             overwrite=True, channels=channels)
+                                             overwrite=True, channels=channels, brain_atlas=self.brain_atlas)
         else:
             channel_upload = False
 

atlaselectrophysiology/load_data_local.py

@@ -6,6 +6,11 @@
 import glob
 import json
 from one.api import ONE
+from atlaselectrophysiology.load_histology import tif2nrrd
+
+# temporarily add this in for neuropixel course until figured out fix to problem on win32
+import ssl
+ssl._create_default_https_context = ssl._create_unverified_context
 
 
 class LoadDataLocal:
@@ -133,15 +138,21 @@ def get_slice_images(self, xyz_channels):
         # First see if the histology file exists before attempting to connect with FlatIron and
         # download
 
-        path_to_rd_image = glob.glob(str(self.folder_path) + '/*RD.nrrd')
-        if path_to_rd_image:
-            hist_path_rd = Path(path_to_rd_image[0])
+        path_to_rd_image_nrrd = glob.glob(str(self.folder_path) + '/*RD.nrrd')
+        path_to_rd_image_tif = glob.glob(str(self.folder_path) + '/*RD.tif')
+        if path_to_rd_image_nrrd:
+            hist_path_rd = Path(path_to_rd_image_nrrd[0])
+        elif path_to_rd_image_tif:
+            hist_path_rd = tif2nrrd(path_to_rd_image_tif[0])
         else:
             hist_path_rd = []
 
-        path_to_gr_image = glob.glob(str(self.folder_path) + '/*GR.nrrd')
-        if path_to_gr_image:
-            hist_path_gr = Path(path_to_gr_image[0])
+        path_to_gr_image_nrrd = glob.glob(str(self.folder_path) + '/*GR.nrrd')
+        path_to_gr_image_tif = glob.glob(str(self.folder_path) + '/*GR.tif')
+        if path_to_gr_image_nrrd:
+            hist_path_gr = Path(path_to_gr_image_nrrd[0])
+        elif path_to_gr_image_tif:
+            hist_path_gr = tif2nrrd(path_to_gr_image_tif[0])
         else:
             hist_path_gr = []
 
@@ -244,11 +255,11 @@ def create_channel_dict(brain_regions):
         channel_dict = {}
         for i in np.arange(brain_regions.id.size):
             channel = {
-                'x': brain_regions.xyz[i, 0] * 1e6,
-                'y': brain_regions.xyz[i, 1] * 1e6,
-                'z': brain_regions.xyz[i, 2] * 1e6,
-                'axial': brain_regions.axial[i],
-                'lateral': brain_regions.lateral[i],
+                'x': np.float64(brain_regions.xyz[i, 0] * 1e6),
+                'y': np.float64(brain_regions.xyz[i, 1] * 1e6),
+                'z': np.float64(brain_regions.xyz[i, 2] * 1e6),
+                'axial': np.float64(brain_regions.axial[i]),
+                'lateral': np.float64(brain_regions.lateral[i]),
                 'brain_region_id': int(brain_regions.id[i]),
                 'brain_region': brain_regions.acronym[i]
             }

atlaselectrophysiology/load_histology.py

@@ -51,7 +51,7 @@ def download_histology_data(subject, lab):
         path_to_image = Path(CACHE_DIR, file)
         if not path_to_image.exists():
             url = (baseurl + '/' + file)
-            http_download_file(url, cache_dir=CACHE_DIR,
+            http_download_file(url, target_dir=CACHE_DIR,
                                username=par.HTTP_DATA_SERVER_LOGIN,
                                password=par.HTTP_DATA_SERVER_PWD)
 
@@ -65,7 +65,7 @@ def download_histology_data(subject, lab):
 
 
 def tif2nrrd(path_to_image):
-    path_to_nrrd = Path(path_to_image.parent, path_to_image.parts[-1][:-3] + 'nrrd')
+    path_to_nrrd = Path(path_to_image).with_suffix('.nrrd')
     if not path_to_nrrd.exists():
         reader = sitk.ImageFileReader()
         reader.SetImageIO("TIFFImageIO")

atlaselectrophysiology/plot_data.py

@@ -3,8 +3,10 @@
 import numpy as np
 import one.alf as alf
 from brainbox.processing import bincount2D
+from brainbox.io.spikeglx import stream
 from brainbox.population.decode import xcorr
 from brainbox.task import passive
+from ibllib.dsp import voltage
 import scipy
 from PyQt5 import QtGui
 
@@ -23,8 +25,9 @@ def __init__(self, probe_path, ephys_path, alf_path, shank_idx):
         self.ephys_path = ephys_path
         self.alf_path = alf_path
 
-        self.chn_coords_all = np.load(Path(self.probe_path, 'channels.localCoordinates.npy'))
-        self.chn_ind_all = np.load(Path(self.probe_path, 'channels.rawInd.npy'))
+        channels = alf.io.load_object(self.probe_path, 'channels')
+        self.chn_coords_all = channels['localCoordinates']
+        self.chn_ind_all = channels['rawInd'].astype(int)
 
         self.chn_min = np.min(self.chn_coords_all[:, 1])
         self.chn_max = np.max(self.chn_coords_all[:, 1])
@@ -71,7 +74,8 @@ def __init__(self, probe_path, ephys_path, alf_path, shank_idx):
             self.cluster_data_status = True
             self.compute_timescales()
 
-        except Exception:
+        except Exception as err:
+            print(err)
             print('cluster data was not found, some plots will not display')
             self.cluster_data_status = False
 
@@ -295,9 +299,11 @@ def get_fr_img(self):
         else:
             T_BIN = 0.05
             D_BIN = 5
+            chn_min = np.min(np.r_[self.chn_min, self.spikes['depths'][self.spike_idx][self.kp_idx]])
+            chn_max = np.max(np.r_[self.chn_max, self.spikes['depths'][self.spike_idx][self.kp_idx]])
             n, times, depths = bincount2D(self.spikes['times'][self.spike_idx][self.kp_idx],
                                           self.spikes['depths'][self.spike_idx][self.kp_idx],
-                                          T_BIN, D_BIN, ylim=[self.chn_min, self.chn_max])
+                                          T_BIN, D_BIN, ylim=[chn_min, chn_max])
             img = n.T / T_BIN
             xscale = (times[-1] - times[0]) / img.shape[0]
             yscale = (depths[-1] - depths[0]) / img.shape[1]
@@ -323,14 +329,16 @@ def get_fr_amp_data_line(self):
         else:
             T_BIN = np.max(self.spikes['times'])
             D_BIN = 10
+            chn_min = np.min(np.r_[self.chn_min, self.spikes['depths'][self.spike_idx][self.kp_idx]])
+            chn_max = np.max(np.r_[self.chn_max, self.spikes['depths'][self.spike_idx][self.kp_idx]])
             nspikes, times, depths = bincount2D(self.spikes['times'][self.spike_idx][self.kp_idx],
                                                 self.spikes['depths'][self.spike_idx][self.kp_idx],
                                                 T_BIN, D_BIN,
-                                                ylim=[self.chn_min, self.chn_max])
+                                                ylim=[chn_min, chn_max])
 
             amp, times, depths = bincount2D(self.spikes['amps'][self.spike_idx][self.kp_idx],
                                             self.spikes['depths'][self.spike_idx][self.kp_idx],
-                                            T_BIN, D_BIN, ylim=[self.chn_min, self.chn_max],
+                                            T_BIN, D_BIN, ylim=[chn_min, chn_max],
                                             weights=self.spikes['amps'][self.spike_idx]
                                             [self.kp_idx])
             mean_fr = nspikes[:, 0] / T_BIN
@@ -363,9 +371,11 @@ def get_correlation_data_img(self):
         else:
             T_BIN = 0.05
             D_BIN = 40
+            chn_min = np.min(np.r_[self.chn_min, self.spikes['depths'][self.spike_idx][self.kp_idx]])
+            chn_max = np.max(np.r_[self.chn_max, self.spikes['depths'][self.spike_idx][self.kp_idx]])
             R, times, depths = bincount2D(self.spikes['times'][self.spike_idx][self.kp_idx],
                                           self.spikes['depths'][self.spike_idx][self.kp_idx],
-                                          T_BIN, D_BIN, ylim=[self.chn_min, self.chn_max])
+                                          T_BIN, D_BIN, ylim=[chn_min, chn_max])
             corr = np.corrcoef(R)
             corr[np.isnan(corr)] = 0
             scale = (np.max(depths) - np.min(depths)) / corr.shape[0]
@@ -464,6 +474,38 @@ def median_subtract(a):
 
         return data_img, data_probe
 
+    # only for IBL sorry
+    def get_raw_data_image(self, pid, t0=(1000, 2000, 3000), one=None):
+
+        def gain2level(gain):
+            return 10 ** (gain / 20) * 4 * np.array([-1, 1])
+        data_img = dict()
+        for t in t0:
+
+            sr, t = stream(pid, t, one=one)
+            raw = sr[:, :-sr.nsync].T
+            channel_labels, channel_features = voltage.detect_bad_channels(raw, sr.fs)
+            raw = voltage.destripe(raw, fs=sr.fs, channel_labels=channel_labels)
+            raw_image = raw[:, int((450 / 1e3) * sr.fs):int((500 / 1e3) * sr.fs)].T
+            x_range = np.array([0, raw_image.shape[0] - 1]) / sr.fs * 1e3
+            levels = gain2level(-90)
+            xscale = (x_range[1] - x_range[0]) / raw_image.shape[0]
+            yscale = (self.chn_max - self.chn_min) / raw_image.shape[1]
+
+            data_raw = {
+                'img': raw_image,
+                'scale': np.array([xscale, yscale]),
+                'levels': levels,
+                'offset': np.array([0, 0]),
+                'cmap': 'bone',
+                'xrange': x_range,
+                'xaxis': 'Time (ms)',
+                'title': 'Power (uV)'
+            }
+            data_img[f'Raw data t={t}'] = data_raw
+
+        return data_img
+
     def get_lfp_spectrum_data(self):
         freq_bands = np.vstack(([0, 4], [4, 10], [10, 30], [30, 80], [80, 200]))
         data_probe = {}