Home · Quickstart · Documentation · Citation · NiPy

# rsfMRI: FSL - CSF regressed

A pipeline example that uses intergrates several interfaces to perform a first and second level analysis on a two-subject data set.

1. Tell python where to find the appropriate functions.

    import numpy as np

    import nipype.interfaces.io as nio           # Data i/o
    import nipype.interfaces.fsl as fsl          # fsl
    import nipype.interfaces.utility as util     # utility
    import nipype.pipeline.engine as pe          # pypeline engine
    import nipype.algorithms.modelgen as model   # model generation
    import os                                    # system functions

    #####################################################################
    # Preliminaries

2. Setup any package specific configuration. The output file format for FSL routines is being set to uncompressed NIFTI and a specific version of matlab is being used. The uncompressed format is required because SPM does not handle compressed NIFTI.

    # Tell fsl to generate all output in compressed nifti format
    print fsl.Info.version()
    fsl.FSLCommand.set_default_output_type('NIFTI_GZ')

    extract_ref = pe.Node(interface=fsl.ExtractROI(t_min=42,
                                                   t_size=1),
                          name = 'extractref')

    # run FSL's bet
    # bet my_structural my_betted_structural

in the provided data set, the nose is behind the head and causes problems for segmentation routines

    nosestrip = pe.Node(interface=fsl.BET(frac=0.3),
                        name = 'nosestrip')
    skullstrip = pe.Node(interface=fsl.BET(mask = True),
                         name = 'stripstruct')

    refskullstrip = pe.Node(interface=fsl.BET(mask = True),
                            name = 'stripref')

    coregister = pe.Node(interface=fsl.FLIRT(dof=6),
                         name = 'coregister')

    # Preprocess functionals
    motion_correct = pe.Node(interface=fsl.MCFLIRT(save_plots = True),
                             name='realign')
                                #iterfield = ['in_file'])

skull strip functional data

    func_skullstrip = pe.Node(interface=fsl.BET(functional = True),
                              name='stripfunc')
                                 #iterfield = ['in_file'])

Run FAST on T1 anatomical image to obtain CSF mask. Create mask for three tissue types.

    getCSFmasks = pe.Node(interface=fsl.FAST(no_pve=True,segments=True),
                          name = 'segment')

Apply registration matrix to CSF segmentation mask.

    applyReg2CSFmask = pe.Node(interface=fsl.ApplyXfm(apply_xfm=True),
                               name = 'applyreg2csfmask')

Threshold CSF segmentation mask from .90 to 1

    threshCSFseg = pe.Node(interface = fsl.ImageMaths(op_string = ' -thr .90 -uthr 1 -bin '),
                           name = 'threshcsfsegmask')

Extract CSF timeseries

    avgCSF = pe.Node(interface = fsl.ImageMeants(), name='extractcsfts')

    def pickfirst(files):
        return files[0]

Create the workflow

    csffilter = pe.Workflow(name='csffilter')
    csffilter.connect([(extract_ref, motion_correct,[('roi_file', 'ref_file')]),
                       (extract_ref, refskullstrip,[('roi_file', 'in_file')]),
                       (nosestrip, skullstrip, [('out_file','in_file')]),
                       (skullstrip, getCSFmasks,[('out_file','in_files')]),
                       (skullstrip, coregister,[('mask_file','in_file')]),
                       (refskullstrip, coregister,[('out_file','reference')]),
                       (motion_correct, func_skullstrip, [('out_file', 'in_file')]),
                       (getCSFmasks, applyReg2CSFmask,[(('tissue_class_files',pickfirst),'in_file')]),
                       (refskullstrip, applyReg2CSFmask,[('out_file','reference')]),
                       (coregister, applyReg2CSFmask,[('out_matrix_file','in_matrix_file')]),
                       (applyReg2CSFmask,threshCSFseg,[('out_file','in_file')]),
                       (func_skullstrip,avgCSF,[('out_file','in_file')]),
                       (threshCSFseg,avgCSF,[('out_file','mask')]),
                       ])

    modelfit = pe.Workflow(name='modelfit')

     c. Use :class:`nipype.interfaces.spm.SpecifyModel` to generate
     SPM-specific design information.

    modelspec = pe.Node(interface=model.SpecifyModel(),  name="modelspec")

     d. Use :class:`nipype.interfaces.fsl.Level1Design` to generate a
     run specific fsf file for analysis

    level1design = pe.Node(interface=fsl.Level1Design(), name="fsfdesign")

     e. Use :class:`nipype.interfaces.fsl.FEATModel` to generate a
     run specific mat file for use by FILMGLS

    modelgen = pe.Node(interface=fsl.FEATModel(), name='modelgen')

     f. Use :class:`nipype.interfaces.fsl.FILMGLS` to estimate a model
     specified by a mat file and a functional run

    modelestimate = pe.Node(interface=fsl.FILMGLS(), name='modelestimate')
                               #iterfield = ['design_file','in_file'])

    modelfit.connect([(modelspec,level1design,[('session_info','session_info')]),
                      (level1design,modelgen,[('fsf_files','fsf_file'),
                                              ('ev_files', 'ev_files')]),
                      (modelgen,modelestimate,[('design_file','design_file')]),
                      ])

