diff --git a/src/constantinople_lab_to_nwb/fiber_photometry/fiber_photometry_notes.md b/src/constantinople_lab_to_nwb/fiber_photometry/fiber_photometry_notes.md index e852e86..2cd0555 100644 --- a/src/constantinople_lab_to_nwb/fiber_photometry/fiber_photometry_notes.md +++ b/src/constantinople_lab_to_nwb/fiber_photometry/fiber_photometry_notes.md @@ -52,16 +52,16 @@ From `J097_rDAgAChDMSDLS_20240820_HJJ_0000.doric`: `LockInAOUT02`: - `Time`: the time in seconds - - `AIN02`: motion corrected green signal (fiber 1) - - `AIN04`: motion corrected green signal (fiber 2) + - `AIN02`: estimated signal for green indicator (fiber 1) + - `AIN04`: estimated signal for green indicator (fiber 2) `LockInAOUT03`: - `Time`: the time in seconds - - `AIN01`: motion corrected mCherry signal (fiber 1) + - `AIN01`: estimated signal for mCherry (fiber 1) `LockInAOUT04`: - `Time`: the time in seconds - - `AIN03`: motion corrected mCherry signal (fiber 2) + - `AIN03`: estimated signal for mCherry (fiber 2) From `J069_ACh_20230809_HJJ_0002.doric`: @@ -72,7 +72,7 @@ TODO: what is the channel mapping here? `AIN01xAOUT02-LockIn`: - `Time`: the time in seconds - - `Values`: motion corrected green signal + - `Values`: ??? `AIN02xAOUT01-LockIn`: - `Time`: the time in seconds @@ -80,7 +80,7 @@ TODO: what is the channel mapping here? `AIN02xAOUT02-LockIn`: - `Time`: the time in seconds - - `Values`: motion corrected signal + - `Values`: ??? ### Fiber photometry metadata @@ -99,21 +99,18 @@ For each `FiberPhotometryResponseSeries` that we add to NWB, we need to specify Example: ```yaml FiberPhotometryResponseSeries: - - name: fiber_photometry_response_series - description: The raw fluorescence signal - channel_column_names: ["AIn-1 - Raw", "AIn-2 - Raw"] + - name: raw_fiber_photometry_signal + description: The raw fiber photometry signal before demodulation. + channel_column_names: ["AIn-1 - Raw", "AIn-2 - Raw", "AIn-3"] fiber_photometry_table_region: [0, 1] - fiber_photometry_table_region_description: The region of the FiberPhotometryTable corresponding to the raw fluorescence signal. - - name: fiber_photometry_response_series_isosbestic - description: The isosbestic signal - channel_column_names: ["AIn-3"] - fiber_photometry_table_region: [0] - fiber_photometry_table_region_description: The region of the FiberPhotometryTable corresponding to the isosbestic signal. - - name: fiber_photometry_response_series_motion_corrected - description: The motion corrected signal + unit: a.u. + fiber_photometry_table_region_description: The region of the FiberPhotometryTable corresponding to the raw signal. + - name: estimated_fiber_photometry_response_series + description: The demodulated (estimated) signal from light stimulation using a proprietary algorithm from Doric. channel_column_names: ["AIn-1 - Dem (AOut-1)", "AIn-2 - Dem (AOut-2)"] + unit: a.u. fiber_photometry_table_region: [0, 1] - fiber_photometry_table_region_description: The region of the FiberPhotometryTable corresponding to the motion corrected signal. + fiber_photometry_table_region_description: The region of the FiberPhotometryTable corresponding to the estimated signal. ``` The metadata file for the .doric files is named `doric_fiber_photometry_metadata.yaml` and it contains the following fields: @@ -128,19 +125,109 @@ For each `FiberPhotometryResponseSeries` that we add to NWB, we need to specify Example: ```yaml FiberPhotometryResponseSeries: - - name: fiber_photometry_response_series - description: The raw fluorescence signal # TBD + - name: raw_fiber_photometry_signal + description: The raw fiber photometry signal before demodulation. stream_names: ["AnalogIn/AIN01", "AnalogIn/AIN02", "AnalogIn/AIN03", "AnalogIn/AIN04"] - fiber_photometry_table_region: [0, 1, 2, 3] #[0, 1] - fiber_photometry_table_region_description: The region of the FiberPhotometryTable corresponding to the raw fluorescence signal. - - name: fiber_photometry_response_series_isosbestic - description: The isosbestic signal # TBD - stream_names: ["LockInAOUT01/AIN02", "LockInAOUT01/AIN04"] - fiber_photometry_table_region: [0, 2] - fiber_photometry_table_region_description: The region of the FiberPhotometryTable corresponding to the isosbestic signal. - - name: fiber_photometry_response_series_motion_corrected - description: The motion corrected signal # TBD + fiber_photometry_table_region: [0, 1, 2, 3] + fiber_photometry_table_region_description: The region of the FiberPhotometryTable corresponding to the raw signal. + unit: a.u. + - name: estimated_fiber_photometry_response_series + description: The demodulated (estimated) signal from light stimulation using a proprietary algorithm from Doric. stream_names: ["LockInAOUT03/AIN01", "LockInAOUT04/AIN03", "LockInAOUT02/AIN02", "LockInAOUT02/AIN04"] fiber_photometry_table_region: [0, 1, 2, 3] - fiber_photometry_table_region_description: The region of the FiberPhotometryTable corresponding to the motion corrected signal. + fiber_photometry_table_region_description: The region of the FiberPhotometryTable corresponding to the estimated signal. + unit: a.u. +``` + +Other example (J069_ACh_20230809_HJJ_0002.doric): + +Example: +```yaml + FiberPhotometryResponseSeries: + - name: raw_fiber_photometry_signal + description: The raw fiber photometry signal from Doric acquisition system before demodulation. + stream_names: ["AnalogIn/AIN01", "AnalogIn/AIN02"] + unit: a.u. + fiber_photometry_table_region: [0, 1] + fiber_photometry_table_region_description: The region of the FiberPhotometryTable corresponding to the raw signal. + - name: estimated_fiber_photometry_response_series + description: The demodulated (estimated) signal from light stimulation using a proprietary algorithm from Doric. + stream_names: ["AIN01xAOUT01-LockIn/Values", "AIN01xAOUT02-LockIn/Values", "AIN02xAOUT01-LockIn/Values", "AIN02xAOUT02-LockIn/Values"] + unit: a.u. + fiber_photometry_table_region: [0, 1, 0, 2] + fiber_photometry_table_region_description: The region of the FiberPhotometryTable corresponding to the estimated signal. + +``` + +### Session start time + +The session start time is the reference time for all timestamps in the NWB file. We are using `session_start_time` from the Bpod output. (The start time of the session in the Bpod data can be accessed from the "Info" struct, with "SessionDate" and "SessionStartTime_UTC" fields.) + +### Bpod trial start time + +We are extracting the trial start times from the Bpod output using the "TrialStartTimestamp" field. + +```python + from pymatreader import read_mat + + bpod_data = read_mat("raw_Bpod/J069/DataFiles/J069_RWTautowait2_20230809_131216.mat")["SessionData"] # should contain "SessionData" named struct + + # The trial start times from the Bpod data + bpod_trial_start_times = bpod_data['TrialStartTimestamp'] + ``` + +```python +bpod_trial_start_times[:7] +>>> [ 11.4261, 104.5276, 146.0112, 203.5646, 211.7232, 215.3226, 224.041 ] ``` + +### Doric trial start time + +The trial start times from the Doric acquisition can be obtained from one of the digital signals ("DigitalIO/DIO02" in .doric file, "DI/O-2" in .csv file). + +```python +import h5py +from neuroconv.tools.signal_processing import get_rising_frames_from_ttl + +doric_file = h5py.File("J069_ACh_20230809_HJJ_0002.doric", mode="r") +ttl_signal = doric_file["/DataAcquisition/FPConsole/Signals/Series0001/DigitalIO/DIO02"][:] +timestamps = doric_file["/DataAcquisition/FPConsole/Signals/Series0001/DigitalIO/Time"][:] + +rising_frames_from_center_port_ttl = get_rising_frames_from_ttl(ttl_signal) +num_trials = len(rising_frames_from_center_port_ttl) +doric_trial_start_times = [timestamps[rising_frames_from_center_port_ttl][i] for i in range(num_trials)] +``` + +```python +doric_trial_start_times[:7] +>>> [17.11626, 110.21736, 151.702835, 209.255035, 217.41393499999998, 221.01406, 229.73321] +``` + +## Alignment + +We are aligning the starting time of the fiber photometry (Doric), video and DLC interfaces to the Bpod interface. + +We are computing the time shift from the Bpod trial start time to the Doric trial start time. + +For example, the computed time shift for this session: +```python +time_shift = bpod_trial_start_times[0] - doric_trial_start_times[0] +>>> -5.6901600000000006 +``` + +We are applying this time_shift to the timestamps for the raw fluorescence signals as: + +```python +doric_timestamps = doric_file["/DataAcquisition/FPConsole/Signals/Series0001/AnalogIn/Time"][:] +aligned_timestamps = doric_timestamps + time_shift +``` + +1) When the time shift is negative and the first aligned timestamp of the doric trace is negative: +- shift back bpod (from every column that has a timestamp they have to be shifted back) +- shift back session start time +- don't have to move doric or video +2) When the time shift is negative and the first aligned timestamp of the doric trace is positive +- we move the doric and video backward +3) When time shift is positive +- we move the doric and video forward + diff --git a/src/constantinople_lab_to_nwb/fiber_photometry/fiber_photometry_nwbconverter.py b/src/constantinople_lab_to_nwb/fiber_photometry/fiber_photometry_nwbconverter.py index 17370f6..fb4b91a 100644 --- a/src/constantinople_lab_to_nwb/fiber_photometry/fiber_photometry_nwbconverter.py +++ b/src/constantinople_lab_to_nwb/fiber_photometry/fiber_photometry_nwbconverter.py @@ -1,3 +1,4 @@ +from datetime import timedelta from pathlib import Path from neuroconv import NWBConverter @@ -8,6 +9,8 @@ DoricCsvFiberPhotometryInterface, ) from constantinople_lab_to_nwb.general_interfaces import BpodBehaviorInterface +from neuroconv.datainterfaces.behavior.deeplabcut._dlc_utils import _get_movie_timestamps +from neuroconv.tools.signal_processing import get_rising_frames_from_ttl class FiberPhotometryNWBConverter(NWBConverter): @@ -20,3 +23,102 @@ class FiberPhotometryNWBConverter(NWBConverter): Video=VideoInterface, Behavior=BpodBehaviorInterface, ) + + def get_time_shift(self) -> float: + if "FiberPhotometryDoric" in self.data_interface_objects: + fiber_photometry_interface = self.data_interface_objects["FiberPhotometryDoric"] + digital_stream_names = ["DigitalIO/DIO01", "DigitalIO/DIO02"] + ttl_signals = fiber_photometry_interface._get_traces(stream_names=digital_stream_names) + timestamps = fiber_photometry_interface.get_original_timestamps(stream_name=digital_stream_names[0]) + else: + fiber_photometry_interface = self.data_interface_objects["FiberPhotometryCsv"] + digital_stream_names = ["DI/O-1", "DI/O-2"] + ttl_signals = fiber_photometry_interface._get_traces(channel_column_names=digital_stream_names) + timestamps = fiber_photometry_interface.get_original_timestamps() + + raw_behavior_interface = self.data_interface_objects["Behavior"] + trial_start_times_from_bpod, trial_stop_times_from_bpod = raw_behavior_interface.get_trial_times() + + # we are aligning the ttl signals to the raw bpod trial start times + rising_frames_from_center_port_ttl = get_rising_frames_from_ttl(ttl_signals[:, 1]) + center_port_aligned_times = timestamps[rising_frames_from_center_port_ttl] + + time_shift = trial_start_times_from_bpod[0] - center_port_aligned_times[0] + return time_shift + + def temporally_align_data_interfaces(self): + if "FiberPhotometryDoric" in self.data_interface_objects: + fiber_photometry_interface = self.data_interface_objects["FiberPhotometryDoric"] + digital_stream_names = ["DigitalIO/DIO01", "DigitalIO/DIO02"] + ttl_signals = fiber_photometry_interface._get_traces(stream_names=digital_stream_names) + timestamps = fiber_photometry_interface.get_original_timestamps(stream_name=digital_stream_names[0]) + + else: + fiber_photometry_interface = self.data_interface_objects["FiberPhotometryCsv"] + digital_stream_names = ["DI/O-1", "DI/O-2"] + ttl_signals = fiber_photometry_interface._get_traces(channel_column_names=digital_stream_names) + timestamps = fiber_photometry_interface.get_original_timestamps() + + raw_behavior_interface = self.data_interface_objects["Behavior"] + trial_start_times_from_bpod, trial_stop_times_from_bpod = raw_behavior_interface.get_trial_times() + + # we are aligning the ttl signals to the raw bpod trial start times + rising_frames_from_center_port_ttl = get_rising_frames_from_ttl(ttl_signals[:, 1]) + center_port_aligned_times = timestamps[rising_frames_from_center_port_ttl] + + time_shift = trial_start_times_from_bpod[0] - center_port_aligned_times[0] + self._time_shift = time_shift + + has_deep_lab_cut = "DeepLabCut" in self.data_interface_objects + has_video = "Video" in self.data_interface_objects + + rising_frames_from_camera_ttl = get_rising_frames_from_ttl(ttl_signals[:, 0]) + timestamps_from_camera_ttl = timestamps[rising_frames_from_camera_ttl] + + # When the time_shift is negative and the first aligned timestamp of the doric trace is negative: + if timestamps[0] + self._time_shift < 0: + aligned_trial_start_times = trial_start_times_from_bpod + self._time_shift + aligned_trial_stop_times = trial_stop_times_from_bpod + self._time_shift + raw_behavior_interface.set_aligned_trial_times( + trial_start_times=aligned_trial_start_times, + trial_stop_times=aligned_trial_stop_times, + ) + if has_deep_lab_cut and has_video: + video_interface = self.data_interface_objects["Video"] + movie_file_path = video_interface.source_data["file_paths"][0] + video_interface.set_aligned_segment_starting_times([timestamps_from_camera_ttl[0]]) + + dlc_interface = self.data_interface_objects["DeepLabCut"] + movie_timestamps = _get_movie_timestamps(movie_file=str(movie_file_path)) + aligned_dlc_timestamps = movie_timestamps + timestamps_from_camera_ttl[0] + dlc_interface.set_aligned_timestamps(aligned_timestamps=aligned_dlc_timestamps) + + elif timestamps[0] + self._time_shift > 0: + fiber_photometry_interface.set_aligned_starting_time(aligned_starting_time=self._time_shift) + + if has_deep_lab_cut and has_video: + video_interface = self.data_interface_objects["Video"] + movie_file_path = video_interface.source_data["file_paths"][0] + aligned_timestamps_camera = timestamps_from_camera_ttl + time_shift + video_interface.set_aligned_segment_starting_times([aligned_timestamps_camera[0]]) + + dlc_interface = self.data_interface_objects["DeepLabCut"] + movie_timestamps = _get_movie_timestamps(movie_file=str(movie_file_path)) + aligned_dlc_timestamps = movie_timestamps + aligned_timestamps_camera[0] + dlc_interface.set_aligned_timestamps(aligned_timestamps=aligned_dlc_timestamps) + + def get_metadata(self): + metadata = super().get_metadata() + + # Explicity set session_start_time to Bpod start time + session_start_time = self.data_interface_objects["Behavior"].get_metadata()["NWBFile"]["session_start_time"] + # When the time_shift is negative and the first aligned timestamp of the doric trace is negative + # we need to add the time_shift to the session_start_time + time_shift = self.get_time_shift() + if time_shift < 0: + # Time shift is provided in seconds, we need to get a timedelta and add it to the session_start_time + time_shift_timedelta = timedelta(seconds=time_shift) + session_start_time += time_shift_timedelta + + metadata["NWBFile"].update(session_start_time=session_start_time) + return metadata diff --git a/src/constantinople_lab_to_nwb/fiber_photometry/interfaces/doric_fiber_photometry_interface.py b/src/constantinople_lab_to_nwb/fiber_photometry/interfaces/doric_fiber_photometry_interface.py index dee3c5b..7080ec9 100644 --- a/src/constantinople_lab_to_nwb/fiber_photometry/interfaces/doric_fiber_photometry_interface.py +++ b/src/constantinople_lab_to_nwb/fiber_photometry/interfaces/doric_fiber_photometry_interface.py @@ -25,6 +25,7 @@ def __init__( ): super().__init__(file_path=file_path, verbose=verbose) self._timestamps = None + self._aligned_starting_time = None self._root_data_path = root_data_path data = self.load() if root_data_path not in data: @@ -46,13 +47,16 @@ def get_original_timestamps(self, stream_name: str) -> np.ndarray: raise ValueError(f"Time column '{self._time_column_name}' not found in '{stream_name}'.") return channel_group[self._time_column_name][:] - def get_timestamps(self, stream_name=str, stub_test: bool = False) -> np.ndarray: - timestamps = ( - self._timestamps if self._timestamps is not None else self.get_original_timestamps(stream_name=stream_name) - ) + def get_timestamps(self, stream_name: str, stub_test: bool = False) -> np.ndarray: + original_timestamps = self.get_original_timestamps(stream_name=stream_name) if stub_test: - return timestamps[:100] - return timestamps + original_timestamps = original_timestamps[:100] + if self._aligned_starting_time is not None: + return original_timestamps + self._aligned_starting_time + return original_timestamps + + def set_aligned_starting_time(self, aligned_starting_time: float) -> None: + self._aligned_starting_time = aligned_starting_time def set_aligned_timestamps(self, aligned_timestamps: np.ndarray) -> None: self._timestamps = np.array(aligned_timestamps)