Open In Binder Open In Colab

FSL’s fMRI Data Analysis via Nipype#

Author: Monika Doerig

Original peper: https://www.frontiersin.org/articles/10.3389/fnimg.2022.953215/full

Original code: https://osf.io/prg53/?view_only=9d7974834a484cdb972bcc3989589b78

Setup Neurodesk#

%%capture
import os
import sys
IN_COLAB = 'google.colab' in sys.modules

if IN_COLAB:
  os.environ["LD_PRELOAD"] = "";
  os.environ["APPTAINER_BINDPATH"] = "/content,/tmp,/cvmfs"
  os.environ["MPLCONFIGDIR"] = "/content/matplotlib-mpldir"
  os.environ["LMOD_CMD"] = "/usr/share/lmod/lmod/libexec/lmod"

  !curl -J -O https://raw.githubusercontent.com/NeuroDesk/neurocommand/main/googlecolab_setup.sh
  !chmod +x googlecolab_setup.sh
  !./googlecolab_setup.sh

  os.environ["MODULEPATH"] = ':'.join(map(str, list(map(lambda x: os.path.join(os.path.abspath('/cvmfs/neurodesk.ardc.edu.au/neurodesk-modules/'), x),os.listdir('/cvmfs/neurodesk.ardc.edu.au/neurodesk-modules/')))))
# Output CPU information:
!cat /proc/cpuinfo | grep 'vendor' | uniq
!cat /proc/cpuinfo | grep 'model name' | uniq
vendor_id	: AuthenticAMD
model name	: AMD EPYC 7742 64-Core Processor
!pip install pandas nilearn matplotlib

First-level GLM using Nipype FSL#

In this notebook, we recreate the first-level GLM and the first level GLM of FSL GUI using nipype code. For each nipype node, we list the corresponding fsl command from the log file. The dataset we use is a Flanker task, which can be downloaded here.

We also borrow some helps from this document.

#load fsl module
import lmod
await lmod.purge(force=True)
await lmod.load('fsl/6.0.5.1') #Original pipeline: FSL 6.0.4, nipype 1.6.1.
await lmod.list()
['fsl/6.0.5.1']
import os
os.environ["FSLDIR"]="/cvmfs/neurodesk.ardc.edu.au/containers/fsl_6.0.5.1_20221016/fsl_6.0.5.1_20221016.simg/opt/fsl-6.0.5.1/"
os.environ["FSLOUTPUTTYPE"]="NIFTI_GZ"
os.environ["SINGULARITY_BINDPATH"]="/data,/neurodesktop-storage,/tmp,/cvmfs"

Preparation#

Import all the relevant libraries needed for the preprocessing stage.

from __future__ import print_function
from __future__ import division
from builtins import str
from builtins import range

from glob import glob

from nipype import Function
from nipype.interfaces import fsl, utility as util, io as nio
import nipype.pipeline.engine as pe  # pypeline engine
import nipype.algorithms.modelgen as model  # model generation
import nipype.algorithms.rapidart as ra  # artifact detection
import matplotlib.pyplot as plt

import nilearn
from nilearn import surface
from nilearn import plotting
# check output type and fsl version
print(fsl.Info.output_type())
print(fsl.Info.version())
NIFTI_GZ
6.0.5.1:57b01774

Set up data path#

# Input: Set the BIDS path
data_dir_first = os.path.join(os.getcwd(), 'ds000102')
# Output: Set path where nipype will store stepwise results
exp_dir_first = os.path.join(os.getcwd(), 'output_level1')

# Input: Set path of first level outputs
data_dir_second = os.path.join(os.getcwd(), 'output_level1/level1_results/')
# Output: Set path where nipype will store stepwise results
exp_dir_second = os.path.join(os.getcwd(), 'output_level2/')

# Input: Set path of second level outputs
data_dir_third = os.path.join(os.getcwd(), 'output_level2/level2_results/')
# Output: Set path where nipype will store stepwise results
exp_dir_third = os.path.join(os.getcwd(), 'output_level3/')

PATTERN = "sub-0*"
!datalad install https://github.com/OpenNeuroDatasets/ds000102.git
!cd ds000102 && datalad get $PATTERN
Clone attempt:   0%|              | 0.00/2.00 [00:00<?, ? Candidate locations/s]
Enumerating: 0.00 Objects [00:00, ? Objects/s]
                                              
Counting:   0%|                               | 0.00/27.0 [00:00<?, ? Objects/s]
                                                                                
Compressing:   0%|                            | 0.00/23.0 [00:00<?, ? Objects/s]
                                                                                
Receiving:   0%|                             | 0.00/2.15k [00:00<?, ? Objects/s]
                                                                                
Resolving:   0%|                                | 0.00/537 [00:00<?, ? Deltas/s]
[INFO   ] Filesystem allows writing to files whose write bit is not set.        
[INFO   ] Detected a crippled filesystem. 
[INFO   ] Disabling core.symlinks. 
[INFO   ] Entering an adjusted branch where files are unlocked as this filesystem does not support locked files. 
[INFO   ] Switched to branch 'adjusted/master(unlocked)' 
[INFO   ] Remote origin not usable by git-annex; setting annex-ignore 
[INFO   ] https://github.com/OpenNeuroDatasets/ds000102.git/config download failed: Not Found 
[INFO   ] access to 1 dataset sibling s3-PRIVATE not auto-enabled, enable with:
| 		datalad siblings -d "/home/jovyan/example-notebooks/functional_imaging/ds000102" enable -s s3-PRIVATE 
install(ok): /home/jovyan/example-notebooks/functional_imaging/ds000102 (dataset)
Total:   0%|                                    | 0.00/618M [00:00<?, ? Bytes/s]
Get sub-07/a .. 7_T1w.nii.gz:   0%|            | 0.00/10.7M [00:00<?, ? Bytes/s]
Get sub-07/a .. 7_T1w.nii.gz:   0%|    | 18.9k/10.7M [00:00<01:28, 121k Bytes/s]
Get sub-07/a .. 7_T1w.nii.gz:   0%|    | 50.7k/10.7M [00:00<01:07, 158k Bytes/s]
Get sub-07/a .. 7_T1w.nii.gz:   1%|     | 138k/10.7M [00:00<00:34, 304k Bytes/s]
Get sub-07/a .. 7_T1w.nii.gz:   3%|▏    | 277k/10.7M [00:00<00:20, 515k Bytes/s]
Get sub-07/a .. 7_T1w.nii.gz:   5%|▎    | 573k/10.7M [00:00<00:10, 947k Bytes/s]
Get sub-07/a .. 7_T1w.nii.gz:  11%|▎  | 1.16M/10.7M [00:01<00:05, 1.79M Bytes/s]
Get sub-07/a .. 7_T1w.nii.gz:  22%|▋  | 2.36M/10.7M [00:01<00:02, 3.45M Bytes/s]
Get sub-07/a .. 7_T1w.nii.gz:  35%|█  | 3.73M/10.7M [00:01<00:01, 4.69M Bytes/s]
Get sub-07/a .. 7_T1w.nii.gz:  54%|█▌ | 5.76M/10.7M [00:01<00:00, 7.82M Bytes/s]
Get sub-07/a .. 7_T1w.nii.gz:  73%|██▏| 7.78M/10.7M [00:01<00:00, 10.6M Bytes/s]
Get sub-07/a .. 7_T1w.nii.gz:  88%|██▋| 9.39M/10.7M [00:01<00:00, 10.6M Bytes/s]
Total:   2%|▍                          | 10.7M/618M [00:03<02:57, 3.42M Bytes/s]
Get sub-07/f .. _bold.nii.gz:   0%|            | 0.00/28.9M [00:00<?, ? Bytes/s]
Get sub-07/f .. _bold.nii.gz:   6%|▏  | 1.88M/28.9M [00:00<00:02, 10.8M Bytes/s]
Get sub-07/f .. _bold.nii.gz:  13%|▍  | 3.80M/28.9M [00:00<00:02, 10.9M Bytes/s]
Get sub-07/f .. _bold.nii.gz:  20%|▌  | 5.75M/28.9M [00:00<00:02, 11.0M Bytes/s]
Get sub-07/f .. _bold.nii.gz:  26%|▊  | 7.52M/28.9M [00:00<00:01, 12.8M Bytes/s]
Get sub-07/f .. _bold.nii.gz:  33%|█  | 9.70M/28.9M [00:00<00:01, 10.8M Bytes/s]
Get sub-07/f .. _bold.nii.gz:  40%|█▏ | 11.7M/28.9M [00:01<00:01, 11.0M Bytes/s]
Get sub-07/f .. _bold.nii.gz:  47%|█▍ | 13.7M/28.9M [00:01<00:01, 11.2M Bytes/s]
Get sub-07/f .. _bold.nii.gz:  54%|█▋ | 15.8M/28.9M [00:01<00:01, 11.4M Bytes/s]
Get sub-07/f .. _bold.nii.gz:  62%|█▊ | 17.8M/28.9M [00:01<00:00, 13.2M Bytes/s]
Get sub-07/f .. _bold.nii.gz:  69%|██ | 19.9M/28.9M [00:01<00:00, 11.4M Bytes/s]
Get sub-07/f .. _bold.nii.gz:  75%|██▏| 21.6M/28.9M [00:01<00:00, 12.6M Bytes/s]
Get sub-07/f .. _bold.nii.gz:  82%|██▍| 23.8M/28.9M [00:02<00:00, 11.9M Bytes/s]
Get sub-07/f .. _bold.nii.gz:  90%|██▋| 26.2M/28.9M [00:02<00:00, 12.0M Bytes/s]
Get sub-07/f .. _bold.nii.gz:  97%|██▉| 28.2M/28.9M [00:02<00:00, 12.6M Bytes/s]
Total:   6%|█▋                         | 39.7M/618M [00:05<01:27, 6.65M Bytes/s]
Get sub-07/f .. _bold.nii.gz:   0%|            | 0.00/29.0M [00:00<?, ? Bytes/s]
Get sub-07/f .. _bold.nii.gz:   9%|▎  | 2.60M/29.0M [00:00<00:02, 13.0M Bytes/s]
Get sub-07/f .. _bold.nii.gz:  14%|▍  | 4.10M/29.0M [00:00<00:02, 11.1M Bytes/s]
Get sub-07/f .. _bold.nii.gz:  21%|▋  | 6.22M/29.0M [00:00<00:02, 9.35M Bytes/s]
Get sub-07/f .. _bold.nii.gz:  27%|▊  | 7.89M/29.0M [00:00<00:01, 11.1M Bytes/s]
Get sub-07/f .. _bold.nii.gz:  36%|█  | 10.4M/29.0M [00:01<00:01, 9.39M Bytes/s]
Get sub-07/f .. _bold.nii.gz:  41%|█▏ | 12.0M/29.0M [00:01<00:01, 9.25M Bytes/s]
Get sub-07/f .. _bold.nii.gz:  47%|█▍ | 13.6M/29.0M [00:01<00:01, 9.21M Bytes/s]
Get sub-07/f .. _bold.nii.gz:  52%|█▌ | 15.2M/29.0M [00:01<00:01, 9.19M Bytes/s]
Get sub-07/f .. _bold.nii.gz:  58%|█▋ | 16.8M/29.0M [00:01<00:01, 9.30M Bytes/s]
Get sub-07/f .. _bold.nii.gz:  64%|█▉ | 18.4M/29.0M [00:01<00:01, 9.24M Bytes/s]
Get sub-07/f .. _bold.nii.gz:  69%|██ | 20.1M/29.0M [00:02<00:00, 9.31M Bytes/s]
Get sub-07/f .. _bold.nii.gz:  75%|██▎| 21.7M/29.0M [00:02<00:00, 9.84M Bytes/s]
Get sub-07/f .. _bold.nii.gz:  80%|██▍| 23.2M/29.0M [00:02<00:00, 10.8M Bytes/s]
Get sub-07/f .. _bold.nii.gz:  87%|██▌| 25.1M/29.0M [00:02<00:00, 9.07M Bytes/s]
Get sub-07/f .. _bold.nii.gz:  92%|██▊| 26.8M/29.0M [00:02<00:00, 10.1M Bytes/s]
Get sub-07/f .. _bold.nii.gz:  98%|██▉| 28.5M/29.0M [00:02<00:00, 9.17M Bytes/s]
Total:  11%|██▉                        | 68.6M/618M [00:09<01:14, 7.36M Bytes/s]
Get sub-02/a .. 2_T1w.nii.gz:   0%|            | 0.00/10.7M [00:00<?, ? Bytes/s]
Get sub-02/a .. 2_T1w.nii.gz:  16%|▍  | 1.73M/10.7M [00:00<00:00, 9.87M Bytes/s]
Get sub-02/a .. 2_T1w.nii.gz:  32%|▉  | 3.46M/10.7M [00:00<00:00, 9.91M Bytes/s]
Get sub-02/a .. 2_T1w.nii.gz:  48%|█▍ | 5.21M/10.7M [00:00<00:00, 9.94M Bytes/s]
Get sub-02/a .. 2_T1w.nii.gz:  65%|█▉ | 6.96M/10.7M [00:00<00:00, 9.98M Bytes/s]
Get sub-02/a .. 2_T1w.nii.gz:  81%|██▍| 8.71M/10.7M [00:00<00:00, 10.0M Bytes/s]
Get sub-02/a .. 2_T1w.nii.gz:  98%|██▉| 10.5M/10.7M [00:01<00:00, 10.1M Bytes/s]
                                                                                
