Processing physio data in BIDS format

[kmwag]

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]

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

[kmwag]

Dear Lars,

Thank you for your reply! I took some time to check the answers to your questions, and here they are:

1. The data was recorded with the CMRR multi-band sequence.
2. The Siemens software release is VE11E.

Additionally, I’ve attached the original log files from the same participant for your reference. To generate the ECG, RESP, and INFO files, we ran a Python script to convert the raw DICOM file. Here’s the script we used:

#!/usr/bin/env` python
# -*- coding: utf-8 -*-
import argparse
import os

import numpy as np
import pydicom as dicom

PHYSIO_DATA_TAG = '7fe11010'


def main(path, outdir=None):
    '''Main function.'''
    if outdir is None:
        outdir, _ = os.path.split(path)
    d = dicom.dcmread(path)
    data = d[int(PHYSIO_DATA_TAG, 16)]
    raw = data.value
    length_bytes = len(raw)
    rows = d.AcquisitionNumber
    columns = length_bytes / rows
    num_files = columns / 1024
    length_file = int(len(raw) / num_files)
    arrs = [raw[i:i + length_file] for i in range(0, len(raw), length_file)]

    for arr in arrs:
        data_length = int(np.array(arr[:4]).view(np.uint32))
        fn_length = int(np.array(arr[4:8]).view(np.uint32))
        fn = arr[8:(8+fn_length)].decode('ascii')
        out_path = os.path.join(outdir, fn)
        print(fn)
        data = arr[1024:(1024 + data_length)]
        with open(out_path, 'wb') as f:
            f.write(data)


def parse_args():
    '''Parse arguments.'''
    parser = argparse.ArgumentParser(add_help=False)
    reqarg = parser.add_argument_group(title='required arguments')
    optarg = parser.add_argument_group(title='optional arguments')
    reqarg.add_argument('--physio', help='path of physio DICOM file',
                        required=True)
    optarg.add_argument('-h', '--help', action='help',
                        default=argparse.SUPPRESS,
                        help='Show this help message and exit')
    optarg.add_argument('--outdir', help='directory for output',
                        required=False)
    return parser.parse_args()


if __name__ == '__main__':
    args = parse_args()
    main(args.physio, args.outdir)

Let me know if there’s anything else I can clarify!

Best regards,

Kristin

siemens_physio_files.zip

[mrikasper]

Dear Kristin,

Thank you, that is very helpful.

1. The CMRR multiband sequence has the modern Siemens_Tics format, which is usually easier to synchronize than the ideaCmdTool version. Great!
2. Then the software version does not matter so much, VE11E is fine.
3. Thank you for sending the raw files in Siemens_Tics format. I reran both the original data and the BIDS-converted one with PhysIO (see attached batches mb_physio_bid_rename2m.txt and
   mb_physio_siemens_rename2m.txt
   ) and checked the synchronization via a zoomed in version of the cutout plot of the time courses that are used for preprocessing and modeling (see jpeg screenshots)
4. Good news: Synchronization (start of black curve cutout region and x-axis) is identical in both approaches, so you can happily use the BIDS-converted files.
5. There is some discrepancy in the respiratory curve

s (green-black curve), but you mentioned you are not interested in those.

• BIDS:
  

• Siemens_Tics:
  

I hope that helps!

All the best,
Lars

[kmwag]

Dear Lars,
Thank you for your efforts and for confirming. That sounds great! I’ll proceed with evaluating the BIDS files accordingly. Much appreciated!

Best regards,
Kristin