translationalneuromodeling / tapas

TAPAS - Translational Algorithms for Psychiatry-Advancing Science
https://translationalneuromodeling.github.io/tapas/
GNU General Public License v3.0
221 stars 90 forks source link

Processing physio data in BIDS format #289

Open kmwag opened 4 days ago

kmwag commented 4 days ago

Hi everyone,

I’m currently using TAPAS PhysIO (R2022a-v8.1.0) to process physiological data recorded with the PERU_098 Physiological EGC/Respiratory Unit (Siemens) in a Siemens 3T Prisma scanner. The raw physiological recordings have already been preprocessed using BIDSmapper (BIDScoin, version 3.7.3), and the data is now in .tsv and .json formats for each participant and sequence (examples attached). My analysis focuses solely on the cardiac data—the respiratory signal is not required.

Below is the batch script I am currently using to extract heart rate data:

matlabbatch{1}.spm.tools.physio.save_dir` = {'C:\Users\wagne\OneDrive\Desktop\Kristin\Uni\Master\Masterarbeit\Auswertung\subs\sub-001\physio_out'};
matlabbatch{1}.spm.tools.physio.log_files.vendor = 'BIDS';
matlabbatch{1}.spm.tools.physio.log_files.cardiac = {'C:\Users\wagne\OneDrive\Desktop\Kristin\Uni\Master\Masterarbeit\Auswertung\subs\sub-001\sub-001_ses-mri01_task-1stScanSTRESSCControl_dir-AP_run-1_physio.tsv'};
matlabbatch{1}.spm.tools.physio.log_files.respiration = {''};
matlabbatch{1}.spm.tools.physio.log_files.scan_timing = {''};
matlabbatch{1}.spm.tools.physio.log_files.sampling_interval = 0.0025;
matlabbatch{1}.spm.tools.physio.log_files.relative_start_acquisition = [];
matlabbatch{1}.spm.tools.physio.log_files.align_scan = 'last';
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Nslices = 78;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.NslicesPerBeat = [];
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.TR = 1.2;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Ndummies = 0;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Nscans = 154;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.onset_slice = 1;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.time_slice_to_slice = [];
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Nprep = [];
matlabbatch{1}.spm.tools.physio.scan_timing.sync.nominal = struct([]);
matlabbatch{1}.spm.tools.physio.preproc.cardiac.modality = 'ECG';
matlabbatch{1}.spm.tools.physio.preproc.cardiac.filter.no = struct([]);
matlabbatch{1}.spm.tools.physio.preproc.cardiac.initial_cpulse_select.auto_matched.min = 0.4;
matlabbatch{1}.spm.tools.physio.preproc.cardiac.initial_cpulse_select.auto_matched.file = 'initial_cpulse_kRpeakfile.mat';
matlabbatch{1}.spm.tools.physio.preproc.cardiac.posthoc_cpulse_select.off = struct([]);
matlabbatch{1}.spm.tools.physio.model.output_multiple_regressors = 'multiple_regressors.txt';
matlabbatch{1}.spm.tools.physio.model.output_physio = 'physio.mat';
matlabbatch{1}.spm.tools.physio.model.orthogonalise = 'none';
matlabbatch{1}.spm.tools.physio.model.censor_unreliable_recording_intervals = false;
matlabbatch{1}.spm.tools.physio.model.retroicor.yes.order.c = 3;
matlabbatch{1}.spm.tools.physio.model.retroicor.yes.order.r = 4;
matlabbatch{1}.spm.tools.physio.model.retroicor.yes.order.cr = 1;
matlabbatch{1}.spm.tools.physio.model.rvt.no = struct([]);
matlabbatch{1}.spm.tools.physio.model.hrv.yes.delays = 0;
matlabbatch{1}.spm.tools.physio.model.noise_rois.no = struct([]);
matlabbatch{1}.spm.tools.physio.model.movement.no = struct([]);
matlabbatch{1}.spm.tools.physio.model.other.no = struct([]);
matlabbatch{1}.spm.tools.physio.verbose.level = 2;
matlabbatch{1}.spm.tools.physio.verbose.fig_output_file = '';
matlabbatch{1}.spm.tools.physio.verbose.use_tabs = false;

When running this script, MATLAB produces the following output:

29-Nov-2024 14:16:23 - Running 'TAPAS PhysIO Toolbox'
TR = 0.003 +/- 0.000 s (Estimated mean +/- std time of repetition for one volume; nTriggers = 156351)
No cardiac R-peak (heartbeat) events provided
    maxscan (incl. dummies) = 154 
    tmin (1st scan start (1st dummy))= 206.50 s
    tmin (1st scan start (after dummies))= 206.50 s
    tmax = 391.28 s 
    mean TR =   1.20 s

Note: Guessed 1 additional cardiac pulse(s) at time series end for phase estimation
29-Nov-2024 14:16:31 - Done    'TAPAS PhysIO Toolbox'
29-Nov-2024 14:16:31 - Done

Additionally, 7 output figures were generated (see attached). These include the filtered signal, R-peak detection diagnostics, and heart rate plots.

Questions:

  1. Does the MATLAB output indicate plausible results for heart rate extraction?
  2. Is this batch script properly configured, or are there adjustments I should make to improve data accuracy (e.g., preprocessing steps like filtering)?
  3. Can I confidently proceed with these outputs for my analysis?

Thanks in advance for your advice! 😊

Kristin

physio_files.zip output_figures.zip

mrikasper commented 1 day ago

Dear Kristin,

Thank you for trying out PhysIO.

I have some questions for clarification before I can answer yours:

  1. Do you know whether the original physiological logfile data was based on Siemens IdeaCmdTool data (e.g., .puls or .ecg data) or Siemens_Tics files or the CMRR multi-band sequence physiological loggin? See our wiki to distinguish the two.
  2. Do you happen to know which Siemens software release (e.g,. VE, XA30, XA60) your scanner uses? That might help answer the previous question?

The reason why I am asking is because you rely on BIDSCoin to convert your data to BIDS, and a critical aspect in that process is the synchronization of your physiological recordings with your image (DICOM/NIfTI) data. I don't know how BIDSCoin manages the synchronization, but in PhysIO, it's one of the most elaborate aspects in the processing.

If the synchronization is fine, the rest of the processing looks good to me, and the output figures (in particular the diagnostics figure 5 with the temporal lag between heartbeats) gives reasonable values (around 60 bpm). Additional filtering would only be necessary if your recordings drifted a lot over time.

If you want, you can send me the original logfiles in Siemens format and I can check whether PhysIO arrives at the same synchronization times.

I hope this helps!

All the best, Lars