Get sub-02/f .. _bold.nii.gz:   0%|            | 0.00/29.2M [00:00<?, ? Bytes/s]
Get sub-02/f .. _bold.nii.gz:   6%|▏  | 1.78M/29.2M [00:00<00:02, 10.2M Bytes/s]
Get sub-02/f .. _bold.nii.gz:  12%|▎  | 3.56M/29.2M [00:00<00:02, 10.2M Bytes/s]
Get sub-02/f .. _bold.nii.gz:  18%|▌  | 5.34M/29.2M [00:00<00:02, 10.2M Bytes/s]
Get sub-02/f .. _bold.nii.gz:  24%|▋  | 7.14M/29.2M [00:00<00:02, 10.2M Bytes/s]
Get sub-02/f .. _bold.nii.gz:  31%|▉  | 8.95M/29.2M [00:00<00:01, 10.3M Bytes/s]
Get sub-02/f .. _bold.nii.gz:  37%|█  | 10.7M/29.2M [00:01<00:01, 10.3M Bytes/s]
Get sub-02/f .. _bold.nii.gz:  43%|█▎ | 12.5M/29.2M [00:01<00:01, 10.3M Bytes/s]
Get sub-02/f .. _bold.nii.gz:  49%|█▍ | 14.4M/29.2M [00:01<00:01, 10.3M Bytes/s]
Get sub-02/f .. _bold.nii.gz:  55%|█▋ | 16.0M/29.2M [00:01<00:01, 11.6M Bytes/s]
Get sub-02/f .. _bold.nii.gz:  62%|█▊ | 18.0M/29.2M [00:01<00:01, 10.1M Bytes/s]
Get sub-02/f .. _bold.nii.gz:  68%|██ | 19.8M/29.2M [00:01<00:00, 10.1M Bytes/s]
Get sub-02/f .. _bold.nii.gz:  74%|██▏| 21.7M/29.2M [00:02<00:00, 10.2M Bytes/s]
Get sub-02/f .. _bold.nii.gz:  80%|██▍| 23.5M/29.2M [00:02<00:00, 10.8M Bytes/s]
Get sub-02/f .. _bold.nii.gz:  85%|██▌| 24.9M/29.2M [00:02<00:00, 11.6M Bytes/s]
Get sub-02/f .. _bold.nii.gz:  93%|██▊| 27.1M/29.2M [00:02<00:00, 9.92M Bytes/s]
Get sub-02/f .. _bold.nii.gz:  99%|██▉| 29.0M/29.2M [00:02<00:00, 11.1M Bytes/s]
Total:  18%|████▉                       | 109M/618M [00:13<01:05, 7.81M Bytes/s]
Get sub-02/f .. _bold.nii.gz:   0%|            | 0.00/29.2M [00:00<?, ? Bytes/s]
Get sub-02/f .. _bold.nii.gz:   6%|▏  | 1.83M/29.2M [00:00<00:02, 10.6M Bytes/s]
Get sub-02/f .. _bold.nii.gz:  13%|▍  | 3.66M/29.2M [00:00<00:02, 10.5M Bytes/s]
Get sub-02/f .. _bold.nii.gz:  19%|▌  | 5.51M/29.2M [00:00<00:02, 10.5M Bytes/s]
Get sub-02/f .. _bold.nii.gz:  25%|▊  | 7.35M/29.2M [00:00<00:02, 10.5M Bytes/s]
Get sub-02/f .. _bold.nii.gz:  31%|▉  | 9.18M/29.2M [00:00<00:01, 10.5M Bytes/s]
Get sub-02/f .. _bold.nii.gz:  38%|█▏ | 11.0M/29.2M [00:01<00:01, 10.5M Bytes/s]
Get sub-02/f .. _bold.nii.gz:  44%|█▎ | 12.9M/29.2M [00:01<00:01, 10.5M Bytes/s]
Get sub-02/f .. _bold.nii.gz:  50%|█▌ | 14.7M/29.2M [00:01<00:01, 10.5M Bytes/s]
Get sub-02/f .. _bold.nii.gz:  57%|█▋ | 16.5M/29.2M [00:01<00:01, 10.5M Bytes/s]
Get sub-02/f .. _bold.nii.gz:  63%|█▉ | 18.4M/29.2M [00:01<00:01, 10.6M Bytes/s]
Get sub-02/f .. _bold.nii.gz:  69%|██ | 20.2M/29.2M [00:01<00:00, 10.5M Bytes/s]
Get sub-02/f .. _bold.nii.gz:  76%|██▎| 22.1M/29.2M [00:02<00:00, 10.5M Bytes/s]
Get sub-02/f .. _bold.nii.gz:  82%|██▍| 23.9M/29.2M [00:02<00:00, 11.0M Bytes/s]
Get sub-02/f .. _bold.nii.gz:  86%|██▌| 25.2M/29.2M [00:02<00:00, 11.6M Bytes/s]
Get sub-02/f .. _bold.nii.gz:  94%|██▊| 27.6M/29.2M [00:02<00:00, 10.2M Bytes/s]
                                                                                
Get sub-02/f .. _bold.nii.gz:   0%|            | 0.00/29.2M [00:00<?, ? Bytes/s]
                                                                                
Get sub-05/a .. 5_T1w.nii.gz:   0%|            | 0.00/10.8M [00:00<?, ? Bytes/s]
Get sub-05/a .. 5_T1w.nii.gz:  17%|▌  | 1.84M/10.8M [00:00<00:00, 10.6M Bytes/s]
Get sub-05/a .. 5_T1w.nii.gz:  34%|█  | 3.69M/10.8M [00:00<00:00, 10.5M Bytes/s]
Get sub-05/a .. 5_T1w.nii.gz:  51%|█▌ | 5.52M/10.8M [00:00<00:00, 10.5M Bytes/s]
Get sub-05/a .. 5_T1w.nii.gz:  68%|██ | 7.36M/10.8M [00:00<00:00, 10.6M Bytes/s]
Get sub-05/a .. 5_T1w.nii.gz:  86%|██▌| 9.21M/10.8M [00:00<00:00, 10.5M Bytes/s]
Total:  24%|██████▋                     | 148M/618M [00:18<00:57, 8.15M Bytes/s]
Get sub-05/f .. _bold.nii.gz:   0%|            | 0.00/29.7M [00:00<?, ? Bytes/s]
Get sub-05/f .. _bold.nii.gz:   6%|▏  | 1.83M/29.7M [00:00<00:02, 10.5M Bytes/s]
Get sub-05/f .. _bold.nii.gz:  12%|▎  | 3.69M/29.7M [00:00<00:02, 10.6M Bytes/s]
Get sub-05/f .. _bold.nii.gz:  19%|▌  | 5.52M/29.7M [00:00<00:02, 10.6M Bytes/s]
Get sub-05/f .. _bold.nii.gz:  25%|▋  | 7.29M/29.7M [00:00<00:01, 12.5M Bytes/s]
Get sub-05/f .. _bold.nii.gz:  31%|▉  | 9.21M/29.7M [00:00<00:02, 10.1M Bytes/s]
Get sub-05/f .. _bold.nii.gz:  37%|█  | 11.1M/29.7M [00:01<00:01, 10.3M Bytes/s]
Get sub-05/f .. _bold.nii.gz:  43%|█▎ | 12.9M/29.7M [00:01<00:01, 10.6M Bytes/s]
Get sub-05/f .. _bold.nii.gz:  50%|█▍ | 14.7M/29.7M [00:01<00:01, 12.2M Bytes/s]
Get sub-05/f .. _bold.nii.gz:  56%|█▋ | 16.6M/29.7M [00:01<00:01, 10.0M Bytes/s]
Get sub-05/f .. _bold.nii.gz:  62%|█▊ | 18.5M/29.7M [00:01<00:01, 10.8M Bytes/s]
Get sub-05/f .. _bold.nii.gz:  67%|██ | 19.8M/29.7M [00:01<00:00, 11.3M Bytes/s]
Get sub-05/f .. _bold.nii.gz:  75%|██▏| 22.2M/29.7M [00:02<00:00, 10.0M Bytes/s]
Get sub-05/f .. _bold.nii.gz:  81%|██▍| 24.0M/29.7M [00:02<00:00, 11.2M Bytes/s]
Get sub-05/f .. _bold.nii.gz:  87%|██▌| 25.9M/29.7M [00:02<00:00, 10.1M Bytes/s]
Get sub-05/f .. _bold.nii.gz:  94%|██▊| 27.8M/29.7M [00:02<00:00, 10.3M Bytes/s]
                                                                                
Get sub-05/f .. _bold.nii.gz:   0%|            | 0.00/29.7M [00:00<?, ? Bytes/s]
Get sub-05/f .. _bold.nii.gz:   6%|▏  | 1.79M/29.7M [00:00<00:02, 10.3M Bytes/s]
Get sub-05/f .. _bold.nii.gz:  12%|▎  | 3.68M/29.7M [00:00<00:02, 10.6M Bytes/s]
Get sub-05/f .. _bold.nii.gz:  19%|▌  | 5.56M/29.7M [00:00<00:02, 10.7M Bytes/s]
Get sub-05/f .. _bold.nii.gz:  25%|▊  | 7.46M/29.7M [00:00<00:02, 10.8M Bytes/s]
Get sub-05/f .. _bold.nii.gz:  32%|▉  | 9.37M/29.7M [00:00<00:01, 10.8M Bytes/s]
Get sub-05/f .. _bold.nii.gz:  38%|█▏ | 11.3M/29.7M [00:01<00:01, 10.8M Bytes/s]
Get sub-05/f .. _bold.nii.gz:  45%|█▎ | 13.2M/29.7M [00:01<00:01, 11.1M Bytes/s]
Get sub-05/f .. _bold.nii.gz:  51%|█▌ | 15.1M/29.7M [00:01<00:01, 12.7M Bytes/s]
Get sub-05/f .. _bold.nii.gz:  58%|█▋ | 17.1M/29.7M [00:01<00:01, 10.5M Bytes/s]
Get sub-05/f .. _bold.nii.gz:  64%|█▉ | 19.0M/29.7M [00:01<00:00, 11.3M Bytes/s]
Get sub-05/f .. _bold.nii.gz:  69%|██ | 20.4M/29.7M [00:01<00:00, 11.9M Bytes/s]
Get sub-05/f .. _bold.nii.gz:  77%|██▎| 23.0M/29.7M [00:02<00:00, 10.6M Bytes/s]
Get sub-05/f .. _bold.nii.gz:  84%|██▌| 24.9M/29.7M [00:02<00:00, 11.9M Bytes/s]
Get sub-05/f .. _bold.nii.gz:  91%|██▋| 26.9M/29.7M [00:02<00:00, 10.7M Bytes/s]
Get sub-05/f .. _bold.nii.gz:  98%|██▉| 29.0M/29.7M [00:02<00:00, 10.9M Bytes/s]
                                                                                
Get sub-05/f .. _bold.nii.gz:   0%|            | 0.00/29.7M [00:00<?, ? Bytes/s]
Total:  34%|█████████▍                  | 208M/618M [00:24<00:48, 8.52M Bytes/s]
Get sub-09/a .. 9_T1w.nii.gz:   0%|            | 0.00/10.8M [00:00<?, ? Bytes/s]
Get sub-09/a .. 9_T1w.nii.gz:  17%|▍  | 1.79M/10.8M [00:00<00:00, 17.8M Bytes/s]
Get sub-09/a .. 9_T1w.nii.gz:  38%|█▏ | 4.05M/10.8M [00:00<00:00, 11.0M Bytes/s]
Get sub-09/a .. 9_T1w.nii.gz:  57%|█▋ | 6.11M/10.8M [00:00<00:00, 11.4M Bytes/s]
Get sub-09/a .. 9_T1w.nii.gz:  76%|██▎| 8.14M/10.8M [00:00<00:00, 13.7M Bytes/s]
Get sub-09/a .. 9_T1w.nii.gz:  95%|██▊| 10.3M/10.8M [00:00<00:00, 11.2M Bytes/s]
                                                                                
Get sub-09/f .. _bold.nii.gz:   0%|            | 0.00/29.2M [00:00<?, ? Bytes/s]
Get sub-09/f .. _bold.nii.gz:   7%|▏  | 2.10M/29.2M [00:00<00:02, 13.3M Bytes/s]
Get sub-09/f .. _bold.nii.gz:  13%|▍  | 3.81M/29.2M [00:00<00:01, 15.1M Bytes/s]
Get sub-09/f .. _bold.nii.gz:  22%|▋  | 6.37M/29.2M [00:00<00:01, 11.6M Bytes/s]
Get sub-09/f .. _bold.nii.gz:  29%|▊  | 8.41M/29.2M [00:00<00:01, 13.9M Bytes/s]
Get sub-09/f .. _bold.nii.gz:  37%|█  | 10.7M/29.2M [00:00<00:01, 11.7M Bytes/s]
Get sub-09/f .. _bold.nii.gz:  44%|█▎ | 12.8M/29.2M [00:00<00:01, 13.7M Bytes/s]
Get sub-09/f .. _bold.nii.gz:  52%|█▌ | 15.1M/29.2M [00:01<00:01, 13.2M Bytes/s]
Get sub-09/f .. _bold.nii.gz:  59%|█▊ | 17.3M/29.2M [00:01<00:00, 13.3M Bytes/s]
Get sub-09/f .. _bold.nii.gz:  67%|██ | 19.7M/29.2M [00:01<00:00, 11.7M Bytes/s]
Get sub-09/f .. _bold.nii.gz:  75%|██▏| 21.8M/29.2M [00:01<00:00, 13.5M Bytes/s]
Get sub-09/f .. _bold.nii.gz:  83%|██▍| 24.3M/29.2M [00:01<00:00, 12.1M Bytes/s]
Get sub-09/f .. _bold.nii.gz:  90%|██▋| 26.4M/29.2M [00:02<00:00, 13.8M Bytes/s]
Get sub-09/f .. _bold.nii.gz:  99%|██▉| 29.0M/29.2M [00:02<00:00, 12.5M Bytes/s]
                                                                                
