Merge branch 'run_kilosort' of https://github.com/sidhulyalkar/element-array-ephys

Author: sidhulyalkar

CHANGELOG.md

@@ -28,6 +28,11 @@ Observes [Semantic Versioning](https://semver.org/spec/v2.0.0.html) standard and
 + Update - Move common functions to `element-interface` package.
 + Add - NWB export function
 
+## [0.1.0b4] - 2021-11-29
+### Added
++ Processing with Kilosort and pyKilosort for Open Ephys and SpikeGLX
+
+
 ## [0.1.0b0] - 2021-05-07
 
 + Update - First beta release

element_array_ephys/__init__.py

@@ -1,2 +1,14 @@
-# ephys_acute as default
-import element_array_ephys.ephys_acute as ephys
+import datajoint as dj
+import logging
+import os
+
+
+dj.config['enable_python_native_blobs'] = True
+
+
+def get_logger(name):
+    log = logging.getLogger(name)
+    log.setLevel(os.getenv('LOGLEVEL', 'INFO'))
+    return log
+
+from . import ephys_acute as ephys

element_array_ephys/ephys_acute.py

@@ -18,7 +18,7 @@ def activate(schema_name, *, create_schema=True, create_tables=True):
     schema.activate(schema_name, create_schema=create_schema, create_tables=create_tables)
 
     # Add neuropixels probes
-    for probe_type in ('neuropixels 1.0 - 3A', 'neuropixels 1.0 - 3B',
+    for probe_type in ('neuropixels 1.0 - 3A', 'neuropixels 1.0 - 3B', 'neuropixels UHD',
                        'neuropixels 2.0 - SS', 'neuropixels 2.0 - MS'):
         ProbeType.create_neuropixels_probe(probe_type)
 
@@ -46,15 +46,38 @@ class Electrode(dj.Part):
     def create_neuropixels_probe(probe_type='neuropixels 1.0 - 3A'):
         """
         Create `ProbeType` and `Electrode` for neuropixels probes:
-         1.0 (3A and 3B), 2.0 (SS and MS)
+        + neuropixels 1.0 - 3A
+        + neuropixels 1.0 - 3B
+        + neuropixels UHD
+        + neuropixels 2.0 - SS
+        + neuropixels 2.0 - MS
+
         For electrode location, the (0, 0) is the
          bottom left corner of the probe (ignore the tip portion)
         Electrode numbering is 1-indexing
         """
 
+        neuropixels_probes_config = {
+            'neuropixels 1.0 - 3A': dict(site_count=960, col_spacing=32, row_spacing=20,
+                                         white_spacing=16, col_count=2,
+                                         shank_count=1, shank_spacing=0),
+            'neuropixels 1.0 - 3B': dict(site_count=960, col_spacing=32, row_spacing=20,
+                                         white_spacing=16, col_count=2,
+                                         shank_count=1, shank_spacing=0),
+            'neuropixels UHD': dict(site_count=384, col_spacing=6, row_spacing=6,
+                                    white_spacing=0, col_count=8,
+                                    shank_count=1, shank_spacing=0),
+            'neuropixels 2.0 - SS': dict(site_count=1280, col_spacing=32, row_spacing=15,
+                                         white_spacing=0, col_count=2,
+                                         shank_count=1, shank_spacing=250),
+            'neuropixels 2.0 - MS': dict(site_count=1280, col_spacing=32, row_spacing=15,
+                                         white_spacing=0, col_count=2,
+                                         shank_count=4, shank_spacing=250)
+        }
+
         def build_electrodes(site_count, col_spacing, row_spacing,
-                             white_spacing, col_count=2,
-                             shank_count=1, shank_spacing=250):
+                             white_spacing, col_count,
+                             shank_count, shank_spacing):
             """
             :param site_count: site count per shank
             :param col_spacing: (um) horrizontal spacing between sites
@@ -66,14 +89,15 @@ def build_electrodes(site_count, col_spacing, row_spacing,
             :return:
             """
             row_count = int(site_count / col_count)
-            x_coords = np.tile([0, 0 + col_spacing], row_count)
-            x_white_spaces = np.tile([white_spacing, white_spacing, 0, 0], int(row_count / 2))
+            x_coords = np.tile(np.arange(0, col_spacing * col_count, col_spacing), row_count)
+            y_coords = np.repeat(np.arange(row_count) * row_spacing, col_count)
 
-            x_coords = x_coords + x_white_spaces
-            y_coords = np.repeat(np.arange(row_count) * row_spacing, 2)
+            if white_spacing:
+                x_white_spaces = np.tile([white_spacing, white_spacing, 0, 0], int(row_count / 2))
+                x_coords = x_coords + x_white_spaces
 
-            shank_cols = np.tile([0, 1], row_count)
-            shank_rows = np.repeat(range(row_count), 2)
+            shank_cols = np.tile(range(col_count), row_count)
+            shank_rows = np.repeat(range(row_count), col_count)
 
             npx_electrodes = []
             for shank_no in range(shank_count):
@@ -88,51 +112,12 @@ def build_electrodes(site_count, col_spacing, row_spacing,
 
             return npx_electrodes
 
-        # ---- 1.0 3A ----
-        if probe_type == 'neuropixels 1.0 - 3A':
-            electrodes = build_electrodes(site_count=960, col_spacing=32, row_spacing=20,
-                                          white_spacing=16, col_count=2)
-
-            probe_type = {'probe_type': 'neuropixels 1.0 - 3A'}
-            with ProbeType.connection.transaction:
-                ProbeType.insert1(probe_type, skip_duplicates=True)
-                ProbeType.Electrode.insert([{**probe_type, **e} for e in electrodes],
-                                           skip_duplicates=True)
-
-        # ---- 1.0 3B ----
-        if probe_type == 'neuropixels 1.0 - 3B':
-            electrodes = build_electrodes(site_count=960, col_spacing=32, row_spacing=20,
-                                          white_spacing=16, col_count=2)
-
-            probe_type = {'probe_type': 'neuropixels 1.0 - 3B'}
-            with ProbeType.connection.transaction:
-                ProbeType.insert1(probe_type, skip_duplicates=True)
-                ProbeType.Electrode.insert([{**probe_type, **e} for e in electrodes],
-                                           skip_duplicates=True)
-
-        # ---- 2.0 Single shank ----
-        if probe_type == 'neuropixels 2.0 - SS':
-            electrodes = build_electrodes(site_count=1280, col_spacing=32, row_spacing=15,
-                                          white_spacing=0, col_count=2,
-                                          shank_count=1, shank_spacing=250)
-
-            probe_type = {'probe_type': 'neuropixels 2.0 - SS'}
-            with ProbeType.connection.transaction:
-                ProbeType.insert1(probe_type, skip_duplicates=True)
-                ProbeType.Electrode.insert([{**probe_type, **e} for e in electrodes],
-                                           skip_duplicates=True)
-
-        # ---- 2.0 Multi shank ----
-        if probe_type == 'neuropixels 2.0 - MS':
-            electrodes = build_electrodes(site_count=1280, col_spacing=32, row_spacing=15,
-                                          white_spacing=0, col_count=2,
-                                          shank_count=4, shank_spacing=250)
-
-            probe_type = {'probe_type': 'neuropixels 2.0 - MS'}
-            with ProbeType.connection.transaction:
-                ProbeType.insert1(probe_type, skip_duplicates=True)
-                ProbeType.Electrode.insert([{**probe_type, **e} for e in electrodes],
-                                           skip_duplicates=True)
+        electrodes = build_electrodes(**neuropixels_probes_config[probe_type])
+        probe_type = {'probe_type': probe_type}
+        with ProbeType.connection.transaction:
+            ProbeType.insert1(probe_type, skip_duplicates=True)
+            ProbeType.Electrode.insert([{**probe_type, **e} for e in electrodes],
+                                       skip_duplicates=True)
 
 
 @schema

element_array_ephys/ephys_chronic.py

@@ -12,7 +12,7 @@
 
 class Kilosort:
 
-    kilosort_files = [
+    _kilosort_core_files = [
         'params.py',
         'amplitudes.npy',
         'channel_map.npy',
@@ -22,47 +22,60 @@ class Kilosort:
         'similar_templates.npy',
         'spike_templates.npy',
         'spike_times.npy',
-        'spike_times_sec.npy',
-        'spike_times_sec_adj.npy',
         'template_features.npy',
         'template_feature_ind.npy',
         'templates.npy',
         'templates_ind.npy',
         'whitening_mat.npy',
         'whitening_mat_inv.npy',
-        'spike_clusters.npy',
+        'spike_clusters.npy'
+    ]
+
+    _kilosort_additional_files = [
+        'spike_times_sec.npy',
+        'spike_times_sec_adj.npy',
         'cluster_groups.csv',
         'cluster_KSLabel.tsv'
     ]
 
-    # keys to self.files, .data are file name e.g. self.data['params'], etc.
-    kilosort_keys = [path.splitext(kilosort_file)[0] for kilosort_file in kilosort_files]
+    kilosort_files = _kilosort_core_files + _kilosort_additional_files
 
     def __init__(self, kilosort_dir):
         self._kilosort_dir = pathlib.Path(kilosort_dir)
         self._files = {}
         self._data = None
         self._clusters = None
 
-        params_filepath = kilosort_dir / 'params.py'
-
-        if not params_filepath.exists():
-            raise FileNotFoundError(f'No Kilosort output found in: {kilosort_dir}')
+        self.validate()
 
+        params_filepath = kilosort_dir / 'params.py'
         self._info = {'time_created': datetime.fromtimestamp(params_filepath.stat().st_ctime),
                       'time_modified': datetime.fromtimestamp(params_filepath.stat().st_mtime)}
 
     @property
     def data(self):
         if self._data is None:
-            self._stat()
+            self._load()
         return self._data
 
     @property
     def info(self):
         return self._info
 
-    def _stat(self):
+    def validate(self):
+        """
+        Check if this is a valid set of kilosort outputs - i.e. all crucial files exist
+        """
+        missing_files = []
+        for f in Kilosort._kilosort_core_files:
+            full_path = self._kilosort_dir / f
+            if not full_path.exists():
+                missing_files.append(f)
+        if missing_files:
+            raise FileNotFoundError(f'Kilosort files missing in ({self._kilosort_dir}):'
+                                    f' {missing_files}')
+
+    def _load(self):
         self._data = {}
         for kilosort_filename in Kilosort.kilosort_files:
             kilosort_filepath = self._kilosort_dir / kilosort_filename
@@ -94,15 +107,15 @@ def _stat(self):
         self._data['channel_map'] = self._data['channel_map'].flatten()
 
         # Read the Cluster Groups
-        for cluster_pattern, cluster_col_name in zip(['cluster_groups.*', 'cluster_KSLabel.*'],
+        for cluster_pattern, cluster_col_name in zip(['cluster_group.*', 'cluster_KSLabel.*'],
                                                      ['group', 'KSLabel']):
             try:
                 cluster_file = next(self._kilosort_dir.glob(cluster_pattern))
             except StopIteration:
                 pass
             else:
                 cluster_file_suffix = cluster_file.suffix
-                assert cluster_file_suffix in ('.csv', '.tsv', '.xlsx')
+                assert cluster_file_suffix in ('.tsv', '.xlsx')
                 break
         else:
             raise FileNotFoundError(