BOLD movement scrubbing#

Why scrub BOLD data based on movement parameters?#

When a subject moves in the scanner during BOLD acquisition, a number of things happen that result in possible movement artifacts-induced changes in signal intensity that are not related to actual changes in BOLD signal due to underlying neural activity, but rather due to movement itself. Some consequences of movement can be corrected using movement correction—individual volumes within the time series are realigned so that each voxel/grayordinate represents the same part of the brain. Regressing movement parameters out of the signal can also reduce the extent of the artifactual, movement-driven change in the signal. Large movements, however, still lead to relatively large artifactual changes in the signal.

When running activation analyses, a large part of the artifactual signal will most probably be averaged out and the final estimate of signal magnitude due to task or event of interest won't be too far off. Movement artifacts will however lead to increased error estimate and reduced power to detect activation within a single session, and might also increase variability of response magnitudes and reduce power in second level analyses.

When running functional connectivity, any noise signal that is shared between voxels or regions of interest will induce positive or negative correlation. The larger the ratio between correlated signal changes due to movement and due to signal of interest, the larger will be the artifactual correlations observed in the final results.

In order to increase the power of activation analyses and reduce the artifactual correlations in functional connectivity analyses, one has to minimize the effect of subject movement on BOLD signal. The best way of achieving that is to have subjects be as still as possible during scanning. Not all is lost, however, if there are just a few large movements. In this case, it is possible to identify the volumes (frames) in the BOLD image that might be affected by movement artifacts and ignore them either in activation analyses or when computing functional connectivity.

What follows is a description of how to go about identifying, visualizing and ignoring problematic volumes using QuNex tools.

The strategy#

The strategy in movement scrubbing is the following:

Compute statistics that reflect movement and its artifactual properties
Identify problematic frames
Ignore or interpolate problematic flagged frames in analyses

The tools described use two sources of information that can be used to identify bad frames. The first one are movement correction parameters. Movement correction is a key part of basic preprocessing of BOLD data. Any package used for movement correction provides an option to report movement correction parameters that describe the amount of shift and rotation that was needed to align each volume of the BOLD data to the reference volume. Tracking the changes in these (most often 6) parameters provides an estimate of frame-to-frame movement of the subject. Having that information, it is possible to combine it into one parameter reflecting the overall amount of movement and use it to identify the frames where movement exceeded an arbitrary threshold. The parameter computed and used by the QuNex tools is frame displacement as described in Power et al. (2012). The default threshold value is 0.5.

The second source of information is the acquired BOLD timeseries itself. By comparing intensity values from frame to frame, one can track the signal changes across the BOLD timeseries and identify frames in which signal change exceeds a chosen threshold. The specific parameter used to measure signal change by QuNex tools is dvars, described by Power et al. (2012). In short, dvars is computed by calculating difference between two timeframes for each voxel within the brain separately. The differences are then squared and a square root of the mean is computed. To normalize it, the computed RMS (root mean square) is divided by mean intensity and multiplied by 100 to reflect mean percent change of the signal. Power et al. (2012) suggest to use 3 as a sensible threshold value for the resulting parameter.

Dvars reflects not only movement effects on the signal but also the effects of signal of interest and overall noise, and quite possibly other properties such as voxel size, TR, brain size, tissue composition and others. For that reason some BOLD timeseries have higher mean dvars then others, reflecting higher baseline frame-to-frame differences. When baseline dvars is higher, relatively smaller changes in the signal can lead to crossing the chosen threshold. To account for the differences in baseline frame-to-frame variability QuNex tools also express dvars at each frame as a ratio to baseline dvars, computed as median across the BOLD timeseries. The resulting dvarsme statistic for each frame is interpreted as the factor of change compared to baseline (e.g. dvarsme of 2 means that the change in that frame was 2 times higher than baseline change between successive frames across the run). The proposed default dvarsme threshold is 1.5.

Using both sources frames with excessive movement and/or signal change are easily identified. As the changes can be smeared in time due to filtering and other manipulations of the BOLD timeseries Power et al. (2012) suggest to reject not only the identified frame, but also one preceding and two following frames. QuNex tools allow specification of how many frames before and after identified artifactual frames to be also marked as bad.