Get sub-09/f .. _bold.nii.gz:   0%|            | 0.00/29.2M [00:00<?, ? Bytes/s]
Get sub-09/f .. _bold.nii.gz:   6%|▏  | 1.70M/29.2M [00:00<00:01, 16.9M Bytes/s]
Get sub-09/f .. _bold.nii.gz:  13%|▍  | 3.85M/29.2M [00:00<00:02, 10.5M Bytes/s]
Get sub-09/f .. _bold.nii.gz:  19%|▌  | 5.58M/29.2M [00:00<00:02, 10.2M Bytes/s]
Get sub-09/f .. _bold.nii.gz:  25%|▊  | 7.31M/29.2M [00:00<00:02, 10.1M Bytes/s]
Get sub-09/f .. _bold.nii.gz:  31%|▉  | 9.13M/29.2M [00:00<00:01, 10.2M Bytes/s]
Get sub-09/f .. _bold.nii.gz:  38%|█▏ | 11.0M/29.2M [00:01<00:01, 10.3M Bytes/s]
Get sub-09/f .. _bold.nii.gz:  44%|█▎ | 12.8M/29.2M [00:01<00:01, 10.9M Bytes/s]
Get sub-09/f .. _bold.nii.gz:  49%|█▍ | 14.4M/29.2M [00:01<00:01, 12.0M Bytes/s]
Get sub-09/f .. _bold.nii.gz:  57%|█▋ | 16.7M/29.2M [00:01<00:01, 10.3M Bytes/s]
Get sub-09/f .. _bold.nii.gz:  64%|█▉ | 18.6M/29.2M [00:01<00:00, 11.4M Bytes/s]
Get sub-09/f .. _bold.nii.gz:  70%|██ | 20.6M/29.2M [00:01<00:00, 10.5M Bytes/s]
Get sub-09/f .. _bold.nii.gz:  77%|██▎| 22.4M/29.2M [00:02<00:00, 12.0M Bytes/s]
Get sub-09/f .. _bold.nii.gz:  84%|██▌| 24.6M/29.2M [00:02<00:00, 12.0M Bytes/s]
Get sub-09/f .. _bold.nii.gz:  91%|██▋| 26.5M/29.2M [00:02<00:00, 11.9M Bytes/s]
Get sub-09/f .. _bold.nii.gz:  98%|██▉| 28.7M/29.2M [00:02<00:00, 10.6M Bytes/s]
Total:  45%|████████████▌               | 277M/618M [00:31<00:38, 8.83M Bytes/s]
Get sub-04/a .. 4_T1w.nii.gz:   0%|            | 0.00/10.7M [00:00<?, ? Bytes/s]
Get sub-04/a .. 4_T1w.nii.gz:  19%|▌  | 2.06M/10.7M [00:00<00:00, 11.8M Bytes/s]
Get sub-04/a .. 4_T1w.nii.gz:  39%|█▏ | 4.19M/10.7M [00:00<00:00, 12.0M Bytes/s]
Get sub-04/a .. 4_T1w.nii.gz:  59%|█▊ | 6.30M/10.7M [00:00<00:00, 15.0M Bytes/s]
Get sub-04/a .. 4_T1w.nii.gz:  79%|██▎| 8.49M/10.7M [00:00<00:00, 11.7M Bytes/s]
Get sub-04/a .. 4_T1w.nii.gz:  98%|██▉| 10.6M/10.7M [00:00<00:00, 13.8M Bytes/s]
                                                                                
Get sub-04/f .. _bold.nii.gz:   0%|            | 0.00/29.1M [00:00<?, ? Bytes/s]
Get sub-04/f .. _bold.nii.gz:   7%|▏  | 2.16M/29.1M [00:00<00:01, 21.6M Bytes/s]
Get sub-04/f .. _bold.nii.gz:  15%|▍  | 4.40M/29.1M [00:00<00:02, 11.6M Bytes/s]
Get sub-04/f .. _bold.nii.gz:  23%|▋  | 6.56M/29.1M [00:00<00:01, 14.7M Bytes/s]
Get sub-04/f .. _bold.nii.gz:  33%|▉  | 9.54M/29.1M [00:00<00:01, 11.0M Bytes/s]
Get sub-04/f .. _bold.nii.gz:  39%|█▏ | 11.4M/29.1M [00:01<00:01, 9.63M Bytes/s]
Get sub-04/f .. _bold.nii.gz:  45%|█▎ | 13.0M/29.1M [00:01<00:01, 9.86M Bytes/s]
Get sub-04/f .. _bold.nii.gz:  50%|█▌ | 14.6M/29.1M [00:01<00:01, 11.0M Bytes/s]
Get sub-04/f .. _bold.nii.gz:  56%|█▋ | 16.3M/29.1M [00:01<00:01, 8.96M Bytes/s]
Get sub-04/f .. _bold.nii.gz:  62%|█▊ | 17.9M/29.1M [00:01<00:01, 9.84M Bytes/s]
Get sub-04/f .. _bold.nii.gz:  68%|██ | 19.6M/29.1M [00:01<00:01, 9.09M Bytes/s]
Get sub-04/f .. _bold.nii.gz:  73%|██▏| 21.3M/29.1M [00:02<00:00, 9.32M Bytes/s]
Get sub-04/f .. _bold.nii.gz:  79%|██▍| 23.0M/29.1M [00:02<00:00, 10.6M Bytes/s]
Get sub-04/f .. _bold.nii.gz:  85%|██▌| 24.8M/29.1M [00:02<00:00, 9.23M Bytes/s]
Get sub-04/f .. _bold.nii.gz:  91%|██▋| 26.5M/29.1M [00:02<00:00, 9.50M Bytes/s]
Get sub-04/f .. _bold.nii.gz:  97%|██▉| 28.3M/29.1M [00:02<00:00, 9.66M Bytes/s]
                                                                                
Get sub-04/f .. _bold.nii.gz:   0%|            | 0.00/29.1M [00:00<?, ? Bytes/s]
Get sub-04/f .. _bold.nii.gz:   6%|▏  | 1.77M/29.1M [00:00<00:02, 10.2M Bytes/s]
Get sub-04/f .. _bold.nii.gz:  12%|▎  | 3.54M/29.1M [00:00<00:02, 10.2M Bytes/s]
Get sub-04/f .. _bold.nii.gz:  18%|▌  | 5.34M/29.1M [00:00<00:02, 10.2M Bytes/s]
Get sub-04/f .. _bold.nii.gz:  25%|▋  | 7.14M/29.1M [00:00<00:02, 10.3M Bytes/s]
Get sub-04/f .. _bold.nii.gz:  31%|▉  | 8.95M/29.1M [00:00<00:01, 10.3M Bytes/s]
Get sub-04/f .. _bold.nii.gz:  37%|█  | 10.8M/29.1M [00:01<00:01, 10.3M Bytes/s]
Get sub-04/f .. _bold.nii.gz:  43%|█▎ | 12.6M/29.1M [00:01<00:01, 10.4M Bytes/s]
Get sub-04/f .. _bold.nii.gz:  50%|█▍ | 14.4M/29.1M [00:01<00:01, 10.4M Bytes/s]
Get sub-04/f .. _bold.nii.gz:  56%|█▋ | 16.3M/29.1M [00:01<00:01, 10.4M Bytes/s]
Get sub-04/f .. _bold.nii.gz:  62%|█▊ | 18.1M/29.1M [00:01<00:01, 10.8M Bytes/s]
Get sub-04/f .. _bold.nii.gz:  68%|██ | 19.9M/29.1M [00:01<00:00, 12.2M Bytes/s]
Get sub-04/f .. _bold.nii.gz:  75%|██▎| 21.9M/29.1M [00:02<00:00, 10.1M Bytes/s]
Get sub-04/f .. _bold.nii.gz:  82%|██▍| 23.7M/29.1M [00:02<00:00, 11.0M Bytes/s]
Get sub-04/f .. _bold.nii.gz:  86%|██▌| 25.0M/29.1M [00:02<00:00, 11.3M Bytes/s]
Get sub-04/f .. _bold.nii.gz:  95%|██▊| 27.5M/29.1M [00:02<00:00, 10.1M Bytes/s]
                                                                                
Get sub-04/f .. _bold.nii.gz:   0%|            | 0.00/29.1M [00:00<?, ? Bytes/s]
                                                                                
Get sub-01/a .. 1_T1w.nii.gz:   0%|            | 0.00/10.6M [00:00<?, ? Bytes/s]
Get sub-01/a .. 1_T1w.nii.gz:  17%|▌  | 1.85M/10.6M [00:00<00:00, 10.7M Bytes/s]
Get sub-01/a .. 1_T1w.nii.gz:  35%|█  | 3.74M/10.6M [00:00<00:00, 10.7M Bytes/s]
Get sub-01/a .. 1_T1w.nii.gz:  53%|█▌ | 5.64M/10.6M [00:00<00:00, 10.8M Bytes/s]
Get sub-01/a .. 1_T1w.nii.gz:  63%|█▉ | 6.72M/10.6M [00:00<00:00, 8.83M Bytes/s]
Get sub-01/a .. 1_T1w.nii.gz:  83%|██▍| 8.81M/10.6M [00:00<00:00, 11.7M Bytes/s]
Get sub-01/a .. 1_T1w.nii.gz:  97%|██▉| 10.3M/10.6M [00:01<00:00, 9.21M Bytes/s]
Total:  58%|████████████████▏           | 356M/618M [00:40<00:29, 8.85M Bytes/s]
Get sub-01/f .. _bold.nii.gz:   0%|            | 0.00/28.1M [00:00<?, ? Bytes/s]
Get sub-01/f .. _bold.nii.gz:   5%|▏  | 1.36M/28.1M [00:00<00:03, 7.99M Bytes/s]
Get sub-01/f .. _bold.nii.gz:  10%|▎  | 2.74M/28.1M [00:00<00:03, 7.94M Bytes/s]
Get sub-01/f .. _bold.nii.gz:  15%|▍  | 4.12M/28.1M [00:00<00:03, 7.93M Bytes/s]
Get sub-01/f .. _bold.nii.gz:  20%|▌  | 5.52M/28.1M [00:00<00:02, 8.31M Bytes/s]
Get sub-01/f .. _bold.nii.gz:  25%|▋  | 6.96M/28.1M [00:00<00:02, 7.95M Bytes/s]
Get sub-01/f .. _bold.nii.gz:  30%|▉  | 8.40M/28.1M [00:01<00:02, 8.04M Bytes/s]
Get sub-01/f .. _bold.nii.gz:  35%|█  | 9.83M/28.1M [00:01<00:02, 8.81M Bytes/s]
Get sub-01/f .. _bold.nii.gz:  39%|█▏ | 11.0M/28.1M [00:01<00:01, 9.40M Bytes/s]
Get sub-01/f .. _bold.nii.gz:  46%|█▎ | 12.8M/28.1M [00:01<00:01, 7.82M Bytes/s]
Get sub-01/f .. _bold.nii.gz:  51%|█▌ | 14.3M/28.1M [00:01<00:01, 9.02M Bytes/s]
Get sub-01/f .. _bold.nii.gz:  56%|█▋ | 15.8M/28.1M [00:01<00:01, 7.92M Bytes/s]
Get sub-01/f .. _bold.nii.gz:  62%|█▊ | 17.3M/28.1M [00:02<00:01, 8.11M Bytes/s]
Get sub-01/f .. _bold.nii.gz:  67%|██ | 18.8M/28.1M [00:02<00:01, 8.28M Bytes/s]
Get sub-01/f .. _bold.nii.gz:  72%|██▏| 20.3M/28.1M [00:02<00:00, 8.41M Bytes/s]
Get sub-01/f .. _bold.nii.gz:  78%|██▎| 21.8M/28.1M [00:02<00:00, 8.54M Bytes/s]
Get sub-01/f .. _bold.nii.gz:  83%|██▍| 23.4M/28.1M [00:02<00:00, 8.62M Bytes/s]
Get sub-01/f .. _bold.nii.gz:  89%|██▋| 24.9M/28.1M [00:02<00:00, 8.68M Bytes/s]
Get sub-01/f .. _bold.nii.gz:  94%|██▊| 26.5M/28.1M [00:03<00:00, 8.72M Bytes/s]
Get sub-01/f .. _bold.nii.gz: 100%|██▉| 28.0M/28.1M [00:03<00:00, 8.81M Bytes/s]
                                                                                
Get sub-01/f .. _bold.nii.gz:   0%|            | 0.00/28.1M [00:00<?, ? Bytes/s]
Get sub-01/f .. _bold.nii.gz:   6%|▏  | 1.55M/28.1M [00:00<00:02, 9.47M Bytes/s]
Get sub-01/f .. _bold.nii.gz:  11%|▎  | 3.14M/28.1M [00:00<00:02, 8.91M Bytes/s]
Get sub-01/f .. _bold.nii.gz:  17%|▌  | 4.71M/28.1M [00:00<00:02, 8.95M Bytes/s]
Get sub-01/f .. _bold.nii.gz:  22%|▋  | 6.28M/28.1M [00:00<00:02, 9.70M Bytes/s]
Get sub-01/f .. _bold.nii.gz:  27%|▊  | 7.57M/28.1M [00:00<00:01, 10.5M Bytes/s]
Get sub-01/f .. _bold.nii.gz:  34%|█  | 9.46M/28.1M [00:01<00:02, 8.56M Bytes/s]
Get sub-01/f .. _bold.nii.gz:  38%|█▏ | 10.6M/28.1M [00:01<00:01, 9.09M Bytes/s]
Get sub-01/f .. _bold.nii.gz:  45%|█▎ | 12.7M/28.1M [00:01<00:01, 8.77M Bytes/s]
Get sub-01/f .. _bold.nii.gz:  51%|█▌ | 14.2M/28.1M [00:01<00:01, 8.89M Bytes/s]
Get sub-01/f .. _bold.nii.gz:  56%|█▋ | 15.8M/28.1M [00:01<00:01, 8.96M Bytes/s]
Get sub-01/f .. _bold.nii.gz:  62%|█▊ | 17.5M/28.1M [00:01<00:01, 9.00M Bytes/s]
Get sub-01/f .. _bold.nii.gz:  68%|██ | 19.1M/28.1M [00:02<00:01, 9.06M Bytes/s]
Get sub-01/f .. _bold.nii.gz:  73%|██▏| 20.7M/28.1M [00:02<00:00, 9.09M Bytes/s]
Get sub-01/f .. _bold.nii.gz:  79%|██▎| 22.3M/28.1M [00:02<00:00, 9.08M Bytes/s]
Get sub-01/f .. _bold.nii.gz:  86%|██▌| 24.2M/28.1M [00:02<00:00, 9.28M Bytes/s]
Get sub-01/f .. _bold.nii.gz:  91%|██▋| 25.5M/28.1M [00:02<00:00, 9.13M Bytes/s]
Get sub-01/f .. _bold.nii.gz:  96%|██▉| 27.1M/28.1M [00:02<00:00, 9.21M Bytes/s]
                                                                                
Get sub-08/a .. 8_T1w.nii.gz:   0%|            | 0.00/10.6M [00:00<?, ? Bytes/s]
Get sub-08/a .. 8_T1w.nii.gz:  15%|▍  | 1.62M/10.6M [00:00<00:00, 9.23M Bytes/s]
Get sub-08/a .. 8_T1w.nii.gz:  31%|▉  | 3.23M/10.6M [00:00<00:00, 9.24M Bytes/s]
Get sub-08/a .. 8_T1w.nii.gz:  46%|█▎ | 4.83M/10.6M [00:00<00:00, 9.23M Bytes/s]
Get sub-08/a .. 8_T1w.nii.gz:  61%|█▊ | 6.44M/10.6M [00:00<00:00, 9.20M Bytes/s]
Get sub-08/a .. 8_T1w.nii.gz:  76%|██▎| 8.06M/10.6M [00:00<00:00, 9.21M Bytes/s]
Get sub-08/a .. 8_T1w.nii.gz:  92%|██▋| 9.67M/10.6M [00:01<00:00, 9.20M Bytes/s]
                                                                                