The nipype tutorial contains data for two subjects. Subject data is in two subdirectories, `s1` and `s2`. Each subject directory contains four functional volumes: f3.nii, f5.nii, f7.nii, f10.nii. And one anatomical volume named struct.nii.

Below we set some variables to inform the `datasource` about the layout of our data. We specify the location of the data, the subject sub-directories and a dictionary that maps each run to a mnemonic (or field) for the run type (`struct` or `func`). These fields become the output fields of the `datasource` node in the pipeline.

In the example below, run ‘f3’ is of type ‘func’ and gets mapped to a nifti filename through a template ‘%s.nii’. So ‘f3’ would become ‘f3.nii’.

    # Specify the location of the data.
    data_dir = os.path.abspath('data')
    # Specify the subject directories
    subject_list = ['s1']
    # Map field names to individual subject runs.
    info = dict(func=[['subject_id', ['f3',]]], #'f5','f7','f10']]],
                struct=[['subject_id','struct']])

    infosource = pe.Node(interface=util.IdentityInterface(fields=['subject_id']),
                         name="infosource")

Here we set up iteration over all the subjects. The following line is a particular example of the flexibility of the system. The `datasource` attribute `iterables` tells the pipeline engine that it should repeat the analysis on each of the items in the `subject_list`. In the current example, the entire first level preprocessing and estimation will be repeated for each subject contained in subject_list.

    infosource.iterables = ('subject_id', subject_list)

## Preprocessing pipeline nodes

Now we create a `nipype.interfaces.io.DataSource` object and fill in the information from above about the layout of our data. The `nipype.pipeline.NodeWrapper` module wraps the interface object and provides additional housekeeping and pipeline specific functionality.

    datasource = pe.Node(interface=nio.DataGrabber(infields=['subject_id'],
                                                   outfields=['func', 'struct']),
                         name = 'datasource')
    datasource.inputs.base_directory = data_dir
    datasource.inputs.template = '%s/%s.nii'
    datasource.inputs.template_args = info

     a. Setup a function that returns subject-specific information about
     the experimental paradigm. This is used by the
     :class:`nipype.modelgen.SpecifyModel` to create the
     information necessary to generate an SPM design matrix. In this
     tutorial, the same paradigm was used for every participant. Other
     examples of this function are available in the `doc/examples`
     folder. Note: Python knowledge required here.

    def subjectinfo(meantsfile):
        import numpy as np
        from nipype.interfaces.base import Bunch
        ts = np.loadtxt(meantsfile)
        output = [Bunch(regressor_names=['MeanIntensity'],
                        regressors=[ts.tolist()])]
        return output

    hpcutoff = np.inf
    TR = 3.

    modelfit.inputs.modelspec.input_units = 'secs'
    modelfit.inputs.modelspec.time_repetition = TR
    modelfit.inputs.modelspec.high_pass_filter_cutoff = hpcutoff

    modelfit.inputs.fsfdesign.interscan_interval = TR
    modelfit.inputs.fsfdesign.bases = {'none': None}
    modelfit.inputs.fsfdesign.model_serial_correlations = False

    modelfit.inputs.modelestimate.autocorr_noestimate = True

Band pass filter the data to remove frequencies below .1 Hz

    bandPassFilterData = pe.Node(interface=fsl.ImageMaths(op_string = ' -bptf 128 12.5 '),
                                 name='bandpassfiltermcdata_fslmaths')

### Set up complete workflow

    l1pipeline = pe.Workflow(name= "resting")
    l1pipeline.base_dir = os.path.abspath('./fslresting/workingdir')
    l1pipeline.connect([(infosource, datasource, [('subject_id', 'subject_id')]),
                        (datasource, csffilter, [('struct','nosestrip.in_file'),
                                                 ('func', 'realign.in_file'),
                                                 #(('func', pickfirst), 'extractref.in_file'),
                                                 ('func', 'extractref.in_file'),
                                                  ]),
                        (csffilter, modelfit, [('stripfunc.out_file', 'modelspec.functional_runs'),
                                               ('realign.par_file', 'modelspec.realignment_parameters'),
                                               (('extractcsfts.out_file', subjectinfo),'modelspec.subject_info'),
                                               ('stripfunc.out_file', 'modelestimate.in_file')
                                               ]),
                        (modelfit, bandPassFilterData, [('modelestimate.residual4d', 'in_file')]),
                        ])

    if __name__ == '__main__':
        l1pipeline.run()
        l1pipeline.write_graph()

Example source code

You can download `the full source code of this example`. This same script is also included in the Nipype source distribution under the `examples` directory.