PhysIO with BIDS physio data as input and cardiac pulses, no analog about tapas HOT 15 CLOSED

Dear Hrvoje,

thank you for the BIDS format suggestion. It does not seem too different from the custom and the Siemens Tics format, so I will see about implementing a reader asap..

Concerning your example and script file I have the following questions inline:

Hello,
I'm trying to use the toolbox with my physio data recorded not by the scanner, but by the additional software called Spike2. The data I get is the trigger times, analog respiratory signal, but for the cardiac data I get only the heart rate 0/1 variable.

Do you mean you get a 1 at each peak of the pulse wave and 0 otherwise. Or do you get indeed a heartrate (averaged?) From the file, it looks like the former, but just wanted to be sure.

I converted all the data to BIDS format, physio data included. This essentially conforms the required inputs to the FSL's PNM (see the manual, section 8.6: pdf). Here are the sample files: physio_test.zip

How did you create the respiration.tsv, cardiac.tsv and trigger.tsv file from the BIDS *_physio.tsv?

Something like this(?)

fname = 'sub-s002_task-fnclearning_run-01_physio.json'
val = jsondecode(fileread(fname));
fid2 = fopen('sub-s002_task-fnclearning_run-01_physio.tsv');
C = textscan(fid2,'%f\t%f\t%f')
for i = 1:numel(val.Columns)
    fileName = sprintf('%s.tsv',val.Columns{i}); 
    X = C{i}; if i == 1, X = [X X]; end; % duplicate cardiac trigger time course
    save(fileName, '-ascii', '-tabs', 'X');
end

Since PNM involves a lot of manual work, I opted for PhysIO first. Do you have support for BIDS input files in the newer versions on github by any chance?

See above. Seems to be straightforward to implement, I will use your example as a template, if that is OK?

Although it seems PhysIO requires analog cardiac signal I gave it a try nevertheless by simply duplicating the cardiac 0/1 column and setting the following parameters

matlabbatch{1}.spm.tools.physio.preproc.cardiac.modality = 'PPU';
matlabbatch{1}.spm.tools.physio.preproc.cardiac.initial_cpulse_select.load_from_logfile = struct([]);
matlabbatch{1}.spm.tools.physio.preproc.cardiac.posthoc_cpulse_select.off = struct([]);
This seems to give sensible results, but I wanted to check with you before including them in fmri GLMs. Would RETROICOR computations be fine? What about HRV?

This seems correct. load_from_logfile will ignore the cardiac input, and if you have the (cardiac) trigger as 2nd column, this should work for both RETROICOR and HRV.

However, when I ran your job after my assumed file conversion above, I could see that the heart beat detection performed by Spike2 is suboptimal in this case. You can see this in the diagnostic plot PhysIO provides, here is the screenshot:

The spikes in the upper subplot (and the warning with the threshold) indicate that the heart cycle (peak-to-peak interval) jumps too drastically for a realistic heart rate. I would assume the peak detection erroneously re-interpreted side peaks as R-peaks (spiky minima) or missed occasional R-peaks (spike maxima).

Therefore, it would be better to retrieve the raw cardiac time series, and then use the auto_matched peak detection algorithm of PhysIO in this case.

Finally, When running a job, the toolbox gives an output: "No scan trigger events provided". I do have trigger data, but unsure where and how to provide it. I have placed it in physio.log_files.scan_timing but perhaps not in a correct format?

It might be that the scan triggers for the custom format are not properly integrated. You don't need them here anyway, if you specify the relative onset of first scan volume and physio logging as relative_start_acquisition field.

