Error when reading physio data in BIDS format

[alexsayal]

Dear all,

I am using PhysIO (version R2021a-v8.0.1) to process cardiac and respiratory data in BIDS format. The two signals have different sampling rates and as such are supplied in BIDS in two separate files (in attachment): sub-01_task-AA_acq-0500_run-01_recording-cardiac_physio.tsv.gz sub-01_task-AA_acq-0500_run-01_recording-respiratory_physio.tsv.gz

However, the code exists with an error related to failed peak detection:

TR = 0.500 +/- 0.000 s (Estimated mean +/- std time of repetition for one volume; nTriggers = 780)
No cardiac R-peak (heartbeat) events provided
18-Nov-2021 17:08:07 - Failed  'TAPAS PhysIO Toolbox'
Error using tapas_physio_log (line 56)
No peaks found in raw cardiac time series. Check raw physiological recordings figure whether there is any non-constantcardiac data
In file "/home/alexandresayal/Documents/MATLAB/spm12/toolbox/code-physio/utils/tapas_physio_log.m" (???), function "tapas_physio_log" at line 56.
In file "/home/alexandresayal/Documents/MATLAB/spm12/toolbox/code-physio/preproc/tapas_physio_get_cardiac_pulse_template.m" (???), function "tapas_physio_get_cardiac_pulse_template" at line 120.
In file "/home/alexandresayal/Documents/MATLAB/spm12/toolbox/code-physio/preproc/tapas_physio_get_cardiac_pulses_auto_matched.m" (???), function "tapas_physio_get_cardiac_pulses_auto_matched" at line 55.
In file "/home/alexandresayal/Documents/MATLAB/spm12/toolbox/code-physio/preproc/tapas_physio_get_cardiac_pulses.m" (???), function "tapas_physio_get_cardiac_pulses" at line 87.
In file "/home/alexandresayal/Documents/MATLAB/spm12/toolbox/code-physio/tapas_physio_main_create_regressors.m" (???), function "tapas_physio_main_create_regressors" at line 169.
In file "/home/alexandresayal/Documents/MATLAB/spm12/toolbox/code-physio/tapas_physio_cfg_matlabbatch.m" (???), function "run_physio" at line 1568.

This is the batch script I am using:

matlabbatch{1}.spm.tools.physio.save_dir = {'/home/alexandresayal/Desktop/testBIDSPhysIO'};
matlabbatch{1}.spm.tools.physio.log_files.vendor = 'BIDS';
matlabbatch{1}.spm.tools.physio.log_files.cardiac = {'/home/alexandresayal/Desktop/testBIDSPhysIO/sub-01_task-AA_acq-0500_run-01_recording-cardiac_physio.tsv.gz'};
matlabbatch{1}.spm.tools.physio.log_files.respiration = {'/home/alexandresayal/Desktop/testBIDSPhysIO/sub-01_task-AA_acq-0500_run-01_recording-respiratory_physio.tsv.gz'};
matlabbatch{1}.spm.tools.physio.log_files.scan_timing = {''};
matlabbatch{1}.spm.tools.physio.log_files.sampling_interval = [];
matlabbatch{1}.spm.tools.physio.log_files.relative_start_acquisition = 0;
matlabbatch{1}.spm.tools.physio.log_files.align_scan = 'first';
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Nslices = 42;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.NslicesPerBeat = [];
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.TR = 0.5;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Ndummies = 0;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Nscans = 780;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.onset_slice = 21;
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 = 'PPU';
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.no = struct([]);
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;

I tried with different files and changing the sampling rates of the signals manually, with no success. Could you please help me detect what is going on?

Thank you in advance.

physio-files.zip

[mrikasper]

Dear Alexandre,

there are NaN values in your cardiac time series (you will see some gaps in the raw data plot, I marked some with yellow in the screenshot). It is easy to fix that in the PhysIO code, but only if one knows whether they are excess values (i.e. remove them) or actual timepoints that were not recorded (i.e., replace them e.g., by a next-neighbours mean).