Get sub-08/f .. _bold.nii.gz:   0%|            | 0.00/28.6M [00:00<?, ? Bytes/s]
Get sub-08/f .. _bold.nii.gz:   5%|▏  | 1.57M/28.6M [00:00<00:03, 9.01M Bytes/s]
Get sub-08/f .. _bold.nii.gz:  11%|▎  | 3.19M/28.6M [00:00<00:02, 9.13M Bytes/s]
Get sub-08/f .. _bold.nii.gz:  17%|▌  | 4.80M/28.6M [00:00<00:02, 9.17M Bytes/s]
Get sub-08/f .. _bold.nii.gz:  22%|▋  | 6.41M/28.6M [00:00<00:02, 9.18M Bytes/s]
Get sub-08/f .. _bold.nii.gz:  28%|▊  | 8.02M/28.6M [00:00<00:02, 9.23M Bytes/s]
Get sub-08/f .. _bold.nii.gz:  34%|█  | 9.64M/28.6M [00:01<00:02, 9.20M Bytes/s]
Get sub-08/f .. _bold.nii.gz:  39%|█▏ | 11.2M/28.6M [00:01<00:01, 9.20M Bytes/s]
Get sub-08/f .. _bold.nii.gz:  45%|█▎ | 12.8M/28.6M [00:01<00:01, 9.24M Bytes/s]
Get sub-08/f .. _bold.nii.gz:  51%|█▌ | 14.5M/28.6M [00:01<00:01, 9.18M Bytes/s]
Get sub-08/f .. _bold.nii.gz:  56%|█▋ | 16.1M/28.6M [00:01<00:01, 9.47M Bytes/s]
Get sub-08/f .. _bold.nii.gz:  61%|█▊ | 17.4M/28.6M [00:01<00:01, 10.2M Bytes/s]
Get sub-08/f .. _bold.nii.gz:  67%|██ | 19.3M/28.6M [00:02<00:01, 8.89M Bytes/s]
Get sub-08/f .. _bold.nii.gz:  73%|██▏| 20.9M/28.6M [00:02<00:00, 9.63M Bytes/s]
Get sub-08/f .. _bold.nii.gz:  77%|██▎| 22.0M/28.6M [00:02<00:00, 9.89M Bytes/s]
Get sub-08/f .. _bold.nii.gz:  84%|██▌| 24.2M/28.6M [00:02<00:00, 8.76M Bytes/s]
Get sub-08/f .. _bold.nii.gz:  90%|██▋| 25.8M/28.6M [00:02<00:00, 9.92M Bytes/s]
Get sub-08/f .. _bold.nii.gz:  96%|██▊| 27.4M/28.6M [00:02<00:00, 8.72M Bytes/s]
Total:  73%|████████████████████▍       | 452M/618M [00:52<00:19, 8.62M Bytes/s]
Get sub-08/f .. _bold.nii.gz:   0%|            | 0.00/28.6M [00:00<?, ? Bytes/s]
Get sub-08/f .. _bold.nii.gz:   6%|▏  | 1.62M/28.6M [00:00<00:02, 9.31M Bytes/s]
Get sub-08/f .. _bold.nii.gz:  11%|▎  | 3.24M/28.6M [00:00<00:02, 9.30M Bytes/s]
Get sub-08/f .. _bold.nii.gz:  17%|▌  | 4.86M/28.6M [00:00<00:02, 9.29M Bytes/s]
Get sub-08/f .. _bold.nii.gz:  23%|▋  | 6.50M/28.6M [00:00<00:02, 9.28M Bytes/s]
Get sub-08/f .. _bold.nii.gz:  28%|▊  | 7.99M/28.6M [00:00<00:01, 10.6M Bytes/s]
Get sub-08/f .. _bold.nii.gz:  34%|█  | 9.76M/28.6M [00:01<00:02, 9.03M Bytes/s]
Get sub-08/f .. _bold.nii.gz:  40%|█▏ | 11.4M/28.6M [00:01<00:01, 9.21M Bytes/s]
Get sub-08/f .. _bold.nii.gz:  46%|█▎ | 13.1M/28.6M [00:01<00:01, 9.21M Bytes/s]
Get sub-08/f .. _bold.nii.gz:  51%|█▌ | 14.7M/28.6M [00:01<00:01, 9.25M Bytes/s]
Get sub-08/f .. _bold.nii.gz:  57%|█▋ | 16.3M/28.6M [00:01<00:01, 9.81M Bytes/s]
Get sub-08/f .. _bold.nii.gz:  61%|█▊ | 17.5M/28.6M [00:01<00:01, 10.2M Bytes/s]
Get sub-08/f .. _bold.nii.gz:  69%|██ | 19.7M/28.6M [00:02<00:00, 9.08M Bytes/s]
Get sub-08/f .. _bold.nii.gz:  75%|██▏| 21.3M/28.6M [00:02<00:00, 10.1M Bytes/s]
Get sub-08/f .. _bold.nii.gz:  80%|██▍| 23.0M/28.6M [00:02<00:00, 9.12M Bytes/s]
Get sub-08/f .. _bold.nii.gz:  86%|██▌| 24.7M/28.6M [00:02<00:00, 9.28M Bytes/s]
Get sub-08/f .. _bold.nii.gz:  92%|██▊| 26.4M/28.6M [00:02<00:00, 9.39M Bytes/s]
Get sub-08/f .. _bold.nii.gz:  98%|██▉| 28.1M/28.6M [00:02<00:00, 9.53M Bytes/s]
Total:  78%|█████████████████████▊      | 481M/618M [00:55<00:15, 8.61M Bytes/s]
Get sub-06/a .. 6_T1w.nii.gz:   0%|            | 0.00/10.6M [00:00<?, ? Bytes/s]
Get sub-06/a .. 6_T1w.nii.gz:  16%|▍  | 1.67M/10.6M [00:00<00:00, 9.64M Bytes/s]
Get sub-06/a .. 6_T1w.nii.gz:  32%|▉  | 3.39M/10.6M [00:00<00:00, 9.76M Bytes/s]
Get sub-06/a .. 6_T1w.nii.gz:  48%|█▍ | 5.13M/10.6M [00:00<00:00, 9.84M Bytes/s]
Get sub-06/a .. 6_T1w.nii.gz:  65%|█▉ | 6.87M/10.6M [00:00<00:00, 10.1M Bytes/s]
Get sub-06/a .. 6_T1w.nii.gz:  81%|██▍| 8.64M/10.6M [00:00<00:00, 9.93M Bytes/s]
Get sub-06/a .. 6_T1w.nii.gz:  98%|██▉| 10.4M/10.6M [00:01<00:00, 9.99M Bytes/s]
                                                                                
Get sub-06/f .. _bold.nii.gz:   0%|            | 0.00/29.4M [00:00<?, ? Bytes/s]
Get sub-06/f .. _bold.nii.gz:   6%|▏  | 1.76M/29.4M [00:00<00:02, 10.1M Bytes/s]
Get sub-06/f .. _bold.nii.gz:  12%|▎  | 3.57M/29.4M [00:00<00:02, 10.3M Bytes/s]
Get sub-06/f .. _bold.nii.gz:  18%|▌  | 5.39M/29.4M [00:00<00:02, 10.3M Bytes/s]
Get sub-06/f .. _bold.nii.gz:  26%|▊  | 7.76M/29.4M [00:00<00:01, 11.0M Bytes/s]
Get sub-06/f .. _bold.nii.gz:  31%|▉  | 8.98M/29.4M [00:01<00:02, 7.39M Bytes/s]
Get sub-06/f .. _bold.nii.gz:  38%|█▏ | 11.1M/29.4M [00:01<00:01, 9.85M Bytes/s]
Get sub-06/f .. _bold.nii.gz:  44%|█▎ | 12.8M/29.4M [00:01<00:01, 8.71M Bytes/s]
Get sub-06/f .. _bold.nii.gz:  48%|█▍ | 14.2M/29.4M [00:01<00:01, 8.39M Bytes/s]
Get sub-06/f .. _bold.nii.gz:  53%|█▌ | 15.5M/29.4M [00:01<00:01, 8.28M Bytes/s]
Get sub-06/f .. _bold.nii.gz:  58%|█▋ | 16.9M/29.4M [00:01<00:01, 8.22M Bytes/s]
Get sub-06/f .. _bold.nii.gz:  63%|█▉ | 18.4M/29.4M [00:02<00:01, 8.20M Bytes/s]
Get sub-06/f .. _bold.nii.gz:  67%|██ | 19.8M/29.4M [00:02<00:01, 8.36M Bytes/s]
Get sub-06/f .. _bold.nii.gz:  73%|██▏| 21.3M/29.4M [00:02<00:00, 8.27M Bytes/s]
Get sub-06/f .. _bold.nii.gz:  78%|██▎| 22.8M/29.4M [00:02<00:00, 8.36M Bytes/s]
Get sub-06/f .. _bold.nii.gz:  83%|██▍| 24.3M/29.4M [00:02<00:00, 8.97M Bytes/s]
Get sub-06/f .. _bold.nii.gz:  87%|██▌| 25.7M/29.4M [00:02<00:00, 9.91M Bytes/s]
Get sub-06/f .. _bold.nii.gz:  93%|██▊| 27.5M/29.4M [00:03<00:00, 8.25M Bytes/s]
Get sub-06/f .. _bold.nii.gz:  99%|██▉| 29.0M/29.4M [00:03<00:00, 9.32M Bytes/s]
Total:  84%|███████████████████████▌    | 521M/618M [01:00<00:11, 8.55M Bytes/s]
Get sub-06/f .. _bold.nii.gz:   0%|            | 0.00/29.4M [00:00<?, ? Bytes/s]
Get sub-06/f .. _bold.nii.gz:   5%|▏  | 1.60M/29.4M [00:00<00:03, 9.19M Bytes/s]
Get sub-06/f .. _bold.nii.gz:  11%|▎  | 3.23M/29.4M [00:00<00:02, 9.25M Bytes/s]
Get sub-06/f .. _bold.nii.gz:  17%|▍  | 4.86M/29.4M [00:00<00:02, 9.31M Bytes/s]
Get sub-06/f .. _bold.nii.gz:  22%|▋  | 6.52M/29.4M [00:00<00:02, 9.37M Bytes/s]
Get sub-06/f .. _bold.nii.gz:  28%|▊  | 8.19M/29.4M [00:00<00:02, 9.44M Bytes/s]
Get sub-06/f .. _bold.nii.gz:  34%|█  | 9.87M/29.4M [00:01<00:02, 9.50M Bytes/s]
Get sub-06/f .. _bold.nii.gz:  39%|█▏ | 11.6M/29.4M [00:01<00:01, 9.65M Bytes/s]
Get sub-06/f .. _bold.nii.gz:  45%|█▎ | 13.3M/29.4M [00:01<00:01, 9.64M Bytes/s]
Get sub-06/f .. _bold.nii.gz:  51%|█▌ | 15.0M/29.4M [00:01<00:01, 9.70M Bytes/s]
Get sub-06/f .. _bold.nii.gz:  57%|█▋ | 16.7M/29.4M [00:01<00:01, 10.2M Bytes/s]
Get sub-06/f .. _bold.nii.gz:  62%|█▊ | 18.1M/29.4M [00:01<00:01, 11.0M Bytes/s]
Get sub-06/f .. _bold.nii.gz:  69%|██ | 20.2M/29.4M [00:02<00:00, 9.49M Bytes/s]
Get sub-06/f .. _bold.nii.gz:  75%|██▏| 22.0M/29.4M [00:02<00:00, 10.6M Bytes/s]
Get sub-06/f .. _bold.nii.gz:  81%|██▍| 23.8M/29.4M [00:02<00:00, 9.58M Bytes/s]
Get sub-06/f .. _bold.nii.gz:  87%|██▌| 25.6M/29.4M [00:02<00:00, 9.74M Bytes/s]
Get sub-06/f .. _bold.nii.gz:  93%|██▊| 27.4M/29.4M [00:02<00:00, 9.93M Bytes/s]
Get sub-06/f .. _bold.nii.gz:  99%|██▉| 29.2M/29.4M [00:02<00:00, 10.0M Bytes/s]
                                                                                
Get sub-03/a .. 3_T1w.nii.gz:   0%|            | 0.00/10.7M [00:00<?, ? Bytes/s]
Get sub-03/a .. 3_T1w.nii.gz:  17%|▍  | 1.78M/10.7M [00:00<00:00, 10.6M Bytes/s]
Get sub-03/a .. 3_T1w.nii.gz:  34%|█  | 3.61M/10.7M [00:00<00:00, 10.3M Bytes/s]
Get sub-03/a .. 3_T1w.nii.gz:  51%|█▌ | 5.42M/10.7M [00:00<00:00, 10.3M Bytes/s]
Get sub-03/a .. 3_T1w.nii.gz:  68%|██ | 7.24M/10.7M [00:00<00:00, 11.1M Bytes/s]
Get sub-03/a .. 3_T1w.nii.gz:  82%|██▍| 8.73M/10.7M [00:00<00:00, 12.0M Bytes/s]
Total:  91%|█████████████████████████▍  | 561M/618M [01:05<00:06, 8.55M Bytes/s]
Get sub-03/f .. _bold.nii.gz:   0%|            | 0.00/28.8M [00:00<?, ? Bytes/s]
Get sub-03/f .. _bold.nii.gz:   6%|▏  | 1.81M/28.8M [00:00<00:02, 10.4M Bytes/s]
Get sub-03/f .. _bold.nii.gz:  13%|▍  | 3.66M/28.8M [00:00<00:02, 10.5M Bytes/s]
Get sub-03/f .. _bold.nii.gz:  19%|▌  | 5.49M/28.8M [00:00<00:02, 10.5M Bytes/s]
Get sub-03/f .. _bold.nii.gz:  26%|▊  | 7.35M/28.8M [00:00<00:02, 10.3M Bytes/s]
Get sub-03/f .. _bold.nii.gz:  32%|▉  | 9.20M/28.8M [00:00<00:01, 10.4M Bytes/s]
Get sub-03/f .. _bold.nii.gz:  38%|█▏ | 10.9M/28.8M [00:01<00:01, 10.1M Bytes/s]
Get sub-03/f .. _bold.nii.gz:  45%|█▎ | 13.0M/28.8M [00:01<00:01, 8.31M Bytes/s]
Get sub-03/f .. _bold.nii.gz:  50%|█▍ | 14.3M/28.8M [00:01<00:01, 8.07M Bytes/s]
Get sub-03/f .. _bold.nii.gz:  54%|█▋ | 15.6M/28.8M [00:01<00:01, 8.09M Bytes/s]
Get sub-03/f .. _bold.nii.gz:  59%|█▊ | 17.0M/28.8M [00:01<00:01, 7.85M Bytes/s]
Get sub-03/f .. _bold.nii.gz:  64%|█▉ | 18.4M/28.8M [00:02<00:01, 7.84M Bytes/s]
Get sub-03/f .. _bold.nii.gz:  69%|██ | 19.7M/28.8M [00:02<00:01, 8.30M Bytes/s]
Get sub-03/f .. _bold.nii.gz:  73%|██▏| 21.1M/28.8M [00:02<00:00, 9.37M Bytes/s]
Get sub-03/f .. _bold.nii.gz:  78%|██▎| 22.5M/28.8M [00:02<00:00, 7.43M Bytes/s]
Get sub-03/f .. _bold.nii.gz:  83%|██▍| 23.9M/28.8M [00:02<00:00, 8.40M Bytes/s]
Get sub-03/f .. _bold.nii.gz:  88%|██▋| 25.4M/28.8M [00:02<00:00, 7.57M Bytes/s]
Get sub-03/f .. _bold.nii.gz:  93%|██▊| 26.8M/28.8M [00:03<00:00, 7.77M Bytes/s]
Get sub-03/f .. _bold.nii.gz:  98%|██▉| 28.3M/28.8M [00:03<00:00, 7.92M Bytes/s]
                                                                                