Here is the script I used:
matlabbatch{1}.spm.tools.physio.save_dir = {'/media/hstojic/dataneuro/fnclearning_fmri/dProcessed/physio_test'};
matlabbatch{1}.spm.tools.physio.log_files.vendor = 'Custom';
matlabbatch{1}.spm.tools.physio.log_files.cardiac = {'/media/hstojic/dataneuro/fnclearning_fmri/dProcessed/physio_test/cardiac.csv'};
matlabbatch{1}.spm.tools.physio.log_files.respiration = {'/media/hstojic/dataneuro/fnclearning_fmri/dProcessed/physio_test/respiratory.csv'};
matlabbatch{1}.spm.tools.physio.log_files.scan_timing = {'/media/hstojic/dataneuro/fnclearning_fmri/dProcessed/physio_test/trigger.csv'};
matlabbatch{1}.spm.tools.physio.log_files.sampling_interval = 0.02;
matlabbatch{1}.spm.tools.physio.log_files.relative_start_acquisition = 255.45;
matlabbatch{1}.spm.tools.physio.log_files.align_scan = 'first';
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Nslices = 16;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.NslicesPerBeat = [];
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.TR = 1.45;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Ndummies = 0;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Nscans = 474;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.onset_slice = 9;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.time_slice_to_slice = 0.090625;
matlabbatch{1}.spm.tools.physio.scan_timing.sqpar.Nprep = 0;
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.load_from_logfile = struct([]);
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.rvt.yes.delays = 0;
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 = '/media/hstojic/dataneuro/fnclearning_fmri/dProcessed/physio_test';
matlabbatch{1}.spm.tools.physio.verbose.use_tabs = false;
P.S. Its a multiband CMRR sequence, and this particular participant exited the scanner after this run so that's why it has a lengthy period with physio recordings after the last trigger

That's OK.

All the best,
Lars

from tapas.

Hi Lars,

support for BIDS format would be great - it does seem like it is becoming a standard for fMRI data sharing.

1. Do you mean you get a 1 at each peak of the pulse wave and 0 otherwise.

• Yes

2. How did you create the respiration.tsv, cardiac.tsv and trigger.tsv file from the BIDS *_physio.tsv?

• I actually use Python for everything, in this case I call matlab from a python script, through shell, so not sure if it makes sense to share my code, instead here are the resulting files:
  individual_files.zip

3. See above. Seems to be straightforward to implement, I will use your example as a template, if that is OK?

• that's fine, but check BIDS specification that I followed (link in the first message), to be on the safe side that I didn't implement things incorrectly
• note that BIDS actually implies analog heart recordings, not 0/1 heart pulse as I have, so you would need to account for that in your code - some might put analog data, some 0/1 like me :)
• note that BIDS specified StartTime in the json file is the opposite of how is specified in the PhysIO

4. Therefore, it would be better to retrieve the raw cardiac time series, and then use the auto_matched peak detection algorithm of PhysIO in this case.

• unfortunately, that data was not recorded :|
• that data is a bit noisy as the unit is attached to participants finger and they are using hands to make responses, so I expected some dropouts in the signal - I attributed these spikes to that
• do you think the data is usable as it is?

5. It might be that the scan triggers for the custom format are not properly integrated. You don't need them here anyway, if you specify the relative onset of first scan volume and physio logging as relative_start_acquisition field.

• it might be beneficial to use the actual triggers recorded in the physio data, some misalignments are I guess possible and actual triggers are probably more reliable than assuming exact TR after the first trigger (I think this is why FSL PNM toolbox requires, per their explanations in the user guide)

from tapas.

Hi Lars,
support for BIDS format would be great - it does seem like it is becoming a standard for fMRI data sharing.

Alright, I am on it :-)

How did you create the respiration.tsv, cardiac.tsv and trigger.tsv file from the BIDS *_physio.tsv?
I actually use Python for everything, in this case I call matlab from a python script, through shell, so not sure if it makes sense to share my code, instead here are the resulting files:
individual_files.zip

Good. Same result to what I did.

See above. Seems to be straightforward to implement, I will use your example as a template, if that is OK?

that's fine, but check BIDS specification that I followed (link in the first message), to be on the safe side that I didn't implement things incorrectly
note that BIDS actually implies analog heart recordings, not 0/1 heart pulse as I have, so you would need to account for that in your code - some might put analog data, some 0/1 like me :)
note that BIDS specified StartTime in the json file is the opposite of how is specified in the PhysIO

Both important points to consider. I will take care of that!

Therefore, it would be better to retrieve the raw cardiac time series, and then use the auto_matched peak detection algorithm of PhysIO in this case.

unfortunately, that data was not recorded :|
that data is a bit noisy as the unit is attached to participants finger and they are using hands to make responses, so I expected some dropouts in the signal - I attributed these spikes to that
do you think the data is usable as it is?

