Example 1: Preprocessing Workflow

This is meant as a very simple example for a preprocessing workflow. In this workflow we will conduct the following steps:

  1. Motion Correction of functional images with FSL's MCFLIRT
  2. Image Resampling of anatomical image to a resolution of 3mmx3mmx3mm voxel size with AFNI's 3dresample
  3. Image Correction of motion corrected functional images to the resampled anatomical image with FSL's FLIRT
  4. Smoothing of coregistrated functional images with FWHM set to 4mm and 8mm

For this example, we have 3 subjects with 2 different runs each. As a short recap, the image properties are:


In [ ]:
!nib-ls /data/ds102/sub-01/*/*.nii.gz


/data/ds102/sub-01/anat/sub-01_T1w.nii.gz                     int16 [176, 256, 256]      1.00x1.00x1.00
/data/ds102/sub-01/func/sub-01_task-flanker_run-1_bold.nii.gz int16 [ 64,  64,  40, 146] 3.00x3.00x4.00x2.00
/data/ds102/sub-01/func/sub-01_task-flanker_run-2_bold.nii.gz int16 [ 64,  64,  40, 146] 3.00x3.00x4.00x2.00

So, let's start!

Imports

First, let's import all modules we later will be needing.


In [ ]:
%pylab inline
from os.path import join as opj
from nipype.interfaces.fsl import MCFLIRT, FLIRT
from nipype.interfaces.afni import Resample
from nipype.interfaces.spm import Smooth
from nipype.interfaces.utility import IdentityInterface
from nipype.interfaces.io import SelectFiles, DataSink
from nipype.pipeline.engine import Workflow, Node


Populating the interactive namespace from numpy and matplotlib

Experiment parameters

It's always a good idea to specify all parameters that might change between experiments at the beginning of your script.


In [ ]:
experiment_dir = '/output'
output_dir = 'datasink'
working_dir = 'workingdir'

# list of subject identifiers
subject_list = ['sub-01', 'sub-02', 'sub-03', 'sub-04', 'sub-05']

# list of session identifiers
session_list = ['run-1', 'run-2']

# Smoothing widths to apply
fwhm = [4, 8]

Specify Nodes

Initiate all the different interfaces (represented as nodes) that you want to use in your workflow.


In [ ]:
# MCFLIRT - motion correction
mcflirt = Node(MCFLIRT(mean_vol=True,
                       save_plots=True,
                       output_type='NIFTI'),
               name="mcflirt")

# Resample - resample anatomy to 3x3x3 voxel resolution
resample = Node(Resample(voxel_size=(3, 3, 3),
                         outputtype='NIFTI'),
                name="resample")

# FLIRT - coregister functional images to anatomical images
coreg_step1 = Node(FLIRT(output_type='NIFTI'), name="coreg_step1")
coreg_step2 = Node(FLIRT(output_type='NIFTI',
                         apply_xfm=True), name="coreg_step2")

# Smooth - image smoothing
smooth = Node(Smooth(), name="smooth")
smooth.iterables = ("fwhm", fwhm)

Specify input & output stream

Specify where the input data can be found & where and how to save the output data.


In [ ]:
# Infosource - a function free node to iterate over the list of subject names
infosource = Node(IdentityInterface(fields=['subject_id', 'session_id']),
                  name="infosource")
infosource.iterables = [('subject_id', subject_list),
                        ('session_id', session_list)]

# SelectFiles - to grab the data (alternativ to DataGrabber)
anat_file = opj('{subject_id}', 'anat', '{subject_id}_T1w.nii.gz')
func_file = opj('{subject_id}', 'func',
                '{subject_id}_task-flanker_{session_id}_bold.nii.gz')

templates = {'anat': anat_file,
             'func': func_file}
selectfiles = Node(SelectFiles(templates,
                               base_directory='/data/ds102'),
                   name="selectfiles")

# Datasink - creates output folder for important outputs
datasink = Node(DataSink(base_directory=experiment_dir,
                         container=output_dir),
                name="datasink")

# Use the following DataSink output substitutions
substitutions = [('_subject_id', ''),
                 ('_session_id_', ''),
                 ('_task-flanker', ''),
                 ('_mcf.nii_mean_reg', '_mean'),
                 ('.nii.par', '.par'),
                 ]
subjFolders = [('%s_%s/' % (sess, sub), '%s/%s' % (sub, sess))
               for sess in session_list
               for sub in subject_list]
subjFolders += [('%s_%s' % (sub, sess), '')
                for sess in session_list
                for sub in subject_list]
subjFolders += [('%s%s_' % (sess, sub), '')
                for sess in session_list
                for sub in subject_list]
substitutions.extend(subjFolders)
datasink.inputs.substitutions = substitutions

Specify Workflow

Create a workflow and connect the interface nodes and the I/O stream to each other.


In [ ]:
# Create a preprocessing workflow
preproc = Workflow(name='preproc')
preproc.base_dir = opj(experiment_dir, working_dir)

# Connect all components of the preprocessing workflow
preproc.connect([(infosource, selectfiles, [('subject_id', 'subject_id'),
                                            ('session_id', 'session_id')]),
                 (selectfiles, mcflirt, [('func', 'in_file')]),
                 (selectfiles, resample, [('anat', 'in_file')]),

                 (mcflirt, coreg_step1, [('mean_img', 'in_file')]),
                 (resample, coreg_step1, [('out_file', 'reference')]),

                 (mcflirt, coreg_step2, [('out_file', 'in_file')]),
                 (resample, coreg_step2, [('out_file', 'reference')]),
                 (coreg_step1, coreg_step2, [('out_matrix_file',
                                              'in_matrix_file')]),

                 (coreg_step2, smooth, [('out_file', 'in_files')]),

                 (mcflirt, datasink, [('par_file', 'preproc.@par')]),
                 (resample, datasink, [('out_file', 'preproc.@resample')]),
                 (coreg_step1, datasink, [('out_file', 'preproc.@coregmean')]),
                 (smooth, datasink, [('smoothed_files', 'preproc.@smooth')]),
                 ])

Visualize the workflow

It always helps to visualize your workflow.


In [ ]:
# Create preproc output graph
preproc.write_graph(graph2use='colored', format='png', simple_form=True)

# Visualize the graph
from IPython.display import Image
Image(filename=opj(preproc.base_dir, 'preproc', 'graph.dot.png'))


170301-21:33:58,238 workflow INFO:
	 Converting dotfile: /output/workingdir/preproc/graph.dot to png format
Out[ ]:

In [ ]:
# Visualize the detailed graph
preproc.write_graph(graph2use='flat', format='png', simple_form=True)
Image(filename=opj(preproc.base_dir, 'preproc', 'graph_detailed.dot.png'))


170301-21:33:58,270 workflow INFO:
	 Creating detailed dot file: /output/workingdir/preproc/graph_detailed.dot
170301-21:33:58,753 workflow INFO:
	 Creating dot file: /output/workingdir/preproc/graph.dot
Out[ ]:

Run the Workflow

Now that everything is ready, we can run the preprocessing workflow. Change n_procs to the number of jobs/cores you want to use.


In [ ]:
preproc.run('MultiProc', plugin_args={'n_procs': 4})

Inspect output

let's check the structure of the output folder, to see if we have everything we wanted to save.


In [ ]:
!tree /output/datasink/preproc/sub-01/


/output/datasink/preproc/sub-01/
├── run-1_bold_mcf.par
├── run-1_bold_mean_flirt.nii
├── run-1_fwhm_4
│   └── s_bold_mcf_flirt.nii
├── run-1_fwhm_8
│   └── s_bold_mcf_flirt.nii
├── run-2_bold_mcf.par
├── run-2_bold_mean_flirt.nii
├── run-2_fwhm_4
│   └── s_bold_mcf_flirt.nii
├── run-2_fwhm_8
│   └── s_bold_mcf_flirt.nii
└── T1w_resample.nii

4 directories, 9 files

Visualize results

let's check the effect of the different smoothing kernels.


In [ ]:
preproc.run('MultiProc', plugin_args={'n_procs': 4})
from nilearn import image, plotting
plotting.plot_epi(
    '/output/datasink/preproc/sub-01/run-1_bold_mean_flirt.nii', title="fwhm = 0mm",
    display_mode='ortho', annotate=False, draw_cross=False, cmap='gray')

mean_img = image.mean_img('/output/datasink/preproc/sub-01/run-1_fwhm_4/s_bold_mcf_flirt.nii')
plotting.plot_epi(mean_img, title="fwhm = 4mm", display_mode='ortho',
                  annotate=False, draw_cross=False, cmap='gray')

mean_img = image.mean_img('/output/datasink/preproc/sub-01/run-1_fwhm_8/s_bold_mcf_flirt.nii')
plotting.plot_epi(mean_img, title="fwhm = 8mm", display_mode='ortho',
                  annotate=False, draw_cross=False, cmap='gray')


Out[ ]:
<nilearn.plotting.displays.OrthoSlicer at 0x7fbebc997210>

How do the motion parameters look like?


In [ ]:
par = np.loadtxt('/output/datasink/preproc/sub-01/run-1_bold_mcf.par')
fig, axes = plt.subplots(2, 1, figsize=(15, 5))
axes[0].set_ylabel('rotation (radians)')
axes[0].plot(par[0:, :3])
axes[1].plot(par[0:, 3:])
axes[1].set_xlabel('time (TR)')
axes[1].set_ylabel('translation (mm)')


Out[ ]:
<matplotlib.text.Text at 0x7fbeb4e2fe10>