To choose the most appropriate threshold and have a feel for your data, it is prudent to visually check the movement parameters and the extent of the suggested frames to be removed. QuNex tools provide means for both visual inspection as well as tabular review of the data. Do make use of them instead of just accepting the defaults.

Once the frames are identified the software for running the relevant analyses needs to be informed about which frames to ignore. The tools provide two venues for specifying which frames to remove or ignore in analyses. The first option is to record the relevant frames in standard .fidl file format. This option creates a .fidl file in which each line records the beginning of contiguous frames to be ignored and the number of frames to be ignored as a negative value. This is the format that fidl as well as QuNex scripts that work with event-related data will accept. To use the information specified in such generated .fidl files, they first need to be combined appropriately with existing .fidl files. This is also accomplished with the use of QuNex tools and existing .conc and .fidl files.

The second option in specifying the frames to be ignored is to generate a scrub.txt file. This is a regular plaintext file in which each frame is marked either with 0 (do not ignore) or 1 (ignore). The information is provided in separate columns for each of the discussed criteria (frame displacement, dvars and dvarsme) as well as their intersection (a frame is considered bad if both criteria are met) or union (a frame is considered bad if either of the criteria is reached).

In the following sections tools for performing each of the listed steps will be presented in more detail. First, the two methods used to compute the statistics and scrubbing will be presented, followed by the functions and commands that use them to process image files in bulk.

Supported methods implemented in MATLAB / Octave#

For computation of BOLD statistics, identification of bad frame img_stats_time and img_compute_scrub methods are used, respectively. Note that these are 'core' functions implemented in MATLAB (or Octave). These functions are called automatically by 'higher-order' QuNex BOLD de-noising tools. However, a user can envoke them specifically within MATLAB.

img_stats_time method#

img_stats_time is a method that computes a desired list of statistics across all specified voxels for each frame of the BOLD image. To use the method run:

[stats, done] = img.img_stats_time(do, mask);

Both parameters are optional. Do should be a cell array of strings that specify what statistics to compute and mask the object that defines a mask for voxel selection. Statistics are computed only across voxels that have a non-zero value in the mask image. If do is absent or empty, all the statistics will be computed by default. If no mask is provided, statistics are computed across all voxels in a volume. The statistics to be computed are specified with the following strings:

n — number of non-NaN values in the mask for each frame
m — mean across all voxels in the mask for each frame
me — median across all voxels in the mask for each frame
max — maximum value within all voxels in the mask for each frame
min — minimum value within all voxels in the mask for each frame
sum — a sum of values of voxels in the mask for each frame
sd — standard deviation computed across all voxels in the mask for each frame
var — variance computed across all voxels in the mask for each frame
dvars — signal change across all voxels in the mask between each frame and the preceding frame

The results are stored in a structure stats that consists of one-dimensional vectors in named fields. The names of the fields are the same as the commands to compute them. The only exception is computation of dvars, which computes two additional statistics: dvarsm is computed by dividing the value of dvars for each frame by the mean for that frame and multiplying by 100; dvarsme is computed by dividing the value of dvarsm for each frame by the median dvarsm computed across all frames. For explanation see the text above. A call:

stats = img.img_stats_time({'m','var','dvars'}, mask);

will return the following fields:

   stats.m
   stats.var
   stats.dvars
   stats.dvarsm
   stats.dvarsme

Besides the structure with the statistics, the method also returns a cell array done that lists the executed commands.

img_compute_scrub method#

img_compute_scrub is a method that evaluates BOLD motion and image statistics to identify bad frames. It is used by calling:

img = img_compute_scrub(img, comm);

img is the image to be evaluated, and comm is the command description of how to compute scrubbing. The command is provided as a string of 'key:value' pairs separated by a pipe (|) character. The following keys can be specified, with the defaults in []:

radius ... head radius in mm [50]
fdt ... frame displacement threshold [0.5]
dvarsmt ... image intensity normalized root mean squared error (RMSE) (dvarsm) threshold [3.0]
dvarsmet ... median normalized RMSE (dvarsme) threshold [1.6]
after ... how many frames after the identified bad ones are to be marked as well [0]
before ... how many frames before the identified bad ones are to be marked as well [0]
reject ... which criteria is to be used for identification of bad frames [udvarsme]

The possible criteria for identification of bad frames are:

mov ... frame displacement threshold (fdt) is exceeded
dvars ... image intensity normalized root mean squared error (RMSE) threshold (dvarsmt) is exceeded
dvarsme ... median normalized RMSE (dvarsmet) threshold is exceeded
idvars ... both fdt and dvarsmt are exceeded (i for intersection)
uvars ... either fdt or dvarsmt are exceeded (u for union)
idvarsme ... both fdt and dvarsmet are exceeded
udvarsme ... either fdt or udvarsmet are exceeded

The method is used on an mrimage img object that has both movement parameters data (img.mov field structure) and per frame statistics data (img.fstats structure) already computed and present. Movement data is read automatically when the image is loaded if the data is present in the relevant folder (images/functional/movement). Frame statistics data is also read if present. It is computed using img_stats_time method and is added and saved when using general_compute_bold_stats command (see below).

Note, the image is not scrubbed, only the scrubbing parameters are computed and the chosen criterium is used to create the img.use mask that in other functions and methods defines what frames to use or not.

Processing BOLD images#

QuNex tools provide multiple ways of processing BOLD for image statistics and scrubbing. One possibility is to use the methods described above directly. The other method is to use QuNex MATLAB and QuNex python functions to process single BOLD image files (general_compute_bold_stats QuNex MATLAB command), a list of BOLD image files (general_compute_bold_list_stats QuNex MATLAB command) or BOLD image files across sessions (qunex compute_bold_stats terminal command).

general_compute_bold_stats function#

When you want to perform complete processing of a BOLD image file, i.e. reading in the relevant data, computing statistics and scrubbing, and saving the results, general_compute_bold_stats is the command to use. The function is run as:

general_compute_bold_stats(img, mask, target, store, scrub, verbose);

with the following parameters:

img ... A path to a BOLD file to process.
mask ... A path to a mask file to use.
target ... A folder to save results into.

If target is set to '' (the default), the results are saved where the BOLD image is, if target is specified as 'none', results are not saved to an external file.
store ... Whether to store the resulting data in the image file.

If store is set to 'same', the computed statistics and scrubbing information are embedded and saved in the original BOLD image itself. If store is set to <ext>, the information is stored embedded in a copy of the original BOLD image, with extension <ext>.

If store is left blank (i.e. '', the default), the results are not embedded and saved in an image file.

NOTE: Due to changes of the default format, this option is currently disabled and might be removed.
scrub ... A string describing whether and how to compute scrubbing information.

See img_compute_scrub above. If 'none' is specified, no scrubbing computation is to be performed.
verbose ... Whether to report the progress or not (the default).

The command initially identifies relevant brain voxels in two manners. First, it identifies voxels with intensity higher than 300 on the first BOLD frame. If there are more than 20000 valid voxels, it then selects those for which the intensity is above the specified threshold throughout the timeseries and selects those for computation of image statistics. Second, if the first method fails (e.g. in the case when images were demeaned), it identifies all the voxels for which the variance across the frames is more than 0. After the voxels are identified, the image is additionally masked if a mask was specified.

To compute the per frame statistics on the identified and masked brain voxels, the command runs img_stats_time method and saves the computed statistics in a simple text file named the same as the BOLD file it relates to with .bstats extension. The structure of the results is rather straightforward. The file starts with a header line that lists all the variables, whereas all the following lines list statistics, one frame per line. All the values are separated by tabs. The columns are provided in the following order: frame (frame number), n, m, min, max, var, sd, dvars, dvarsm, dvarsme.

If scrubbing was specified it also runs img_compute_scrub method and creates a .scrub file with detailed scrubbing information, as well as .use file with information on which frames are to be used, and which not. Both are simple text files. .scrub file starts with the header line that specifies the columns: frame, mov, dvars, dvarsme, idvars, idvarsme, udvars, udvarsme, use. The following lines provide the information on frame number, whether the frame is marked as bad (1) or not (0) by each specified criterium, and whether the frame is to be used (1) or not (0) based on the criterium specified when running the function.

general_compute_bold_list_stats function#

When you want to compute BOLD statistics for a larger number of images using MATLAB, general_compute_bold_list_stats is the right command to use. To use the command execute:

general_compute_bold_list_stats(flist, target, store, scrub, verbose);

with the following parameters:

flist ... .list text file providing a list of session's BOLD image or conc files to process.