It could very well be the hand movements. Do subjects press with the same hand where the PPU is attached? I would try to avoid that with new data.

About usability of the data: Hard to say, I would give it a try, but in particular if you really think phys recording artifacts are correlated to subject response, you definitely have to check the F-contrasts (e.g., with tapas_physio_report_contrasts) for the cardiac regressors, so that they don't look, e.g., like motor activation.

One thing you could try is to remove the presumable double peaks that create the local minima peaks (e.g., 0.3 s at 361s) by removing the 1 there in the cardiac file. Not very principled w/o the raw data though. PhysIO has also the ability to flag and remove intervals of low reliability from the regressors. It's physio.model.censor_unreliable_recording_intervals = true; and see whether it detects the right regions. Not sure whether it works in your case though, I usually use it for respiratory belt detachment.

It might be that the scan triggers for the custom format are not properly integrated. You don't need them here anyway, if you specify the relative onset of first scan volume and physio logging as relative_start_acquisition field.

it might be beneficial to use the actual triggers recorded in the physio data, some misalignments are I guess possible and actual triggers are probably more reliable than assuming exact TR after the first trigger (I think this is why FSL PNM toolbox requires, per their explanations in the user guide)

This is true in general. For our Philips scanner here, I haven't seen much of a relevant drift in TR, but this might differ between systems. I will see to fix that for your data so you can use the triggers.

All the best,
Lars

from tapas.

Thanks for quick responses and suggestions how to deal with the data quality. We did try to put it PPU on the hand they didn't use to make responses, but they often used both hands. I will know more whether dropouts are movement related when I examine the physio data of more subjects, and I will definitely have to check the F-contrasts.

A detail regarding the triggers: we were manually stopping each run so we actually most of the time have 1 trigger more than the actual number of volumes (your scan_timing.sqpar.Nscans variable) - this is the reason why I was also specifiying "first" in log_files.align_scan - essentially Nscans should have priority over triggers in the files. Though, perhaps this is a specific aspect of my data and I should remove those extra triggers.

Perhaps I should open a separate issue for this, but I ran into an error when data in cardiac file would begin with heart beat (1, and then zeros). I artificially convert those cases into 0s, assuming it shouldnt be a big influence. Perhaps this caused an error only in my data as both columns in cardiac file were the same?

On a sidenote, is there a way to instruct the toolbox to not open windows for drawing figures, just saving figures instead? I'm doing it on a server but would like to save all the figures at the same time, so that I can later go through the figures to identify subjects/runs that need more scrutiny, without rerunning the code to generate the figures. At the moment I have to specify physio.verbose.level = 0; to run it on a server but then no figures are produced.

from tapas.

Dear Hrvoje,

• Nscans as priorities over triggers...I am not sure what you mean there exactly. I think the trigger time would be used from the trigger file, but the counting of scans should be according to Nscans, either starting from beginning or end of the log file.
• I will have a look at the error with the heart beat in the beginning...maybe you could indeed create a new issue and provide an example file?
• I think there are ways to create figures on a server without opening them, but it really depends on whether there is a proper graphics renderer on that machine. Instead, I would propose you try out tapas_physio_review() which recreates the most relevant figures post-hoc from any saved physio.mat structure (physio.model.output_physio = 'physio.mat';). The advantage there is that you can interactively browse and zoom in the figures, unlike jpegs etc.

All the best,
Lars

from tapas.

Nscans as priorities over triggers...I am not sure what you mean there exactly. I think the trigger time would be used from the trigger file, but the counting of scans should be according to Nscans, either starting from beginning or end of the log file.

This was referring to possibility to use the actual triggers from the file (as in these BIDS formatted files). But yes, you seemed to implement it that way anyway - trigger times from the trigger files and actual count according to the Nscans.

I will have a look at the error with the heart beat in the beginning...maybe you could indeed create a new issue and provide an example file?

Ok, will do.

from tapas.

Hello,

I have a quick additional question. I'm trying to verify whether physio regressors produce sensible effects with an F-contrast, per suggestion in your paper. I'm not using the functions provided in the toolbox as I'm using SPM from Nipype, so I would like to know a bit more details about it..

