v0.7.2 (1) updated USV summary to classify as noise spectrograms <.3 …

…correlation, (2) updated binary file conncatenation to add session start to every tracking start/end, (3) added ransac as Anipose default, (4) updated calibration board lengths to correct values, (5) updated e-phys sync and concatenation to read META file correctly (channel, LF, and AP), (6) added more waiting time at ethernet reconnection in behavioral_experiments

README.md

@@ -1,12 +1,12 @@
-# usv-playpen v0.7.1
+# usv-playpen v0.7.2
 
 <br>
 <div align="center">
   <img src="src/img/usv_playpen_gui.png">
 </div>
 <hr>
 
-GUI to facilitate conducting experiments with multi-probe e-phys (Neuropixels), multi-channel audio (Avisoft) and multi-camera video (Loopbio) acquisition. Developed for behavioral recording purposes at the [Princeton Neuroscience Institute](https://pni.princeton.edu/) 2021-24 (Falkner/Murthy labs). Due to proprietary software design and limitations, recordings can only be performed on OS Windows. The data processing branch of the GUI is platform-independent.
+GUI to facilitate conducting experiments with multi-probe e-phys (Neuropixels), multi-channel audio (Avisoft) and multi-camera video (Loopbio) acquisition. Developed for behavioral recording purposes at the [Princeton Neuroscience Institute](https://pni.princeton.edu/) 2021-25 (Falkner/Murthy labs). Due to proprietary software design and limitations, recordings can only be performed on OS Windows. The data processing branch of the GUI is platform-independent.
 
 [![Python version](https://img.shields.io/badge/Python-3.10-blue)](https://img.shields.io/badge/Python-3.10-blue)
 [![DOI](https://zenodo.org/badge/566588932.svg)](https://zenodo.org/badge/latestdoi/566588932)
@@ -87,4 +87,4 @@ Run the GUI.
 python usv_playpen_gui.py
 ```
 
-Developed and tested in PyCharm Pro 2024.2.3, on macOS Sonoma 14.4 / Pop!_OS 22.04 / Windows 11.
+Developed and tested in PyCharm Pro 2024.3, on macOS Sequoia 15.1 / Pop!_OS 22.04 / Windows 11.

setup.py

@@ -5,7 +5,7 @@
 
 setup(
     name='usv-playpen',
-    version='0.7.1',
+    version='0.7.2',
     author='@bartulem',
     author_email='mimica.bartul@gmail.com',
     classifiers=[

src/anipose_operations.py

@@ -313,6 +313,7 @@ def conduct_anipose_triangulation(self):
                                       calib=calibration_toml_file,
                                       frames=tuple(self.input_parameter_dict['conduct_anipose_triangulation']['frame_restriction']),
                                       excluded_views=tuple(self.input_parameter_dict['conduct_anipose_triangulation']['excluded_views']),
+                                      ransac=True,
                                       fname=f"{self.session_root_joint_date_dir}{os.sep}{self.session_root_name}_points3d.h5",
                                       disp_progress=self.input_parameter_dict['conduct_anipose_triangulation']['display_progress_bool'],
                                       constraints=self.input_parameter_dict['conduct_anipose_triangulation']['rigid_body_constraints'],
@@ -431,21 +432,21 @@ def translate_rotate_metric(self):
         z_theta_extra = -math.pi / 4
         arena_data = rotate_z(arena_data_temp, z_theta_extra)
 
-        # correct corners to be the middle of rail instead of corner
-        arena_data[0, 0, arena_nodes.index('East'), :] = [arena_data[0, 0, arena_nodes.index('East'), 0] + .0125,
-                                                          arena_data[0, 0, arena_nodes.index('East'), 1] + .0125,
+        # correct corners to be the outer edge of rail instead of inside corner
+        arena_data[0, 0, arena_nodes.index('East'), :] = [arena_data[0, 0, arena_nodes.index('East'), 0] + .025,
+                                                          arena_data[0, 0, arena_nodes.index('East'), 1] + .025,
                                                           0]
 
-        arena_data[0, 0, arena_nodes.index('West'), :] = [arena_data[0, 0, arena_nodes.index('West'), 0] - .0125,
-                                                          arena_data[0, 0, arena_nodes.index('West'), 1] - .0125,
+        arena_data[0, 0, arena_nodes.index('West'), :] = [arena_data[0, 0, arena_nodes.index('West'), 0] - .025,
+                                                          arena_data[0, 0, arena_nodes.index('West'), 1] - .025,
                                                           0]
 
-        arena_data[0, 0, arena_nodes.index('North'), :] = [arena_data[0, 0, arena_nodes.index('North'), 0] - .0125,
-                                                           arena_data[0, 0, arena_nodes.index('North'), 1] + .0125,
+        arena_data[0, 0, arena_nodes.index('North'), :] = [arena_data[0, 0, arena_nodes.index('North'), 0] - .025,
+                                                           arena_data[0, 0, arena_nodes.index('North'), 1] + .025,
                                                            0]
 
-        arena_data[0, 0, arena_nodes.index('South'), :] = [arena_data[0, 0, arena_nodes.index('South'), 0] + .0125,
-                                                           arena_data[0, 0, arena_nodes.index('South'), 1] - .0125,
+        arena_data[0, 0, arena_nodes.index('South'), :] = [arena_data[0, 0, arena_nodes.index('South'), 0] + .025,
+                                                           arena_data[0, 0, arena_nodes.index('South'), 1] - .025,
                                                            0]
 
         if self.input_parameter_dict['translate_rotate_metric']['save_arena_data_bool']:
@@ -475,6 +476,11 @@ def translate_rotate_metric(self):
             mouse_data_temp = mouse_data.copy()
             mouse_data = rotate_z(mouse_data_temp, z_theta_extra)
 
+            # set all negative mouse z-coordinates to 0
+            negative_z_indices = np.where(mouse_data[:, :, :, 2] < 0)
+            if negative_z_indices[0].size > 0:
+                mouse_data[negative_z_indices[0], negative_z_indices[1], negative_z_indices[2], 2] = 0
+
             mouse_track_names = find_mouse_names(root_directory=self.root_directory)
 
             with h5py.File(name=f"{self.session_root_joint_date_dir}{os.sep}{self.session_root_name}_points3d_translated_rotated_metric.h5", mode='w') as h5_file_write:

src/behavioral_experiments.py

@@ -123,12 +123,12 @@ def check_camera_vitals(self, camera_fr=None):
 
         try:
             api = motifapi.MotifApi(ip_address, api_key)
+            available_cameras = api.call('cameras')['cameras']
         except motifapi.api.MotifError:
             self.message_output('Motif not running or reachable. Check hardware and connections.')
             QTest.qWait(10000)
             sys.exit()
         else:
-            available_cameras = api.call('cameras')['cameras']
             self.camera_serial_num = [camera_dict['serial'] for camera_dict in available_cameras]
             self.message_output(f"The system is running Motif v{api.call('version')['software']} "
                                 f"and {len(available_cameras)} camera(s) is/are online: {self.camera_serial_num}")
@@ -460,7 +460,7 @@ def conduct_behavioral_recording(self):
                              stdout=subprocess.DEVNULL,
                              stderr=subprocess.STDOUT).wait()
             self.message_output(f"Ethernet RECONNECTED at {datetime.datetime.now().hour:02d}:{datetime.datetime.now().minute:02d}.{datetime.datetime.now().second:02d}.")
-            QTest.qWait(10000)
+            QTest.qWait(20000)
 
         self.message_output(f"Transferring audio/video files started at: {datetime.datetime.now().hour:02d}:{datetime.datetime.now().minute:02d}.{datetime.datetime.now().second:02d}")
 

src/das_inference.py

@@ -6,6 +6,8 @@
 from PyQt6.QtTest import QTest
 import glob
 import json
+import librosa
+import matplotlib.pyplot as plt
 import numpy as np
 import os
 import pathlib
@@ -247,17 +249,42 @@ def summarize_das_findings(self):
             usv_summary['duration'] = usv_summary['stop'] - usv_summary['start']
             usv_summary.sort_values(by='start', ascending=True, inplace=True)
 
-            # find peak and mean amplitude channels
+            # find peak and mean amplitude channels and filter out noise
+            mean_signal_correlations = np.zeros(usv_summary.shape[0])
+            mean_signal_correlations[:] = np.nan
             audio_file_loc = sorted(glob.glob(f"{self.root_directory}{os.sep}audio{os.sep}hpss_filtered{os.sep}*.mmap"))[0]
             audio_file_name = os.path.basename(audio_file_loc)
             data_type, channel_num, sample_num, audio_sampling_rate = audio_file_name.split("_")[-1][:-5], int(audio_file_name.split("_")[-2]), int(audio_file_name.split("_")[-3]), int(audio_file_name.split("_")[-4])
             audio_file_data = np.memmap(filename=audio_file_loc, mode='r', dtype=data_type, shape=(sample_num, channel_num))
+
+            len_win_signal = 512
+            low_freq_cutoff = 30000
+            noise_corr_cutoff = 0.3
+            frequency_resolution = audio_sampling_rate / len_win_signal
+            lower_bin = int(np.floor(low_freq_cutoff / frequency_resolution))
+
             for index, row in tqdm(usv_summary.iterrows(), desc="USV clean-up progress", total=usv_summary.shape[0], position=0, leave=True):
                 start_usv = int(np.floor(row['start'] * audio_sampling_rate))
                 stop_usv = int(np.ceil(row['stop'] * audio_sampling_rate))
                 peak_amp_ch = np.unravel_index(np.argmax(audio_file_data[start_usv:stop_usv, :]), audio_file_data.shape)[1]
                 mean_amp_ch = np.argmax(np.abs(audio_file_data[start_usv:stop_usv, :]).mean(axis=0))
-                usv_detected_chs = list(row['chs_detected'])
+                usv_detected_chs = row['chs_detected']
+
+                # the following section computes channel-wise signal correlations in the frequency domain
+                # if the signal has a mean channel-wise correlation of less than 0.3 (or 0.4 for <4 channels!), it is likely noise
+                if len(usv_detected_chs) > 1:
+                    usv_data_selected_ch = audio_file_data[start_usv:stop_usv, usv_detected_chs].astype('float32').T
+                    spectrogram_data_selected_ch = np.abs(librosa.stft(usv_data_selected_ch, n_fft=len_win_signal))
+                    reshaped_spectrogram = spectrogram_data_selected_ch[:, lower_bin:, :].reshape(len(usv_detected_chs), -1)
+                    correlation_matrix = np.corrcoef(reshaped_spectrogram)
+                    unique_correlations = correlation_matrix[np.triu_indices(n=len(usv_detected_chs), k=1)]
+                    mean_signal_correlations[index] = np.mean(unique_correlations)
+                    if len(usv_detected_chs) > 3:
+                        condition_4 = np.mean(unique_correlations) < noise_corr_cutoff
+                    else:
+                        condition_4 = np.mean(unique_correlations) < (noise_corr_cutoff + 0.1)
+                else:
+                    condition_4 = False
 
                 # remove USV segments if they appear only on one channel; this gets rid of some true positives, but false positives to a larger extent
                 condition_1 = len(usv_detected_chs) < 2
@@ -266,12 +293,25 @@ def summarize_das_findings(self):
                 # if the USV is detected on two channels, but the peak amplitude channel is not the same as the mean amplitude channel, it is likely noise
                 condition_3 = len(usv_detected_chs) == 2 and peak_amp_ch != mean_amp_ch
 
-                if condition_1 or condition_2 or condition_3:
+                if condition_1 or condition_2 or condition_3 or condition_4:
                     usv_summary.drop(labels=index, inplace=True)
                 else:
                     usv_summary.at[index, 'peak_amp_ch'] = peak_amp_ch
                     usv_summary.at[index, 'mean_amp_ch'] = mean_amp_ch
 
+            fig, ax = plt.subplots(nrows=1, ncols=1)
+            ax.hist(x=mean_signal_correlations,
+                    bins=20,
+                    histtype='stepfilled',
+                    color='#BBD5E8',
+                    ec='#000000',
+                    alpha=.5)
+            ax.set_xlabel('Mean signal/spectral correlation')
+            ax.set_ylabel('Number of putative USVs')
+            ax.axvline(x=noise_corr_cutoff, ls='-.', lw=1.2, c='#000000')
+            fig.savefig(f"{self.root_directory}{os.sep}audio{os.sep}{session_id}_usv_signal_correlation_histogram.svg", dpi=300)
+            plt.close()
+
             # give ID number to each USV
             usv_summary['usv_id'] = [f"{_num:04d}" for _num in range(usv_summary.shape[0])]
 

src/file_manipulation.py

@@ -143,7 +143,8 @@ def split_clusters_to_sessions(self):
                 spike_times = np.load(f"{ephys_dir}{os.sep}kilosort{self.input_parameter_dict['get_spike_times']['kilosort_version']}{os.sep}spike_times.npy")
 
                 if phy_curation_bool:
-                    cluster_info = pd.read_csv(f"{ephys_dir}{os.sep}kilosort{self.input_parameter_dict['get_spike_times']['kilosort_version']}{os.sep}cluster_info.tsv", sep='\t')
+                    cluster_info = pd.read_csv(filepath_or_buffer=f"{ephys_dir}{os.sep}kilosort{self.input_parameter_dict['get_spike_times']['kilosort_version']}{os.sep}cluster_info.tsv",
+                                               sep='\t')
                 else:
                     self.message_output("Phy2 curation has not been done for this session, no cluster_info file exists.")
                     continue
@@ -241,7 +242,7 @@ def concatenate_binary_files(self):
         for probe_idx, probe_id in enumerate(available_probes):
             binary_files_info = {}
             changepoints = [0]
-            concatenation_command = 'copy /b '
+            concatenation_command = 'copy /b ' if os.name == 'nt' else 'cat '
             for ord_idx, one_root_dir in enumerate(self.root_directory):
                 if os.path.isdir(f'{one_root_dir}{os.sep}ephys{os.sep}{probe_id}'):
                     for one_file, one_meta in zip(sorted(list(pathlib.Path(f'{one_root_dir}{os.sep}ephys{os.sep}{probe_id}').glob(f"*{npx_file_type}.bin*"))),
@@ -253,7 +254,7 @@ def concatenate_binary_files(self):
                                 for line in meta_data_file:
                                     key, value = line.strip().split("=")
                                     if key == 'acqApLfSy':
-                                        total_num_channels = int(value.split(',')[-1]) + int(value.split(',')[-2])
+                                        total_num_channels = int(value.split(',')[0]) + int(value.split(',')[-1])
                                     elif key == 'imDatHs_sn':
                                         headstage_sn = value
                                         spike_glx_sr = float(calibrated_sr_config['CalibratedHeadStages'][headstage_sn])
@@ -285,9 +286,15 @@ def concatenate_binary_files(self):
                                 binary_files_info[binary_file_info_id]['session_start_end'][0] = changepoints[-1]
                                 binary_files_info[binary_file_info_id]['session_start_end'][1] = int(one_recording.shape[0] // total_num_channels) + changepoints[-1]
                                 changepoints.append(int(one_recording.shape[0] // total_num_channels) + changepoints[-1])
-                                concatenation_command += '+ {} '.format(one_file)
+                                if os.name == 'nt':
+                                    concatenation_command += '+ {} '.format(one_file)
+                                else:
+                                    concatenation_command += '{} '.format(one_file)
 
-            concatenation_command += f'"{concat_save_dir[probe_idx]}{os.sep}concatenated_{concat_save_dir[probe_idx].split(os.sep)[-1]}.{npx_file_type}.bin"'
+            if os.name == 'nt':
+                concatenation_command += f'"{concat_save_dir[probe_idx]}{os.sep}concatenated_{concat_save_dir[probe_idx].split(os.sep)[-1]}.{npx_file_type}.bin"'
+            else:
+                concatenation_command += f'> {concat_save_dir[probe_idx]}{os.sep}concatenated_{concat_save_dir[probe_idx].split(os.sep)[-1]}.{npx_file_type}.bin'
 
             # create save directory if one doesn't exist already
             pathlib.Path(concat_save_dir[probe_idx]).mkdir(parents=True, exist_ok=True)
@@ -307,6 +314,9 @@ def concatenate_binary_files(self):
                         for component_key in changepoint_info_data[file_key].keys():
                             if component_key != 'tracking_start_end' and component_key != 'root_directory' and changepoint_info_data[file_key][component_key] != binary_files_info[file_key][component_key]:
                                 changepoint_info_data[file_key][component_key] = binary_files_info[file_key][component_key]
+                            elif component_key == 'tracking_start_end' and changepoint_info_data[file_key][component_key] != [np.nan, np.nan]:
+                                changepoint_info_data[file_key][component_key][0] = changepoint_info_data[file_key][component_key][0] + binary_files_info[file_key]['session_start_end'][0]
+                                changepoint_info_data[file_key][component_key][1] = changepoint_info_data[file_key][component_key][1] + binary_files_info[file_key]['session_start_end'][0]
 
                 with open(f'{concat_save_dir[probe_idx]}{os.sep}changepoints_info_{concat_save_dir[probe_idx].split(os.sep)[-1]}.json', 'w') as binary_info_output_file:
                     json.dump(changepoint_info_data, binary_info_output_file, indent=4)

src/preprocessing_plot.py

@@ -324,5 +324,6 @@ def preprocessing_summary(self, ipi_discrepancy_dict, phidget_data_dictionary):
                 axin6.set_yticks([])
                 axin6.set_xlabel('time (s)')
 
-        fig.savefig(f"{self.root_directory}{os.sep}sync{os.sep}{self.root_directory.split(os.sep)[-1]}_summary.png",
+        fig.savefig(f"{self.root_directory}{os.sep}sync{os.sep}{self.root_directory.split(os.sep)[-1]}_summary.svg",
                     dpi=300)
+        plt.close()