Get sub-03/f .. _bold.nii.gz:   0%|            | 0.00/28.8M [00:00<?, ? Bytes/s]
Get sub-03/f .. _bold.nii.gz:   5%|▏  | 1.45M/28.8M [00:00<00:03, 8.35M Bytes/s]
Get sub-03/f .. _bold.nii.gz:  10%|▎  | 2.93M/28.8M [00:00<00:03, 8.40M Bytes/s]
Get sub-03/f .. _bold.nii.gz:  15%|▍  | 4.40M/28.8M [00:00<00:02, 8.43M Bytes/s]
Get sub-03/f .. _bold.nii.gz:  20%|▌  | 5.89M/28.8M [00:00<00:02, 8.47M Bytes/s]
Get sub-03/f .. _bold.nii.gz:  26%|▊  | 7.40M/28.8M [00:00<00:02, 8.51M Bytes/s]
Get sub-03/f .. _bold.nii.gz:  31%|▉  | 8.91M/28.8M [00:01<00:02, 8.53M Bytes/s]
Get sub-03/f .. _bold.nii.gz:  36%|█  | 10.4M/28.8M [00:01<00:02, 8.61M Bytes/s]
Get sub-03/f .. _bold.nii.gz:  41%|█▏ | 11.9M/28.8M [00:01<00:01, 8.63M Bytes/s]
Get sub-03/f .. _bold.nii.gz:  47%|█▍ | 13.5M/28.8M [00:01<00:01, 8.66M Bytes/s]
Get sub-03/f .. _bold.nii.gz:  52%|█▌ | 15.0M/28.8M [00:01<00:01, 8.69M Bytes/s]
Get sub-03/f .. _bold.nii.gz:  57%|█▋ | 16.5M/28.8M [00:01<00:01, 8.73M Bytes/s]
Get sub-03/f .. _bold.nii.gz:  63%|█▉ | 18.1M/28.8M [00:02<00:01, 8.76M Bytes/s]
Get sub-03/f .. _bold.nii.gz:  68%|██ | 19.6M/28.8M [00:02<00:01, 8.86M Bytes/s]
Get sub-03/f .. _bold.nii.gz:  74%|██▏| 21.2M/28.8M [00:02<00:00, 8.80M Bytes/s]
Get sub-03/f .. _bold.nii.gz:  79%|██▎| 22.7M/28.8M [00:02<00:00, 8.82M Bytes/s]
Get sub-03/f .. _bold.nii.gz:  84%|██▌| 24.3M/28.8M [00:02<00:00, 9.32M Bytes/s]
Get sub-03/f .. _bold.nii.gz:  88%|██▋| 25.3M/28.8M [00:02<00:00, 9.45M Bytes/s]
Get sub-03/f .. _bold.nii.gz:  95%|██▊| 27.4M/28.8M [00:03<00:00, 8.65M Bytes/s]
                                                                                
Get sub-03/f .. _bold.nii.gz:   0%|            | 0.00/28.8M [00:00<?, ? Bytes/s]
get(ok): sub-07/anat/sub-07_T1w.nii.gz (file) [from s3-PUBLIC...]
get(ok): sub-07/func/sub-07_task-flanker_run-1_bold.nii.gz (file) [from s3-PUBLIC...]
get(ok): sub-07/func/sub-07_task-flanker_run-2_bold.nii.gz (file) [from s3-PUBLIC...]
get(ok): sub-02/anat/sub-02_T1w.nii.gz (file) [from s3-PUBLIC...]
get(ok): sub-02/func/sub-02_task-flanker_run-1_bold.nii.gz (file) [from s3-PUBLIC...]
get(ok): sub-02/func/sub-02_task-flanker_run-2_bold.nii.gz (file) [from s3-PUBLIC...]
get(ok): sub-05/anat/sub-05_T1w.nii.gz (file) [from s3-PUBLIC...]
get(ok): sub-05/func/sub-05_task-flanker_run-1_bold.nii.gz (file) [from s3-PUBLIC...]
get(ok): sub-05/func/sub-05_task-flanker_run-2_bold.nii.gz (file) [from s3-PUBLIC...]
get(ok): sub-09/anat/sub-09_T1w.nii.gz (file) [from s3-PUBLIC...]
  [17 similar messages have been suppressed; disable with datalad.ui.suppress-similar-results=off]
get(ok): sub-07 (directory)
get(ok): sub-02 (directory)
get(ok): sub-05 (directory)
get(ok): sub-09 (directory)
get(ok): sub-04 (directory)
get(ok): sub-01 (directory)
get(ok): sub-08 (directory)
get(ok): sub-06 (directory)
get(ok): sub-03 (directory)
action summary:
  get (ok: 36)

Start the workflow#

wf_first = pe.Workflow(name='level1', base_dir=exp_dir_first)
wf_first.config["execution"]["crashfile_format"] = "txt"

The following two nodes (infosource & dg) together define all inputs required for the preprocessing workflow

# get subject_id list - we only do 9 subjects here so the notebook runs in a reasonable time
subj_list = [x.split('-')[-1] for x in glob(data_dir_first+"/"+PATTERN)]
subj_list.sort()
print(subj_list)
['01', '02', '03', '04', '05', '06', '07', '08', '09']
infosource_first = pe.Node(util.IdentityInterface(fields=["subject_id"]),
                  name="infosource")
infosource_first.iterables = [("subject_id", subj_list)]
dg_first = pe.Node(
    interface=nio.DataGrabber(
        infields=["subject_id","run_id"], outfields=["struct", "func", "events"]
    ),
    name="dg_first"
)

# Specify run_ids and return a sorted filelist to ensure we match files to correct runs/tasks
dg_first.inputs.run_id = [1,2]
dg_first.inputs.sort_filelist = True
dg_first.inputs.template = "*"
dg_first.inputs.base_directory = data_dir_first

# Define arguments fill the wildcards in the below paths 
dg_first.inputs.template_args = dict(
    struct=[["subject_id","subject_id"]],
    func=[["subject_id","subject_id","run_id"]],
    events=[["subject_id","subject_id","run_id"]]
)

dg_first.inputs.field_template = dict(
    struct = "sub-%s/anat/sub-%s_T1w.nii.gz",
    func="sub-%s/func/sub-%s_task-flanker_run-%d_bold.nii.gz",
    events="sub-%s/func/sub-%s_task-flanker_run-%d_events.tsv", 
)

wf_first.connect([
        (infosource_first, dg_first, [("subject_id", "subject_id")])
])

Initialisation#

Convert functional images to float representation. Since there can be more than one functional run we use a MapNode to convert each run.

Corresponding FSL command:

/usr/local/fsl/bin/fslmaths ../sub-11/func/sub-11_task-flanker_run-1_bold prefiltered_func_data -odt float
img2float = pe.MapNode(
    interface=fsl.ImageMaths(out_data_type='float', op_string='', suffix='_dtype'),
    name='img2float',
    iterfield=['in_file'])
wf_first.connect(dg_first, 'func', img2float, 'in_file')

Extract the middle volume of the first run as the reference

(Head movement, motion-correction)

Corresponding FSL command:

/usr/local/fsl/bin/fslroi prefiltered_func_data example_func 73 1
# Extract one roi
extract_ref = pe.MapNode(interface=fsl.ExtractROI(t_size=1), 
                         name='extractref',
                         iterfield=['in_file'])
wf_first.connect(img2float, 'out_file', extract_ref, 'in_file')

# Define a function to return the 1 based index of the middle volume
def getmiddlevolume(func):
    from nibabel import load
    funcfile = func
    if isinstance(func, list):
        funcfile = func[0]
    _, _, _, timepoints = load(funcfile).shape
    return int(timepoints/2) 

wf_first.connect(img2float, ('out_file', getmiddlevolume), extract_ref, 't_min')

Preprocessing#

Motion Correction#

Realign the functional runs to the middle volume of each run

Corresponding FSL command:

/usr/local/fsl/bin/mcflirt -in prefiltered_func_data -out prefiltered_func_data_mcf -mats -plots -reffile example_func -rmsrel -rmsabs -spline_final
save_mats (a boolean) – Save transformation matrices. Maps to a command-line argument: -mats.
save_plots (a boolean) – Save transformation parameters. Maps to a command-line argument: -plots.
save_rms (a boolean) – Save rms displacement parameters. Maps to a command-line argument: -rmsabs -rmsrel.
Interpolation (‘spline’ or ‘nn’ or ‘sinc’) – Interpolation method for transformation. Maps to a command-line argument: -%s_final.
motion_correct = pe.MapNode(
    interface=fsl.MCFLIRT(save_mats=True, save_plots=True, save_rms=True, interpolation='spline'),
    name='realign',
    iterfield=['in_file','ref_file'])
wf_first.connect(img2float, 'out_file', motion_correct, 'in_file')
wf_first.connect(extract_ref, 'roi_file', motion_correct, 'ref_file')

Plot the estimated motion parameters

Corresponding FSL command:

/usr/local/fsl/bin/fsl_tsplot -i prefiltered_func_data_mcf.par -t 'MCFLIRT estimated rotations (radians)' -u 1 --start=1 --finish=3 -a x,y,z -w 640 -h 144 -o rot.png 

/usr/local/fsl/bin/fsl_tsplot -i prefiltered_func_data_mcf.par -t 'MCFLIRT estimated translations (mm)' -u 1 --start=4 --finish=6 -a x,y,z -w 640 -h 144 -o trans.png 

/usr/local/fsl/bin/fsl_tsplot -i prefiltered_func_data_mcf_abs.rms,prefiltered_func_data_mcf_rel.rms -t 'MCFLIRT estimated mean displacement (mm)' -u 1 -w 640 -h 144 -a absolute,relative -o disp.png
plot_motion = pe.MapNode(
    interface=fsl.PlotMotionParams(in_source='fsl'),
    name='plot_motion',
    iterfield=['in_file'])
plot_motion.iterables = ('plot_type', ['rotations', 'translations','displacement'])
wf_first.connect(motion_correct, 'par_file', plot_motion, 'in_file')

Functionally Masking#

Passing the reference volume to the FSL command-line tool bet to generate a binary brain mask and afterward multiplying the processed functional time series by the brain mask using the fslmaths command to produce a skull-stripped time series.

See Ciric et al.(2018)

Extract the mean volume of each functional run

Corresponding FSL command:

/usr/local/fsl/bin/fslmaths prefiltered_func_data_mcf -Tmean mean_func
meanfunc = pe.MapNode(
    interface=fsl.ImageMaths(op_string='-Tmean', suffix='_mean'),
    name='meanfunc',
    iterfield=['in_file'])
wf_first.connect(motion_correct, 'out_file', meanfunc, 'in_file')

Strip the skull from the mean functional to generate a mask

Corresponding FSL command:

/usr/local/fsl/bin/bet2 mean_func mask -f 0.3 -n -m; /usr/local/fsl/bin/immv mask_mask mask
meanfuncmask = pe.MapNode(
    interface=fsl.BET(mask=True, no_output=True, frac=0.3),
    name='meanfuncmask',
    iterfield=['in_file'])
wf_first.connect(meanfunc, 'out_file', meanfuncmask, 'in_file')

Mask the functional data with the extracted mask

Corresponding FSL command:

/usr/local/fsl/bin/fslmaths prefiltered_func_data_mcf -mas mask prefiltered_func_data_bet
maskfunc = pe.MapNode(
    interface=fsl.ImageMaths(suffix='_bet', op_string='-mas'),
    name='maskfunc',
    iterfield=['in_file','in_file2'])
wf_first.connect(motion_correct, 'out_file', maskfunc, 'in_file')
wf_first.connect(meanfuncmask, 'mask_file', maskfunc, 'in_file2')

Grand Mean Scaling#

Determine the 2nd and 98th percentile intensities

Corresponding FSL command:

/usr/local/fsl/bin/fslstats prefiltered_func_data_bet -p 2 -p 98
0.000000 873.492249 (these numbers are for subject-11 run-01)

More info here

getthresh = pe.MapNode(
    interface=fsl.ImageStats(op_string='-p 2 -p 98'),
    name='getthreshold',
    iterfield=['in_file'])
wf_first.connect(maskfunc, 'out_file', getthresh, 'in_file')

Threshold the first TR of the functional data at 10% of the 98th percentile

Corresponding FSL command:

/usr/local/fsl/bin/fslmaths prefiltered_func_data_bet -thr 87.3492249 -Tmin -bin mask -odt char
threshold = pe.MapNode(
    interface=fsl.ImageMaths(out_data_type='char', suffix='_thresh'),
    name='threshold',
    iterfield=['in_file','op_string'])
wf_first.connect(maskfunc, 'out_file', threshold, 'in_file')

# Define a function to get 10% of the intensity
def getthreshop(thresh):
    return ['-thr %.10f -Tmin -bin' % (0.1 * th[1]) for th in thresh]

wf_first.connect(getthresh, ('out_stat', getthreshop), threshold, 'op_string')

Determine the median value of the TRs using the mask

Corresponding FSL command:

/usr/local/fsl/bin/fslstats prefiltered_func_data_mcf -k mask -p 50
728.800232 (this number is for subject-11 run-01)
medianval = pe.MapNode(
    interface=fsl.ImageStats(op_string='-k %s -p 50'),
    name='medianval',
    iterfield=['in_file','mask_file'])