See list files for more information.
target ... A folder (relative to where BOLD image is stored) to save results into.

If target is set to '' (the default), the results are saved where the BOLD image is,

if target is specified as 'none', results are not saved to an external file.
store ... Whether to store the resulting data in the image file.

If store is set to 'same', the computed statistics and scrubbing information are embedded and saved in the original BOLD image itself. If store is set to <ext>, the information is stored embedded in a copy of the original BOLD image, with extension <ext>.

If store is left blank (i.e. '', the default), the results are not embedded and saved in an image file.

NOTE: Due to changes of the default format, this option is currently disabled and might be removed.
scrub ... A string describing whether and how to compute scrubbing information.

See img_compute_scrub above. If 'none' is specified, no scrubbing computation is to be performed.
verbose ... Whether to report the progress or not (the default).

The command runs general_compute_bold_stats on each BOLD file specified in the provided .list file. The only parameter that has to be provided is filelist, whereas target and verbose are optional. If they are not provided, the general_compute_bold_stats defaults will be used. If filelist specifies an ROI file for a session (see list file format, that ROI will be used as a mask for computing statistics. If an ROI file is not specified, only the general_compute_bold_stats internal method for identifying brain signal will be used.

QuNex compute_bold_stats command#

Note that the compute_bold_stats is a 'higher-level QuNex command that wraps the 'core' functions described above. Thus, the easiest way to compute BOLD image statistics for all the sessions in a study is to use the compute_bold_stats QuNex command. This command will process the specified BOLD image files for all the sessions in the specified batch file, taking into account parameters for filtering or explicitly specifying which sessions and BOLD files to work with. Further, all the specific scrubbing parameters described above can be specified. For detailed information about the use of the command, and the results, please consult online help by running qunex compute_bold_stats. In brief, this is a description of key parameters:

general parameters

When running the command, the following general processing parameters are taken into account:

--batchfile       ... The batch.txt file with all the sessions information [batch.txt].
--studyfolder     ... The path to the study/sessions folder, where the imaging data is supposed to go [.].
--parsessions     ... How many sessions to run in parallel [1].
--overwrite       ... Whether to overwrite existing data (yes) or not (no) [no].
--bold_preprocess ... Which BOLD images (as they are specified in the batch.txt file) to copy over. 
                      It can be a single type (e.g. 'task'), a pipe separated list (e.g. 
                      'WM|Control|rest') or 'all' to work on all [rest].
--boldname        ... The default name of the bold files in the images folder [bold].

specific parameters

In addition the following specific parameters define the actual results:

--mov_radius  ... Estimated head radius (in mm) for computing frame displacement statistics [50].
--mov_fd      ... Frame displacement threshold (in mm) to use for identifying bad frames [0.5]
--mov_dvars   ... The (mean normalized) dvars threshold to use for identifying bad frames [3.0].
--mov_dvarsme ... The (median normalized) dvarsm threshold to use for identifying bad frames [1.5].
--mov_after   ... How many frames after each frame identified as bad to also exclude from further 
                  processing and analysis [0].
--mov_before  ... How many frames before each frame identified as bad to also exclude from further 
                  processing and analysis [0].
--mov_bad     ... Which criteria to use for identification of bad frames [udvarsme].

This command is best used within the general preprocessing workflow of BOLD data for purposes of functional connectivity preprocessing or preparation for task activation analysis.

Examples#

Using the defaults:

qunex compute_bold_stats \
    --batchfile="<path to batch file>" \
    --studyfolder="<path to study sessions folder>" \
    --overwrite="<whether to overwrite existing files (yes) or not (no)>" \
    --bold_preprocess="<which BOLD files to process, e.g. 'all' or 'rest'>"

Specifying additional parameters for identification of bad frames:

qunex compute_bold_stats \
    --batchfile="<path to batch file>" \
    --studyfolder="<path to study sessions folder>" \
    --overwrite="<whether to overwrite existing files (yes) or not (no)>" \
    --bold_preprocess="<which BOLD files to process, e.g. 'all' or 'rest'>" \
    --mov_fd="0.9" \
    --mov_dvarsme="1.6" \
    --mov_before="1" \
    --mov_after="2"

The breakdown of the parameters is as follows:

mov_fd="0.9" ... sets frame displacement threshold to 0.9
mov_dvarsme="1.6" ... sets dvarsme threshold to 1.6
mov_before="1" ... exclude also one frame before identified bad ones
mov_after="2" ... exclude also one frame after identified bad ones

g_BoldMovement R script#

g_BoldMovement.R is a R script that works on .bstat files and generates .scrub and ._scrub.fidl files to be used in further preprocessing, as well as graphical representation of the statistics and resulting scrubbing. To create graphs the script depends on ggplot2 library that should be installed on the local R installation.

The simplest way to run the script is to execute g_BoldMovement from the command line. In this case the script will run on the .dat and .bstat files that it finds in the folder in which it was ran, with the default values and settings. For more advanced use, a number of parameters and options can be listed when running the script. The general format for specifying parameters is:

g_BoldMovement.R --args -option=value -option=value -option=value

The available options and parameters, and their default values (in brackets) are:

folder or f — the folder in which to look for .dat and .bstat files [current working folder]
mreport or mr — the file to add movement report to [none]
sreport or sr — the file to add scrubbing report to [none]
subject or s — the subject id to use in plots and reports [none]
dvars or d — the threshold to use for dvars based rejection [3]
dvarsme or e — the threshold to use for dvarsme based rejection [1.5]
movement or m — the threshold to use for frame-to-frame movement rejection [0.5]
radius or rd — radius (in mm) from center of head to cortex to estimate rotation based movement [50]
tr — TR to be used when generating _scrub.fidl files [2.5]
scrub — whether to output scrub file (yes) or not (no) [no]
fidl — whether to output and what to base _scrub.fidl files on, possible values are:
- fd — use frame-to-frame displacement only
- dvars — use dvars only
- dvarsme — use dvarsme only
- u — use union of frame-to-frame displacement and dvars
- ume — use union of frame-to-frame displacement and dvarsme
- i — use intersection of frame-to-frame displacement and dvars
- ime — use intersection of frame-to-frame displacement and dvarsme
- none — do not output a _scrub.fidl file [the default]
plot — whether to create plots (yes) or not (no) [no]

The following example:

g_BoldMovement.R --args -f=/data/study1/s01/movement -mr=/data/study1/mov_report.txt -tr=3 -fidl=ume -plot=yes

would look for data in folder /data/study1/s01/movement, identify bad frames as those in which either the frame-to-frame displacement or dvarsme exceeded the default threshold, save the to-be-ignored frames for each processed run in file ending with _scrub.fidl, and also create and save graphical reports. A short overview of movement parameters across all frames for each BOLD run will be appended to file /data/study1/mov_report.txt.

When specified, three graphical report files are saved:

mov-cor-report.pdf contains plots of 6 movement corrections for each run separately,
move-dvars-report.pdf contains plots of frame-to-frame displacement, dvars and to-be-ignored frames based on them for each run separately,
move-dvarsme-report.pdf contains plots of frame-to-frame displacement, dvarsme and to-be-ignored frames based on them for each run separately.

QuNex create_stats_report#

To ease running g_BoldMovement, QuNex provides create_stats_report command that runs g_BoldMovement on all the specified sessions and BOLD images in the study. Additionally, the command will ease generation of group report files providing an easy way to review group level statistics related to movement scrubbing. The same as compute_bold_stats commands, create_stats_report accepts both general and specific parameters. For detailed information, run qunex ?create_stats_report. In brief, this is a list of key parameters.

general parameters

When running the command, the following general processing parameters are taken into account:

--batchfile       ... The batch.txt file with all the session information [batch.txt].
--studyfolder     ... The path to the study/sessions folder, where the imaging data is supposed to go [.].
--parsessions     ... How many sessions to run in parallel [1].
--overwrite       ... Whether to overwrite existing data (yes) or not (no) [no].
--bold_preprocess ... Which BOLD images (as they are specified in the batch.txt file) to copy over. 
                      It can be a single type (e.g. 'task'), a pipe separated list (e.g. 
                      'WM|Control|rest') or 'all' to work on all [rest].
--boldname        ... The default name of the bold files in the images folder [bold].

scrubbing specific options

In addition the following scrubbing specific parameters define computation of results:

--mov_radius  ... Estimated head radius (in mm) for computing frame displacement statistics [50].
--mov_fd      ... Frame displacement threshold (in mm) to use for identifying bad frames [0.5]
--mov_dvars   ... The (mean normalized) dvars threshold to use for identifying bad frames [3.0].
--mov_dvarsme ... The (median normalized) dvarsm threshold to use for identifying bad frames [1.5].
--mov_after   ... How many frames after each frame identified as bad to also exclude from further 
                  processing and analysis [0].
--mov_before  ... How many frames before each frame identified as bad to also exclude from further 
                  processing and analysis [0].
--mov_bad     ... Which criteria to use for identification of bad frames [udvarsme].

reportng specific options

Finally, the following options determine how the results will be reported and saved:

--TR          ... TR of the BOLD files [2.5].
--mov_pref    ... The prefix to be used for the group reports [].
--mov_plot    ... The base name of the plot files. If set to empty no plots are generated [mov_report].
--mov_mreport ... The name of the group movement report file. If set to an empty string, no file is 
                  generated [movement_report.txt].
--mov_sreport ... The name of the group scrubbing report file. If set to an empty string, no file is 
                  generated [movement_scrubbing_report.txt].
--mov_preport ... The name of group report file with stats computed with frames identified as bad 
                  exluded from analysis. If set to an empty string, no file is generated
                  [movement_report_post.txt].
--mov_post    ... The criterium for identification of bad frames that is used when generating a post 
                  scrubbing statistics group report (fd/dvars/dvars/dvarsme/idvars/idvarsme/udvars/
                  udvarsme/none) [udvarsme].
--mov_fidl    ... Whether to create fidl file snippets with listed bad frames, and what criterium to 
                  use for the definition of bad frames (fd/dvars/dvars/dvarsme/idvars/idvarsme/udvars/
                  udvarsme). Set to none to not generate them [udvarsme].
--mov_pdf     ... The name of the folder in sessions/QC/movement in which to copy the individuals' 
                  movement plots [movement_plots].

Again this command is best used within the general preprocessing workflow of BOLD data for purposes of functional connectivity preprocessing or preparation for task activation analysis.

Example#

qunex create_stats_report \
    --batchfile="<path to batch file>" \
    --sessionsfolder="<path to study sessions folder>" \
    --overwrite="<whether to overwrite existing reports and files (yes) or no (no)>" \
    --bold_preprocess="<which BOLD files to process, e.g. 'all' or 'rest'>" \
    --parsessions="<how many sessions to run in parallel, e.g. 1>"

Additional supporting functions#

QuNex also provides additional supporting functions.

meltmovfidl#

To be able to use the generated _scrub.fidl files they possibly need to be combined and melted with existing .fidl files. Combination ensures that when working with multiple runs, as specified by .conc file, _scrub.fidl that are generated for each run are combined in the correct order and their timing information is adjusted to correspond to the right frame within the concatenated run. Melting enables scrubbing information to be combined with existing event information in a fidl file. Melting is accomplished using:

qunex meltmovfidl <concfile> <ignore_fidl_pattern> <input_fidl_file> <output_fidl_file>

All the listed parameters have to be provided. --concfile provides a path to the .conc file that describes the sequence of BOLD runs to be concatenated in analyses. --ignore_fidl_pattern specifies the ending for the fidl files to be melted with the original fidl file. Typically the ending would be scrub.fidl. --input_fidl_file lists the path to the original fidl file that the scrubbing information is to be added to. --output_fidl_file gives the path to the file in which the melted data is to be saved.

In the following example:

qunex /data/study1/s01/wm.conc _scrub.fidl /data/study1/s01/wm.fidl /data/study1/s01/wm_scrubbed.fidl

The command will generate a new fidl file, based on /data/study1/s01/wm.fidl and /data/study1/s01/wm.conc, that will include scrubbing information in *_scrub.fidl files and save it to /data/study1/s01/wm_scrubbed.fidl file.

Warning on file names and paths

It needs to be taken into account that to correctly match BOLD images listed in the conc file with the corresponding per-run fidl files, both have to include in their name a sequence b*N where N corresponds to a BOLD number. Examples would be s01_bold2.4dfp.img, ad7325_b04_x3rd_atl.nii, id285_brun3_scrub.fidl. In addition, per-run fidl files to be used are expected to reside in the movement folder that corresponds to the session from which the BOLD run is taken from.

QuNex documentation

BOLD movement scrubbing

Contents