Do you have any idea where they could come from? Since they are not at extremal points (min/max) of the timeseries, I have no immediate idea what would cause such a behavior. What device did you use for the cardiac recordings?

All the best, Lars

[alexsayal]

Hello Lars,

Thank you for spotting this. It seems like it is the .log to .tsv converter I am using that is introducing the NaNs. If I use the .log files directly no errors occur. The recording device is the PPU from Siemens.

I will check this converter code first and then get back to you.

All the best,

[mrikasper]

Dear Alexandre,

if this is Siemens data, you can also read it in directly with PhysIO, instead of converting it to BIDS first. We can read most of the common Siemens software releases like VB, VD, VE, and I am working on some bugfixes for XA30 at the moment.

If the only reason to use a different tool for conversion is that you have to store BIDS data for archiving later on, please let me know. It's on my todo-list for a while now to write out the preprocessed raw timeseries data in BIDS-format within PhysIO as well. Maybe I get some time (or help :-)) during the upcoming Brainhack Global to do so.

All the best, Lars

[alexsayal]

Hello Lars,

I am using BIDS here since I am building on the preprocessing of fMRIPrep. You are right, I can always use the .log files directly if this approach does not work, but I wanted to try with BIDS-formatted files first.

Nevertheless, the NaNs have a reason - the .log files have in fact some missing data (missing tics): see that for instance between these two consecutive samples of the PULS file

     21583117     PULS   2060 
     21583121     PULS   2030 

the ACQ_TIME_TICS increases by 4 and not 2 as expected. The code seems to be handling this when we provide the .log but not when we provide the .tsv.

I attach the original .log files. physio_files_log.zip

I don't know why this happened, but these seem like actual timepoints that were not recorded and could hence be replaced by meaningful values. Is this easy to implement?

All the best,

[mrikasper]

Dear Alexandre,

it's great to hear that the PhysIO reader for the native Siemens VE format works!

It is very common for this format to skip certain time points (ticks), I presume it's a buffering issue. The scanner is a real time system with a lot of safety measures, so if anything is more urgent than writing out physiological recordings (eg estimating SAR or executing gradient waveforms), this operation probably gets demoted. But I honestly didn't find any good documentation on that, there was a lot of reverse engineering involved on our side here to implement this reader.

What you encountered, is unfortunately a limitation of the BIDS format for physiological logfiles, in that it assumes constant sampling rates (instead of eg saving a time vector). So there is inherent information loss here, unless they update their specifications.

Internally, what PhysIO does in these cases for Siemens files is interpolation, because these are very short time windows compared to the bandwidth of physiological signals. I can mimic the same behavior for NaNs in the BIDS, and update the code, but it might take some time.

Which tool do you use for the BIDS conversion? You could probably let them know that the NaNs are unexpected behavior, maybe they are open to change it to interpolation or at least produce a warning.

To understand your workflow better, are the Siemens logfiles stored in different places than the BIDS ones, i.e., you cannot use them in your FMRIprep workflow with PhysIO, because you want to keep the data folders strictly BIDS? If that's the case, all the more reason for me to write out BIDS physio files within our toolbox, because then the conversion happens in a consistent fashion.

All the best, Lars

[alexsayal]

Hello Lars,

I require a valid BIDS directory for fmriPrep but not for the next stages, so no, I don't necessarily need the physio files in BIDS format for now. And I can still use the derivatives folder to organize these files and move on.

I am using bidsphysio (https://github.com/cbinyu/bidsphysio) for the conversion. I recently opened an issue regarding the NaNs (https://github.com/cbinyu/bidsphysio/issues/7) and agreed that it was a good option here since they allow the user to deal with those values as wanted.

I dig into the function tapas_physio_read_physlogfiles_bids.m and made the import of the .tsv files possible (with different sampling rates for cardiac and resp signals) and the output confound matrix is very similar to what the toolbox returns with the .log files.

tapas_physio_read_physlogfiles_bids.m.txt

I used linear interpolation to replace the NaNs and upscaled the respiratory signal, replicating what happens when using the .log files.

What you do think, would this be valid?

Have a nice weekend!