wf_first.connect(motion_correct, 'out_file', medianval, 'in_file')
wf_first.connect(threshold, 'out_file', medianval, 'mask_file')

Dilate the mask

The brain mask is “dilated” slightly before being used. Because it is normally important that masking be liberal (ie that there be little risk of cutting out valid brain voxels)

Corresponding FSL command:

/usr/local/fsl/bin/fslmaths mask -dilF mask

The output of dilatemask (i.e., /srv/scratch/yc/fsl/hw2/level1/_subject_id_26/dilatemask/mapflow/_dilatemask0/sub-26_task-flanker_run-1_bold_dtype_mcf_bet_thresh_dil.nii.gz) is equivalent to mask.nii.gz from FSL GUI.

dilatemask = pe.MapNode(
    interface=fsl.ImageMaths(suffix='_dil', op_string='-dilF'),
    name='dilatemask',
    iterfield=['in_file'])
wf_first.connect(threshold, 'out_file', dilatemask, 'in_file')

Mask the motion corrected functional runs with the dilated mask

Corresponding FSL command:

/usr/local/fsl/bin/fslmaths prefiltered_func_data_mcf -mas mask prefiltered_func_data_thresh
maskfunc2 = pe.MapNode(
    interface=fsl.ImageMaths(suffix='_thresh', op_string='-mas'),
    iterfield=['in_file','in_file2'],
    name='maskfunc2')
wf_first.connect(motion_correct, 'out_file', maskfunc2, 'in_file')
wf_first.connect(dilatemask, 'out_file', maskfunc2, 'in_file2')

SUSAN Noise Reduction#

Determine the mean image from each TR

Corresponding FSL command:

/usr/local/fsl/bin/fslmaths prefiltered_func_data_thresh -Tmean mean_func
meanfunc2 = pe.MapNode(
    interface=fsl.ImageMaths(op_string='-Tmean', suffix='_mean'),
    name='meanfunc2',
    iterfield=['in_file'])
wf_first.connect(maskfunc2, 'out_file', meanfunc2, 'in_file')

Merge the median values with the mean functional images into a coupled list

The output of this merge node will go into susan as usans,

# why here use node not mapnode?
mergenode = pe.Node(interface=util.Merge(2, axis='hstack'), 
                       name='merge')
wf_first.connect(meanfunc2, 'out_file', mergenode, 'in1')
wf_first.connect(medianval, 'out_stat', mergenode, 'in2')

Smooth each run using SUSAN with the brightness threshold set to 75% of the median value for each run and a mask constituting the mean functional

Usage: susan <input> <bt> <dt> <dim> <use_median> <n_usans> [<usan1> <bt1> [<usan2> <bt2>]] <output>
<bt> is brightness threshold and should be greater than noise level and less than contrast of edges to be preserved.
<dt> is spatial size (sigma, i.e., half-width) of smoothing, in mm.
<dim> is dimensionality (2 or 3), depending on whether smoothing is to be within-plane (2) or fully 3D (3).
<use_median> determines whether to use a local median filter in the cases where single-point noise is detected (0 or 1).
<n_usans> determines whether the smoothing area (USAN) is to be found from secondary images (0, 1 or 2).
A negative value for any brightness threshold will auto-set the threshold at 10% of the robust range

Corresponding FSL command:

/usr/local/fsl/bin/susan prefiltered_func_data_thresh 546.600174 2.12314225053 3 1 1 mean_func 546.600174 prefiltered_func_data_smooth

Note:

for <bt>, Nipype uses a different algorithm to calculate it -> float(fwhm) / np.sqrt(8 * np.log(2)). Therefore, to get 2.12314225053, fwhm should be 4.9996179300001655 instead of 5

fwhm_thr = 4.9996179300001655

smooth = pe.MapNode(
    interface=fsl.SUSAN(fwhm = fwhm_thr),
    name='smooth',
    iterfield=['in_file', 'brightness_threshold', 'usans'])

# Define a function to get the brightness threshold for SUSAN
def getbtthresh(medianvals):
    return [0.75 * val for val in medianvals]

def getusans(x):
    return [[tuple([val[0], 0.75 * val[1]])] for val in x]

wf_first.connect(maskfunc2, 'out_file', smooth, 'in_file')
wf_first.connect(medianval, ('out_stat', getbtthresh), smooth,
                'brightness_threshold')
wf_first.connect(mergenode, ('out', getusans), smooth, 'usans')

Mask the smoothed data with the dilated mask

Corresponding FSL command:

/usr/local/fsl/bin/fslmaths prefiltered_func_data_smooth -mas mask prefiltered_func_data_smooth
maskfunc3 = pe.MapNode(
    interface=fsl.ImageMaths(op_string='-mas'),
    name='maskfunc3',
    iterfield=['in_file','in_file2'])
wf_first.connect(smooth, 'smoothed_file', maskfunc3, 'in_file')
wf_first.connect(dilatemask, 'out_file', maskfunc3, 'in_file2')

Scale each volume of the TR so that the median value of the TR is set to 10000

Corresponding FSL command:

/usr/local/fsl/bin/fslmaths prefiltered_func_data_smooth -mul 13.7211811425 prefiltered_func_data_intnorm 
(this number is for subject-11 run-01)
intnorm = pe.MapNode(
    interface=fsl.ImageMaths(suffix='_intnorm'),
    iterfield=['in_file', 'op_string'],
    name='intnorm')
wf_first.connect(maskfunc3, 'out_file', intnorm, 'in_file')

# Define a function to get the scaling factor for intensity normalization
def getinormscale(medianvals):
    return ['-mul %.10f' % (10000. / val) for val in medianvals]

wf_first.connect(medianval, ('out_stat', getinormscale), intnorm, 'op_string')

Temporal Filtering#

Perform temporal highpass filtering on the data

Corresponding FSL command:

/usr/local/fsl/bin/fslmaths prefiltered_func_data_intnorm -Tmean tempMean
/usr/local/fsl/bin/fslmaths prefiltered_func_data_intnorm -bptf 25.0 -1 -add tempMean prefiltered_func_data_tempfilt

The output of highpass (i.e., sub-11_task-flanker_run-1_bold_dtype_mcf_mask_smooth_mask_intnorm_tempfilt.nii.gz) is equivalent to filtered_func_data.nii.gz from FSL GUI.

# Generate a mean functional image from the scaled data and this mean func will be used in performing temporal highpass filtering
meanfunc3 = pe.MapNode(
    interface=fsl.ImageMaths(op_string='-Tmean'),
    iterfield=['in_file'],
    name='meanfunc3')
wf_first.connect(intnorm, 'out_file', meanfunc3, 'in_file')
# Perform temporal highpass filtering on the data
hpcutoff = 100
TR = 2.  # ensure float

highpass = pe.MapNode(
    interface=fsl.ImageMaths(suffix='_tempfilt'),
    name='highpass',
    iterfield=['in_file','op_string'])

# 25 = (hpcutoff / 2*TR) not (hpcutoff / TR)
def gethpstring(tempMean):
    return ['-bptf 25 -1 -add %s' % (tm) for tm in tempMean]

wf_first.connect(intnorm, 'out_file', highpass, 'in_file')
wf_first.connect(meanfunc3, ('out_file',gethpstring), highpass, 'op_string')

Generate a mean functional image from the functional run

Corresponding FSL command:

/usr/local/fsl/bin/fslmaths prefiltered_func_data_tempfilt filtered_func_data
/usr/local/fsl/bin/fslmaths filtered_func_data -Tmean mean_func
meanfunc4 = pe.MapNode(
    interface=fsl.ImageMaths(op_string='-Tmean', suffix='_mean'),
    iterfield=['in_file'],
    name='meanfunc4')
wf_first.connect(highpass, 'out_file', meanfunc4, 'in_file')

First-Level GLM#

Get events#

def subjinfo(events):
    from nipype.interfaces.base import Bunch
    import pandas as pd
    import numpy as np

    subject_info = []

    ev = pd.read_csv(events, sep="\t")
    ev = ev[ev['correctness']=='correct']
    ev['new_type'] = ev['trial_type'].apply(lambda x: str(x).split('_')[0])

    run_info = Bunch(onsets=[], 
                     durations=[])

    run_info.set(conditions=[g[0] for g in ev.groupby("new_type")])

    for group in ev.groupby("new_type"):
        run_info.onsets.append(group[1].onset.tolist())
        run_info.durations.append(group[1].duration.tolist())
    subject_info.append(run_info)

    return subject_info

get_sub_info = pe.MapNode(
    Function(
        function=subjinfo, input_names=["events"], output_names="subject_info"
    ),
    name="get_sub_info", iterfield=["events"]
)

# Connect to workflow
wf_first.connect(dg_first, 'events', get_sub_info, 'events')

Set-up contrasts#

condition_names = ["congruent", "incongruent"]

# Activation Contrasts (similar to https://direct.mit.edu/jocn/article/23/10/3162/5327/Is-Morality-Unified-Evidence-that-Distinct-Neural)
## From FEAT: "The correct way to tell whether two conditions are significantly different is to run a differential contrast like [1 -1] between them"
cont01 = ("congruent", "T", condition_names,  [1, 0])
cont02 = ("incongruent", "T", condition_names, [0, 1])
cont03 = ("congruent-incongruent", "T", condition_names, [1, -1])
cont04 = ("incongruent-congruent", "T", condition_names,  [-1, 1])

contrast_list = [
    cont01,
    cont02,
    cont03,
    cont04
]
TR = 2.0  # Repetition Time
hpcutoff = 100  # highpass filter cutoff in seconds! 

modelspec = pe.MapNode(interface=model.SpecifyModel(), name="modelspec", iterfield=["subject_info","functional_runs"])
modelspec.inputs.input_units = "secs"
modelspec.inputs.high_pass_filter_cutoff = hpcutoff
modelspec.inputs.time_repetition = TR

wf_first.connect(get_sub_info, 'subject_info', modelspec,'subject_info')
wf_first.connect(highpass, 'out_file',modelspec,'functional_runs')
level1design = pe.MapNode(interface=fsl.Level1Design(), name="level1design", iterfield=["session_info"])
level1design.inputs.interscan_interval = TR
# Set HRF bases functions
level1design.inputs.bases = {"gamma": {"derivs": False, 'gammasigma':3, 'gammadelay':6}}
level1design.inputs.model_serial_correlations = True
level1design.inputs.contrasts = contrast_list

wf_first.connect(modelspec,'session_info', level1design, 'session_info')
modelgen = pe.MapNode(
    interface=fsl.FEATModel(),
    name='modelgen',
    iterfield=['fsf_file', 'ev_files'])

wf_first.connect([(level1design, modelgen, [('fsf_files', 'fsf_file'), 
                                    ('ev_files','ev_files')])])

Corresponding fsl command:

/usr/local/fsl/bin/film_gls --in=filtered_func_data --rn=stats --pd=design.mat --thr=1000.0 --sa --ms=5 --con=design.con  

--thr: threshold
--sa: smooth_autocorr
--ms: mask_size
level1estimate = pe.MapNode(
    interface=fsl.FILMGLS(smooth_autocorr=True, mask_size=5, threshold=1000),
    name='level1estimate',
    iterfield=['design_file', 'in_file', 'tcon_file'])

wf_first.connect([
    (highpass,level1estimate,[('out_file','in_file')]),
    (modelgen,level1estimate,[('design_file','design_file'),
                               ('con_file','tcon_file')])
])

no need to use ContrastMgr,

In interface mode this file assumes that all the required inputs are in the same location. This has deprecated for FSL versions 5.0.7+ as the necessary corrections file is no longer generated by FILMGLS.

see this link

Registration#

According to FSL UserGuide, the different sessions need to be registered to each other before any multi-session or multi-subject analyses can be carried out. Registration inside FEAT uses FLIRT and is a two-statge process:

  1. An example FMRI low resolution image is registered to an example high resolution image (normally the same subject’s T1-weighted structural). The transformation for this is saved into the FEAT directory. Then the high res image is registered to a standard image (normally a T1-weighted image in standard space, such as the MNI 152 average image).

  2. The two transformations are combined into a third, which will take the low resolution FMRI images (and the statistic images derived from the first-level analyses) straight into standard space, when applied later, during group analysis.

Step 1#

Corresponding FSL command:

/usr/local/fsl/bin/flirt -in example_func -ref standard -out example_func2standard -omat example_func2standard.mat -cost corratio -dof 12 -searchrx -90 90 -searchry -90 90 -searchrz -90 90 -interp trilinear 

/usr/local/fsl/bin/convert_xfm -inverse -omat standard2example_func.mat example_func2standard.mat
flt = pe.MapNode(interface=fsl.FLIRT(cost='corratio', dof=12,
                                    searchr_x = [-90,90],
                                    searchr_y = [-90,90],
                                    searchr_z = [-90,90],
                                    interp = 'trilinear'), 
               name='example_func2standard',
               iterfield=['in_file'])

flt.inputs.reference = fsl.Info.standard_image('MNI152_T1_2mm_brain.nii.gz')
wf_first.connect(extract_ref, 'roi_file', flt, 'in_file')

convertxfm = pe.MapNode(interface=fsl.ConvertXFM(invert_xfm = True), 
                         name='convertxfm',
                         iterfield=['in_file'])
wf_first.connect(flt, 'out_matrix_file', convertxfm, 'in_file')

Step 2#

Corresponding fsl command:

/usr/local/fsl/bin/flirt -ref reg/standard -in stats/cope1 -out /srv/scratch/yc/fsl/nipype_fsl_comp/gui/sub-01/run1.feat/frgrot_chksugax -applyxfm -init reg/example_func2standard.mat -interp trilinear -datatype float
def warp_files(copes, varcopes, masks, mat):
    # need to reimport here, otherwise errors come out
    import nipype.interfaces.fsl as fsl 
    
    out_copes = []
    out_varcopes = []
    out_masks = []
    
    # register mask, same function, different parameters
    warp_mask = fsl.FLIRT(apply_xfm = True, 
                     interp = 'nearestneighbour')
    warp_mask.inputs.reference = fsl.Info.standard_image('MNI152_T1_2mm_brain.nii.gz')
    warp_mask.inputs.in_matrix_file = mat
    warp_mask.inputs.output_type = "NIFTI_GZ"
    warp_mask.inputs.in_file = masks
    res_mask = warp_mask.run()
    out_masks.append(str(res_mask.outputs.out_file))
    
    # register copes & varcopes using same function, different parameters
    warp = fsl.FLIRT(apply_xfm = True, 
                     interp = 'trilinear')
    warp.inputs.reference = fsl.Info.standard_image('MNI152_T1_2mm_brain.nii.gz')
    warp.inputs.in_matrix_file = mat
    warp.inputs.output_type = "NIFTI_GZ"
    
    # register copes
    for cope in copes:
        warp.inputs.in_file = cope
        res = warp.run()
        out_copes.append(str(res.outputs.out_file))
        
     # register varcopes
    for varcope in varcopes:
        warp.inputs.in_file = varcope
        res = warp.run()
        out_varcopes.append(str(res.outputs.out_file))

    return out_copes, out_varcopes, out_masks