What is the GLM that you would run - do you include only the physio regressors in the 1st level GLM, and specify F contrast over all of the RETROICOR, HRV, RVT variables, perhaps also separate F over cardiac RETROICOR, respiratory, interactions etc? Or is it OK to include these contrasts in a GLM together with events and other nuisance regressors (motion etc)?

I actually already tried the second approach, and I get results that do not look like those reported in your paper - I get very little activation in the midbrain and around ventricles, so now I'm trying to figure out whether the GLM is wrong or I have a bug somewhere in my code producing the physio regressors.

Finally, is there a way to bring those F contrast to the group level, in a 2nd level GLM?

Thanks!

from tapas.

Dear Hrvoje,

I created a new issue to answer this question to keep the discussion thread clean.

from tapas.

Thanks, that's ok of course.

I see that you closed this issue, do you have a new version of the toolbox with the BIDS format support that I can try out?

from tapas.

Dear Hrvoje,

I updated the TAPAS development branch with the BIDS reader for PhysIO now. Have a look at

https://github.com/translationalneuromodeling/tapas/blob/development/PhysIO

in particular

https://github.com/translationalneuromodeling/tapas/blob/development/PhysIO/code/readin/tapas_physio_read_physlogfiles_bids.m

and the example in

https://github.com/translationalneuromodeling/tapas/blob/development/PhysIO/examples/tapas_physio_example_bids_cpulse3t_spm_job.m

Have a nice weekend,
Lars

from tapas.

That's great to hear, thanks a lot for a speedy development!

Hrvoje

from tapas.

I tried the new code and I find an error:

Error using rmdir
No directories were removed.

Error in tapas_physio_read_physlogfiles_bids (line 195)
    rmdir(fileName, 's');

Error in tapas_physio_read_physlogfiles (line 62)
            tapas_physio_read_physlogfiles_bids(log_files,
            cardiac_modality, verbose);

Error in tapas_physio_main_create_regressors (line 120)
                verbose] = tapas_physio_read_physlogfiles(...

fileName in rmdir is not really a directory name and hence the failure. It's a quick fix with fileparts(fileName) instead, no?

from tapas.

Dear Hrvoje,

thank you for spotting this! This was indeed a bug, and when I checked it on my Windows machine, I also got another access problem, since I didn't properly close the unzipped .tsv file with fclose(). Both should be fixed by the last commit in the development branch.

Best,
Lars

from tapas.

Great, I made the correction locally and it worked fine.

I didn't check that part yet, but its great that you added automatic extraction of parameters from the json file!

from tapas.

Hi Lars,

I got new data, and this time I collected analog cardiac measurements. If I feed this into PhysIO, with BIDS option I get the following error:

TR = 1.450 +/- 0.010 s (Estimated mean +/- std time of repetition for one volume)
No cardiac R-peak (heartbeat) events provided
    maxscan (incl. dummies) = 420 
    tmin (1st scan start (1st dummy))=   0.00 s
    tmin (1st scan start (after dummies))=   0.00 s
    tmax = 608.90 s 
    mean TR =   1.45 s

Error using stem (line 43)
X and Y must be the same size.

Error in tapas_physio_plot_cropped_phys_to_acqwindow
(line 68)
    stem(cpulse, ampc*ones(length(cpulse),1), 'r--')
    ;

Error in tapas_physio_main_create_regressors (line
222)
                    tapas_physio_plot_cropped_phys_to_acqwindow(ons_secs,
                    sqpar);

Instead of option

physio.preproc.cardiac.initial_cpulse_select.method = 'load_from_logfile';

I tried using

physio.preproc.cardiac.initial_cpulse_select.auto_matched.min = 0.4;
physio.preproc.cardiac.initial_cpulse_select.auto_matched.file = 'initial_cpulse_kRpeakfile.mat';

as in tapas_physio_example_siemens_hcp_ppu3t_spm_job.m.

BIDS format allows analog cardiac responses, perhaps you set it up according to my data which had heart beat data instead?

Here are the data files:
sub-s080.zip

Hrvoje

from tapas.