Skip to content

Commit

Permalink
Merge pull request #22 from catalystneuro/FiberPhotometry/time-alignment
Browse files Browse the repository at this point in the history
[Fiber photometry] Add fiber photometry time alignment
  • Loading branch information
weiglszonja authored Dec 5, 2024
2 parents df0b12b + 7d43275 commit c9ec88a
Show file tree
Hide file tree
Showing 3 changed files with 229 additions and 36 deletions.
147 changes: 117 additions & 30 deletions src/constantinople_lab_to_nwb/fiber_photometry/fiber_photometry_notes.md
Original file line number Diff line number Diff line change
Expand Up @@ -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`:

Expand All @@ -72,15 +72,15 @@ 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
- `Values`: the isosbestic signal

`AIN02xAOUT02-LockIn`:
- `Time`: the time in seconds
- `Values`: motion corrected signal
- `Values`: ???

### Fiber photometry metadata

Expand All @@ -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:
Expand All @@ -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

Original file line number Diff line number Diff line change
@@ -1,3 +1,4 @@
from datetime import timedelta
from pathlib import Path

from neuroconv import NWBConverter
Expand All @@ -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):
Expand All @@ -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
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand All @@ -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)
Expand Down

0 comments on commit c9ec88a

Please sign in to comment.