warpfunc = pe.MapNode(util.Function(input_names=['copes', 'varcopes', 'masks', 'mat'],
                                    output_names=['out_copes', 'out_varcopes', 'out_masks'],
                                    function=warp_files),
                      iterfield=['copes', 'varcopes', 'masks','mat'],
                      name='warpfunc')

wf_first.connect(level1estimate, 'copes', warpfunc, 'copes')
wf_first.connect(level1estimate, 'varcopes', warpfunc, 'varcopes')
wf_first.connect(dilatemask, 'out_file', warpfunc, 'masks')
wf_first.connect(flt, 'out_matrix_file', warpfunc, 'mat')
# save all results into single/separate folder for further analysis use
datasink_first = pe.Node(nio.DataSink(), name='sinker')
datasink_first.inputs.base_directory=os.path.join(exp_dir_first, "level1_results")

wf_first.connect(infosource_first, 'subject_id', datasink_first, 'container')
wf_first.connect([(level1estimate, datasink_first, [('results_dir', 'results_dir')])])
wf_first.connect([(warpfunc, datasink_first, [('out_copes', 'reg_copes')])])
wf_first.connect([(warpfunc, datasink_first, [('out_varcopes', 'reg_varcopes')])])
wf_first.connect([(warpfunc, datasink_first, [('out_masks', 'reg_masks')])])
# Create 1st-level analysis output graph
wf_first.write_graph(graph2use='colored', format='png', simple_form=True)

# Visualize the graph
from IPython.display import Image
Image(filename=os.path.join(wf_first.base_dir, wf_first.name, 'graph.png'))
230622-04:04:13,752 nipype.workflow INFO:
	 Generated workflow graph: /home/jovyan/example-notebooks/functional_imaging/output_level1/level1/graph.png (graph2use=colored, simple_form=True).
../_images/a062279807265ff39c961ecbe84a428502dcece4ef8ede59a95af6b44c5290a8.png
import logging
logging.getLogger('nipype.workflow').setLevel(0) #set to 1 if nipype output is wanted
# Run Workflow
wf_first.run(plugin="MultiProc", plugin_args={"n_procs": 8})
<networkx.classes.digraph.DiGraph at 0x7f21f9399810>

The first-level GLM#

Here we randomly choose the four copes from subject-09 run-1

fsaverage = nilearn.datasets.fetch_surf_fsaverage()
plotting.show() 
nipype_cope1 = './output_level1/level1_results/09/results_dir/_subject_id_09/_level1estimate0/results/cope1.nii.gz'
nipype_cope2 = './output_level1/level1_results/09/results_dir/_subject_id_09/_level1estimate0/results/cope2.nii.gz'
nipype_cope3 = './output_level1/level1_results/09/results_dir/_subject_id_09/_level1estimate0/results/cope3.nii.gz'
nipype_cope4 = './output_level1/level1_results/09/results_dir/_subject_id_09/_level1estimate0/results/cope4.nii.gz'
plotting.plot_stat_map(nipype_cope1, bg_img=nipype_cope1, title = 'Nipype - COPE1', cmap = 'bwr', colorbar = False)
plt.show()
plotting.plot_stat_map(nipype_cope2, bg_img=nipype_cope1, title = 'Nipype - COPE2', cmap = 'bwr', colorbar = False)
plt.show()
plotting.plot_stat_map(nipype_cope3, bg_img=nipype_cope1, title = 'Nipype - COPE3', cmap = 'bwr', colorbar = False)
plt.show()
plotting.plot_stat_map(nipype_cope4, bg_img=nipype_cope1, title = 'Nipype - COPE4', cmap = 'bwr', colorbar = False)
plt.show()
../_images/931ab0f4c71eda49ae5e6167a2672dd59539a5e9f6a9fbe1ed2d8c757f724a8e.png ../_images/ca87416da6220b551ec0094ece084427f186273593d763cba989ddd2a9eec6e3.png ../_images/966f776d79bec4829951b87747c9d1785ba05de3fb87b91dd1c897d4e75a6b93.png ../_images/c656f64454b3cf78de946e554d846c9985e937a5333556973768d421a209572e.png

Second level GLM using Nipype FSL#

wf_second = pe.Workflow(name='level2', base_dir=exp_dir_second)
wf_second.config["execution"]["crashfile_format"] = "txt"

The following two nodes (infosource & dg) together define all inputs required for the preprocessing workflow

# we want to group the outcome by contrast not subject
contr_list = [1,2,3,4]
infosource_second = pe.Node(util.IdentityInterface(fields=["contr_id"]),
                  name="infosource_second")
infosource_second.iterables = [("contr_id", contr_list)]
# here we use SelectFiles, instead of DataGrabber, because the former is more flexible with formatting syntax
templates_second = {
    "reg_copes":"*/reg_copes/*/*/cope{contr_id}_flirt.nii.gz",
    "reg_varcopes":"*/reg_varcopes/*/*/varcope{contr_id}_flirt.nii.gz",
    "reg_masks":"*/reg_masks/*/*/*.nii.gz"
}
dg_second = pe.Node(interface=nio.SelectFiles(templates_second),
             name="selectfiles_second")
dg_second.inputs.base_directory = data_dir_second

wf_second.connect([
        (infosource_second, dg_second, [("contr_id", "contr_id")])
])

Second-level GLM#

Combining results from multiple runs of one subject into one

Higher-level input files preparation#

Step 1: Merge registered copes & varcopes & masks#

Corresponding FSL command:

/usr/local/fsl/bin/fslmerge -t mask (masks from all 52 inputs)
/usr/local/fsl/bin/fslmerge -t cope (copes from all 52 inputs)
/usr/local/fsl/bin/fslmerge -t varcop (varcopes from all 52 inputs)
copemerge_second = pe.Node(
    interface=fsl.Merge(dimension='t'),
    iterfield=['in_files'],
    name="copemerge_second")

varcopemerge_second = pe.Node(
    interface=fsl.Merge(dimension='t'),
    iterfield=['in_files'],
    name="varcopemerge_second")

maskmerge_second = pe.Node(
    interface=fsl.Merge(dimension='t'),
    iterfield=['in_files'],
    name="maskmerge_second")


def sort_copes(files):
    numelements = len(files[0])
    outfiles = []
    for i in range(numelements):
        outfiles.insert(i, [])
        for j, elements in enumerate(files):
            outfiles[i].append(elements[i])
    return outfiles

wf_second.connect(dg_second, 'reg_varcopes', varcopemerge_second, 'in_files')
wf_second.connect(dg_second, 'reg_copes', copemerge_second, 'in_files')
wf_second.connect(dg_second, 'reg_masks', maskmerge_second, 'in_files')

Step 2: Making mask#

In FSL, there are many commands about maskunique, which is unless for the second level. We can ignore it.

Corresponding FSL command:

/usr/local/fsl/bin/fslmaths mask -Tmin mask
minmask_second = pe.Node(
    interface=fsl.ImageMaths(op_string='-Tmin'),
    iterfield=['in_file'],
    name='minmask_second')

wf_second.connect(maskmerge_second, 'merged_file', minmask_second, 'in_file')

Step 3: Masking copes & varcopes#

Corresponding FSL command:

we have four contrasts so the following commands repeat four times

/usr/local/fsl/bin/fslmaths cope1 -mas mask cope1
/usr/local/fsl/bin/fslmaths varcope1 -mas mask varcope1
maskcope_second = pe.Node(
    interface=fsl.ImageMaths(op_string='-mas'),
    iterfield=['in_file', 'in_file2'],
    name='maskcope_second')

maskvarcope_second = pe.Node(
    interface=fsl.ImageMaths(op_string='-mas'),
    iterfield=['in_file', 'in_file2'],
    name='maskvarcope_second')

wf_second.connect(copemerge_second, 'merged_file', maskcope_second, 'in_file')
wf_second.connect(minmask_second, 'out_file', maskcope_second, 'in_file2')
wf_second.connect(varcopemerge_second, 'merged_file', maskvarcope_second, 'in_file')
wf_second.connect(minmask_second, 'out_file', maskvarcope_second, 'in_file2')

Set up second-level contrasts and fixed-effects#

def get_contrasts_l2(in_files):
    import numpy as np
    import glob
    total = len(in_files)
    print(in_files)
    print(total)
    n_sub = 9  
    # n_sub = 26
    ev_list = ['ev'+str(x) for x in range(1,n_sub+1)]
    weight_mtx = np.zeros((n_sub,n_sub))
    weight_mtx = weight_mtx.astype(np.float64)
    np.fill_diagonal(weight_mtx,1.)
    #contr = ['','T',ev_list,list(weight_mtx[0])] --> original paper
    contr = np.array(['','T', ev_list, list(weight_mtx[0])], dtype=object) #adaptation so the construction of the array works
    contr_lst = np.tile(contr, (n_sub,n_sub))
    contr_lst = [tuple(x) for x in contr_lst] 
    for i in range(n_sub):
        # Instead of trying to change the tuple, create a new tuple with the changed value.
        contr_lst[i] = (contr_lst[i][0], contr_lst[i][1], contr_lst[i][2], list(weight_mtx[i]))

    reg_dict = {k:None for k in ev_list}
    for k in reg_dict.keys():
        start_lst = [0.0] * total
        idx = ev_list.index(k)
        start_idx = idx*2
        end_idx = idx*2 + 1
        start_lst[start_idx] = 1.
        start_lst[end_idx] = 1.
        reg_dict[k] = start_lst    
    return contr_lst, reg_dict

contrastgen_l2 = pe.Node(util.Function(input_names=['in_files'],
                                    output_names=['contr_lst', 'reg_dict'],
                                    function=get_contrasts_l2),
                      iterfield=['in_files'],
                      name='contrastgen_l2')

wf_second.connect(dg_second,'reg_copes', contrastgen_l2, 'in_files')

Nipype recommends using L2Model, which only works for the single subject. It takes the number of runs (copes at the first level) as input and does estimations for subject one by one. Instead, we use MultipleRegressDesign. As it’s name indicates, this one can deal with multiple predictors (subjects) at the same time.

level2model = pe.Node(interface=fsl.MultipleRegressDesign(),
                      name='l2model')

wf_second.connect([(contrastgen_l2, level2model, [('contr_lst','contrasts'),
                                        ('reg_dict','regressors')])])

level2estimate = pe.Node(
    interface=fsl.FLAMEO(run_mode='fe'),
    name="level2estimate",
    iterfield=['cope_file', 'var_cope_file'])

wf_second.connect([
    (maskcope_second, level2estimate, [('out_file', 'cope_file')]),
    (maskvarcope_second, level2estimate, [('out_file', 'var_cope_file')]),
    (minmask_second, level2estimate, [('out_file', 'mask_file')]),
    (level2model, level2estimate, [('design_mat', 'design_file'),
                                   ('design_con', 't_con_file'), 
                                   ('design_grp', 'cov_split_file')]),
])

datasink_second = pe.Node(nio.DataSink(), name='sinker')
datasink_second.inputs.base_directory=os.path.join(exp_dir_second, "level2_results")

int2string = lambda x: 'contrast_'+str(x)
    
wf_second.connect(infosource_second, ('contr_id',int2string), datasink_second, 'container')
wf_second.connect([(level2estimate, datasink_second, [('stats_dir', 'stats_dir')])])

# Create 1st-level analysis output graph
wf_second.write_graph(graph2use='colored', format='png', simple_form=True)

# Visualize the graph
from IPython.display import Image
Image(filename=os.path.join(wf_second.base_dir, wf_second.name, 'graph.png'))
../_images/31f70591725b61c51f29bfe0b898025ee6fe0b079f5a48028dcbed27f371e8c8.png
# Run Workflow
wf_second.run(plugin="MultiProc", plugin_args={"n_procs": 8})
<networkx.classes.digraph.DiGraph at 0x7f21dd61b130>

The second-level GLM#

Here we choose the four copes from subject-01

Regapply#

This happens before the second-level GLM, all results from the first level are needed to be registered into the standard space

nipype_cope1 = './output_level1/level1_results/01/reg_copes/_subject_id_01/_warpfunc1/cope1_flirt.nii.gz'
nipype_cope2 = './output_level1/level1_results/01/reg_copes/_subject_id_01/_warpfunc1/cope2_flirt.nii.gz'
nipype_cope3 = './output_level1/level1_results/01/reg_copes/_subject_id_01/_warpfunc1/cope3_flirt.nii.gz'
nipype_cope4 = './output_level1/level1_results/01/reg_copes/_subject_id_01/_warpfunc1/cope4_flirt.nii.gz'
plotting.plot_stat_map(nipype_cope1, bg_img=nipype_cope1, title = 'Nipype - COPE1', cmap = 'bwr', colorbar = False)
plotting.plot_stat_map(nipype_cope2, bg_img=nipype_cope2, title = 'Nipype - COPE2', cmap = 'bwr', colorbar = False)
plotting.plot_stat_map(nipype_cope3, bg_img=nipype_cope3, title = 'Nipype - COPE3', cmap = 'bwr', colorbar = False)
plotting.plot_stat_map(nipype_cope4, bg_img=nipype_cope4, title = 'Nipype - COPE4', cmap = 'bwr', colorbar = False)
plt.show()
../_images/15240f937ff62d363ff0d2f0038aebf84d60e8551bf63e704cde17ba82e39145.png ../_images/dbf4b2ce846953767b9826a8f70fc6fe013b94f24c97cb33dffad322ba2d75e8.png ../_images/9bec52a10b3f6b83cf56f5cc5110b80ce1c4eb00ecb3abc21af7b7387a179654.png ../_images/d79995b6da6a54db5754d771c5f1a34a7f5d228ba507343c6a124d0650231715.png

Copes from the second-level GLM#

copes from subject-01

nipype_cope1 = './output_level2/level2_results/contrast_1/stats_dir/_contr_id_1/stats/cope1.nii.gz'
nipype_cope2 = './output_level2/level2_results/contrast_2/stats_dir/_contr_id_2/stats/cope1.nii.gz'
nipype_cope3 = './output_level2/level2_results/contrast_3/stats_dir/_contr_id_3/stats/cope1.nii.gz'
nipype_cope4 = './output_level2/level2_results/contrast_4/stats_dir/_contr_id_4/stats/cope1.nii.gz'
plotting.plot_glass_brain(nipype_cope1, bg_img=nipype_cope1, title = 'Nipype - COPE1')
plotting.plot_glass_brain(nipype_cope2, bg_img=nipype_cope2, title = 'Nipype - COPE2')
plotting.plot_glass_brain(nipype_cope3, bg_img=nipype_cope3, title = 'Nipype - COPE3')
plotting.plot_glass_brain(nipype_cope4, bg_img=nipype_cope4, title = 'Nipype - COPE4')
plt.show()
../_images/c2cc4b07dd32b076d70696f9bd3f5176241f60e653777dcb51941e7fea47b06e.png ../_images/428afe70713871101005cdeececcdbbe7db8120a6ea622ec7a52fca1359f1cea.png ../_images/5223962992efee919878975c9a19ed881113629466f73bee5b4231c671aa0e2b.png ../_images/1113b0104e03314bf24b341821e092f139d046be77d7aba255734292741b66b2.png

Third-level GLM using Nipype FSL#

Start the workflow

wf_third = pe.Workflow(name='level3', base_dir=exp_dir_third)
wf_third.config["execution"]["crashfile_format"] = "txt"

The following two nodes (infosource & dg) together define all inputs required for the preprocessing workflow

# we want to group the outcome by contrast not subject
contr_list = [1,2,3,4]
infosource_third = pe.Node(util.IdentityInterface(fields=["contr_id"]),
                  name="infosource_third")
infosource_third.iterables = [("contr_id", contr_list)]
# here we use SelectFiles, instead of DataGrabber, because the former is more flexible with formatting syntax
templates_third = {
    "copes":"contrast_{contr_id}/*/_contr_id_{contr_id}/*/cope*.nii.gz",
    "varcopes":"contrast_{contr_id}/*/_contr_id_{contr_id}/*/varcope*.nii.gz",
    "masks":"contrast_{contr_id}/*/_contr_id_{contr_id}/*/mask.nii.gz"
}
templates = {
    "reg_copes":"*/reg_copes/*/*/cope{contr_id}_flirt.nii.gz",
    "reg_varcopes":"*/reg_varcopes/*/*/varcope{contr_id}_flirt.nii.gz",
    "reg_masks":"*/reg_masks/*/*/*.nii.gz"
}

dg_third = pe.Node(interface=nio.SelectFiles(templates_third),
             name="selectfiles_third")
dg_third.inputs.base_directory = data_dir_third

wf_third.connect([
        (infosource_third, dg_third, [("contr_id", "contr_id")])
])

Third-level GLM#

Combining results from multiple runs of one subject into one

Higher-level input files preparation#

Step 1: Merge registered copes & varcopes & masks#

Corresponding FSL command:

/usr/local/fsl/bin/fslmerge -t mask (masks from all 26 inputs)
/usr/local/fsl/bin/fslmerge -t cope (copes from all 26 inputs)
/usr/local/fsl/bin/fslmerge -t varcop (varcopes from all 26 inputs)
copemerge_third = pe.Node(
    interface=fsl.Merge(dimension='t'),
    iterfield=['in_files'],
    name="copemerge_third")

varcopemerge_third = pe.Node(
    interface=fsl.Merge(dimension='t'),
    iterfield=['in_files'],
    name="varcopemerge_third")

maskmerge_third = pe.Node(
    interface=fsl.Merge(dimension='t'),
    iterfield=['in_files'],
    name="maskmerge_third")

def repeat_mask(file):
    import numpy as np
    import glob
    # n_sub = len(glob.glob(os.path.join(os.getcwd(), 'ds000102')+"/sub-0*"))    
    n_sub = 9
    mask_lst = [file]
    repeated = np.repeat(mask_lst,n_sub)
    return list(repeated)

wf_third.connect(dg_third, 'copes', copemerge_third, 'in_files')
wf_third.connect(dg_third, 'varcopes', varcopemerge_third, 'in_files')
wf_third.connect(dg_third, ('masks', repeat_mask), maskmerge_third, 'in_files')

Step 2: Making mask#

In FSL, there are many commands about maskunique, which is unless for the second level. We can ignore it.

Corresponding FSL command:

/usr/local/fsl/bin/fslmaths mask -Tmin mask
# /usr/local/fsl/bin/fslmaths mask -Tmin mask
minmask_third = pe.Node(
    interface=fsl.ImageMaths(op_string='-Tmin'),
    iterfield=['in_file'],
    name='minmask_third')

wf_third.connect(maskmerge_third, 'merged_file', minmask_third, 'in_file')

Step 3: Masking copes & varcopes#

Corresponding FSL command:

we have four contrasts so the following commands repeat four times

/usr/local/fsl/bin/fslmaths cope1 -mas mask cope1
/usr/local/fsl/bin/fslmaths varcope1 -mas mask varcope1
maskcope_third = pe.Node(
    interface=fsl.ImageMaths(op_string='-mas'),
    iterfield=['in_file', 'in_file2'],
    name='maskcope_third')

maskvarcope_third = pe.Node(
    interface=fsl.ImageMaths(op_string='-mas'),
    iterfield=['in_file', 'in_file2'],
    name='maskvarcope_third')

wf_third.connect(copemerge_third, 'merged_file', maskcope_third, 'in_file')
wf_third.connect(minmask_third, 'out_file', maskcope_third, 'in_file2')
wf_third.connect(varcopemerge_third, 'merged_file', maskvarcope_third, 'in_file')
wf_third.connect(minmask_third, 'out_file', maskvarcope_third, 'in_file2')

Set up third-level contrasts#

Since we only have a single-group set-up, we can actually use L2Model. If we have an ANOVA-like design, using MultipleRegressDesign would be a better option.

def num_copes(files):
    return len(files)

level3model = pe.Node(interface=fsl.L2Model(), name='l3model')
wf_third.connect(dg_third, ('copes',num_copes), level3model, 'num_copes')
level3estimate = pe.Node(
    interface=fsl.FLAMEO(run_mode='flame1'),
    name="level3estimate",
    iterfield=['cope_file', 'var_cope_file'])

wf_third.connect([
    (maskcope_third, level3estimate, [('out_file', 'cope_file')]),
    (maskvarcope_third, level3estimate, [('out_file', 'var_cope_file')]),
    (minmask_third, level3estimate, [('out_file', 'mask_file')]),
    (level3model, level3estimate, [('design_mat', 'design_file'),
                                   ('design_con', 't_con_file'), 
                                   ('design_grp', 'cov_split_file')]),
])

Post-Stats#

Smoothness estimation#

to get dlh and volume for thresholding Corresponding FSL command:

/usr/local/fsl/bin/smoothest -d 25 -m mask -r stats/res4d > stats/smoothness
smoothest = pe.Node(
    interface=fsl.SmoothEstimate(dof=25),
    name="smoothest",
    iterfield=['residual_fit_file', 'mask_file'])

wf_third.connect(minmask_third, 'out_file', smoothest, 'mask_file')
wf_third.connect(level3estimate, 'res4d', smoothest, 'residual_fit_file')

Mask zstats file#

preparation for thresholding

Corresponding FSL command:

/usr/local/fsl/bin/fslmaths stats/zstat1 -mas mask thresh_zstat1
level3mask = pe.Node(
    interface=fsl.ImageMaths(op_string='-mas'),
    iterfield=['in_file','in_file2'],
    name='level3mask')

wf_third.connect([
    (level3estimate,level3mask,[('zstats', 'in_file')]),
    (minmask_third, level3mask, [('out_file', 'in_file2')]),
])

Cluster-wise thresholding#

Corresponding FSL command:

/usr/local/fsl/bin/cluster -i thresh_zstat1 -t 3.1 --othresh=thresh_zstat1 -o cluster_mask_zstat1 --connectivity=26 --mm --olmax=lmax_zstat1_std.txt --scalarname=Z -p 0.05 -d 0.0595781 --volume=254734 -c stats/cope1 > cluster_zstat1_std.txt
level3threshold = pe.Node(
        interface=fsl.Cluster(threshold = 3.1,
                              pthreshold = 0.05,
                              use_mm = True,
                              out_threshold_file =True,
                              out_index_file = True,
                              out_localmax_txt_file = True),
    iterfield=['in_file','cope_file'],
    name='level3threshold')

wf_third.connect([
    (level3mask, level3threshold, [('out_file', 'in_file')]),
    (level3estimate, level3threshold, [('copes', 'cope_file')]),
    (smoothest, level3threshold, [('dlh', 'dlh'),
                                   ('volume', 'volume')]),
])

Save the output#

datasink_third = pe.Node(nio.DataSink(), name='sinker')
datasink_third.inputs.base_directory=os.path.join(exp_dir_third, "level3_results")

int2string = lambda x: 'contrast_'+str(x)
    
wf_third.connect(infosource_third, ('contr_id',int2string), datasink_third, 'container')
wf_third.connect([
    (level3estimate, datasink_third, [('stats_dir', 'stats_dir')]),
    (level3threshold, datasink_third, [('threshold_file', 'thresholded'),
                                ('localmax_txt_file', 'localmax_txt'),
                                ('index_file', 'index')]),
           ])
# Create 1st-level analysis output graph
wf_third.write_graph(graph2use='colored', format='png', simple_form=True)

# Visualize the graph
from IPython.display import Image
Image(filename=os.path.join(wf_third.base_dir, wf_third.name, 'graph.png'))
../_images/f6d933650d8d90cff1ad80495d24bee603aaad2d482c8980c12a08ab8c126985.png
%%capture
# Run Workflow
wf_third.run(plugin="MultiProc", plugin_args={"n_procs": 8})

Unthresholded COPEs#

nipype_cope1 = './output_level3/level3_results/contrast_1/stats_dir/_contr_id_1/stats/cope1.nii.gz'
nipype_cope2 = './output_level3/level3_results/contrast_2/stats_dir/_contr_id_2/stats/cope1.nii.gz'
nipype_cope3 = './output_level3/level3_results/contrast_3/stats_dir/_contr_id_3/stats/cope1.nii.gz'
nipype_cope4 = './output_level3/level3_results/contrast_4/stats_dir/_contr_id_4/stats/cope1.nii.gz'
t_nipype_cope1 = surface.vol_to_surf(nipype_cope1, fsaverage.pial_right)
plotting.plot_surf_stat_map(
        fsaverage.infl_right, t_nipype_cope1, hemi='right',
        title="Nipype - COPE1", colorbar=False,
        cmap = 'bwr', bg_map=fsaverage.sulc_right)
plt.show()


t_nipype_cope2 = surface.vol_to_surf(nipype_cope2, fsaverage.pial_right)
plotting.plot_surf_stat_map(
        fsaverage.infl_right, t_nipype_cope2, hemi='right',
        title="Nipype - COPE2", colorbar=False,
        cmap = 'bwr', bg_map=fsaverage.sulc_right)
plt.show()


t_nipype_cope3 = surface.vol_to_surf(nipype_cope3, fsaverage.pial_right)
plotting.plot_surf_stat_map(
        fsaverage.infl_right, t_nipype_cope3, hemi='right',
        title="Nipype - COPE3", colorbar=False,
        cmap = 'bwr', bg_map=fsaverage.sulc_right)
plt.show()


t_nipype_cope4 = surface.vol_to_surf(nipype_cope4, fsaverage.pial_right)
plotting.plot_surf_stat_map(
        fsaverage.infl_right, t_nipype_cope4, hemi='right',
        title="Nipype - COPE4", colorbar=False,
        cmap = 'bwr', bg_map=fsaverage.sulc_right)
plt.show()
../_images/9d99f0e9955b0dd05bd8d559317160c66a6de377e9359720f72c5bb4f4097918.png ../_images/d7597b8a3c27323bc61b31d52aa0de87050572be410b321db85b3bd4b3bcd6c7.png ../_images/adddb5a70e683d8dd2c61d0019fe2562e1f6ac10cae42a4e3c8ad20c8b975a20.png ../_images/5651903fdea6b4f5ffdf0559b002e441a00e392a76b731d15f48b84997347ba8.png

Thresholded (p < 0.05) zstats#

nipype_cope1 = './output_level3/level3_results/contrast_1/thresholded/_contr_id_1/zstat1_maths_threshold.nii.gz'
nipype_cope2 = './output_level3/level3_results/contrast_2/thresholded/_contr_id_2/zstat1_maths_threshold.nii.gz'
nipype_cope3 = './output_level3/level3_results/contrast_3/thresholded/_contr_id_3/zstat1_maths_threshold.nii.gz'
nipype_cope4 = './output_level3/level3_results/contrast_4/thresholded/_contr_id_4/zstat1_maths_threshold.nii.gz'
plotting.plot_img_on_surf(nipype_cope1,
                          views=['lateral'],
                          hemispheres=['right'],
                          title="Nipype - Thresholded-zstats1",
                          cmap = 'hot', colorbar=True)
plotting.show()

plotting.plot_img_on_surf(nipype_cope2,
                          views=['lateral'],
                          hemispheres=['right'],
                          title="Nipype - Thresholded-zstats2",
                          cmap = 'hot', colorbar=True)
plotting.show()

plotting.plot_img_on_surf(nipype_cope3,
                          views=['lateral'],
                          hemispheres=['right'],
                          title="Nipype - Thresholded-zstats3",
                          cmap = 'hot', colorbar=True)
plotting.show()

plotting.plot_img_on_surf(nipype_cope4,
                          views=['lateral'],
                          hemispheres=['right'],
                          title="Nipype - Thresholded-zstats4",
                          cmap = 'hot', colorbar=True)
plotting.show()
../_images/fb44e1fce4b4aab248a9c23a766496b0ae3d3a104c9054a833deabfe36ca12fb.png ../_images/8bf8c54cf56c7f18bb0a200f880e178dee70af1d7c477399c23b256b93613fe2.png
/opt/conda/lib/python3.10/site-packages/nilearn/plotting/surf_plotting.py:508: RuntimeWarning: invalid value encountered in divide
  data_copy /= (vmax - vmin)
../_images/4aadca3f2e660e46196def100e021dc4db585c30a8b4e71524abd670845dfda0.png
/opt/conda/lib/python3.10/site-packages/nilearn/plotting/surf_plotting.py:508: RuntimeWarning: invalid value encountered in divide
  data_copy /= (vmax - vmin)
../_images/9d70c80cfc750d580eb4474d3170e9b10a99cba2169c59df8bfa05f407a61f98.png