fMRI preprocessing with Nipype

Hi, I have followd Michael Notters tutorial on fmri preprocessing found here http://miykael.github.io/nipype-beginner-s-guide/firstLevel.html.

Although, I stop before he continues on to do a first-level.
I want to get the preprocessing working first. I tried with the code below.
I have just added the last few lines, trying to run it with preproc.run( ).

I seem to have a problem with nodes not running. Its something with ‘iterfield’.
Now it is specified as [‘in_file’].
Should I replace this with something?

import os
from nipype.interfaces.afni import Despike
from nipype.interfaces.freesurfer import (BBRegister, ApplyVolTransform,
Binarize, MRIConvert, FSCommand)
from nipype.interfaces.spm import (SliceTiming, Realign, Smooth, Level1Design,
EstimateModel, EstimateContrast)
from nipype.interfaces.utility import Function, IdentityInterface
from nipype.interfaces.io import FreeSurferSource, SelectFiles, DataSink
from nipype.algorithms.rapidart import ArtifactDetect
from nipype.algorithms.misc import TSNR, Gunzip
from nipype.algorithms.modelgen import SpecifySPMModel
from nipype.pipeline.engine import Workflow, Node, MapNode

MATLAB - Specify path to current SPM and the MATLAB’s default mode

from nipype.interfaces.matlab import MatlabCommand
MatlabCommand.set_default_paths(’//usr/local/MATLAB/R2012b/SPM/spm12b’)
MatlabCommand.set_default_matlab_cmd(“matlab -nodesktop -nosplash”)

FreeSurfer - Specify the location of the freesurfer folder

fs_dir = ‘/usr/local/freesurfer’
FSCommand.set_default_subjects_dir(fs_dir)

participantDir = ‘/disk2/home/granit/Link_to_TwinMR/Twins/Granits/Nipype/Output/twin_9389_20170211_102958’
os.chdir(participantDir)

#########################################

Define experiment specific parameters

#########################################
experiment_dir = ‘~/nipype_tutorial’ # location of experiment folder

subject_list = os.listdir(participantDir) # list of subject identifiers
output_dir = ‘output_fMRI_example_1st’ # name of 1st-level output folder
working_dir = ‘workingdir_fMRI_example_1st’ # name of 1st-level working directory

number_of_slices = 47 # number of slices in volume
TR = 2.396 # time repetition of volume
fwhm_size = 6 # size of FWHM in mm

#################################

Create preprocessing pipeline

#################################

Despike - Removes ‘spikes’ from the 3D+time input dataset

despike = MapNode(Despike(outputtype=‘NIFTI’),
name=“despike”, iterfield=[‘in_file’])

Slicetiming - correct for slice wise acquisition

interleaved_order = range(1,number_of_slices+1,2) + range(2,number_of_slices+1,2)
sliceTiming = Node(SliceTiming(num_slices=number_of_slices,
time_repetition=TR,
time_acquisition=TR-TR/number_of_slices,
slice_order=interleaved_order,
ref_slice=2),
name=“sliceTiming”)

Realign - correct for motion

realign = Node(Realign(register_to_mean=True),
name=“realign”)

TSNR - remove polynomials 2nd order

tsnr = MapNode(TSNR(regress_poly=2),
name=‘tsnr’, iterfield=[‘in_file’])

Artifact Detection - determine which of the images in the functional series

are outliers. This is based on deviation in intensity or movement.

art = Node(ArtifactDetect(norm_threshold=1,
zintensity_threshold=3,
mask_type=‘file’,
parameter_source=‘SPM’,
use_differences=[True, False]),
name=“art”)

Gunzip - unzip functional

gunzip = MapNode(Gunzip(), name=“gunzip”, iterfield=[‘in_file’])

Smooth - to smooth the images with a given kernel

smooth = Node(Smooth(fwhm=fwhm_size),
name=“smooth”)

FreeSurferSource - Data grabber specific for FreeSurfer data

fssource = Node(FreeSurferSource(subjects_dir=fs_dir),
run_without_submitting=True,
name=‘fssource’)

BBRegister - coregister a volume to the Freesurfer anatomical

bbregister = Node(BBRegister(init=‘header’,
contrast_type=‘t2’,
out_fsl_file=True),
name=‘bbregister’)

Volume Transformation - transform the brainmask into functional space

applyVolTrans = Node(ApplyVolTransform(inverse=True),
name=‘applyVolTrans’)

Binarize - binarize and dilate an image to create a brainmask

binarize = Node(Binarize(min=0.5,
dilate=1,
out_type=‘nii’),
name=‘binarize’)

#########################################################################

Create the preprocessing subworkflow and add all those nodes to it and#

connect them to each other

#########################################################################

Create a preprocessing workflow

preproc = Workflow(name=‘preproc’)

Connect all components of the preprocessing workflow

preproc.connect([(despike, sliceTiming, [(‘out_file’, ‘in_files’)]),
(sliceTiming, realign, [(‘timecorrected_files’, ‘in_files’)]),
(realign, tsnr, [(‘realigned_files’, ‘in_file’)]),
(tsnr, art, [(‘detrended_file’, ‘realigned_files’)]),
(realign, art, [(‘mean_image’, ‘mask_file’),
(‘realignment_parameters’,
‘realignment_parameters’)]),
(tsnr, gunzip, [(‘detrended_file’, ‘in_file’)]),
(gunzip, smooth, [(‘out_file’, ‘in_files’)]),
(realign, bbregister, [(‘mean_image’, ‘source_file’)]),
(fssource, applyVolTrans, [(‘brainmask’, ‘target_file’)]),
(bbregister, applyVolTrans, [(‘out_reg_file’, ‘reg_file’)]),
(realign, applyVolTrans, [(‘mean_image’, ‘source_file’)]),
(applyVolTrans, binarize, [(‘transformed_file’, ‘in_file’)]),
])

#########################

Run the preprocessing

#########################

preproc.base_dir = os.path.abspath(’/disk2/home/granit/Link_to_TwinMR/Twins/Granits/Nipype/Output/twin_9389_20170211_102958’)

preproc.inputs.inputspec.func = os.path.abspath(’/disk2/home/granit/Output/fp1/prox1.nii’)

preproc.inputs.inputspec.struct = os.path.abspath(’/disk2/home/granit/Output/fp1/t1.nii’)

preproc.write_graph(graph2use = ‘colored’)
preproc.run()

From what I can see, you have no input to despike.

That is my biggest concern. How do I get input in.
I have googled like crazy. After several hours, I gave up.
Tried all sorts of things.

If you have a list of input files in_list:

despike.inputs.in_file = in_list

You can do this at any point before preproc.run().

Hi,
Thanx for your reply:

I got a different output now. Among other things:

Exception: Subnodes of node: despike failed:
Subnode 0 failed
Error:
command ‘3dDespike’ could not be found on host xxx
Interface Despike failed to run.

It seems you don’t have AFNI installed, or at least not in your path. If you’re using Debian/Ubuntu, you can install it through NeuroDebian.

I have it. And added it to the part via the terminal.
Still get this message :frowning:

That’s strange. I can only guess that the path isn’t making it into your nipype instance. I’d suggest setting it in your /etc/profile (in NeuroDebian, you can add source /etc/afni/afni.sh). Then open a new terminal and verify that you can run 3dDespike.

If that works, you should be able to run your script from that terminal.

Im working on this. Not easy. Should I set export PATH=/usr/lib/afni/bin:$PATH
in /etc/profile?
And when you say in NeuroDebian. Is that a program. (My apologies for the stupid questions).
I installed NeuroDebian through the terminal way back then (6 months ago).

And last but not least, how can I run 3dDespike through the terminal. Because I am trying to run it through Python :slight_smile:

NeuroDebian is a repository for Debian/Ubuntu systems that maintains up-to-date versions of neuroscience tools, because the normal release schedules of Debian/Ubuntu often don’t match the needs of developers and users.

If you’ve set up NeuroDebian, then installing AFNI via sudo apt-get install afni will create a script in /etc/afni/afni.sh. You can edit /etc/profile to include the following line:

source /etc/afni/afni.sh

Then each time you load a terminal, that script will be run, which will update the PATH for you.

If you didn’t install AFNI through NeuroDebian, you can still update your PATH in /etc/profile:

PATH="$PATH:/usr/lib/afni/bin"

You can check your PATH in a given terminal with echo $PATH. If you see /usr/lib/afni/bin somewhere in there (or wherever you’ve installed AFNI), you should be good to go. You can again test this by running the command 3dDespike (literally type “3dDespike” in your terminal and press Enter), and you should see output like the following:

Usage: 3dDespike [options] dataset
Removes 'spikes' from the 3D+time input dataset and writes
a new dataset with the spike values replaced by something
more pleasing to the eye.

Method:
 * L1 fit a smooth-ish curve to each voxel time series
    [see -corder option for description of the curve]
    [see -NEW option for a different & faster fitting method]
 * Compute the MAD of the difference between the curve and
    the data time series (the residuals).
 * Estimate the standard deviation 'sigma' of the residuals
    as sqrt(PI/2)*MAD.
 * For each voxel value, define s = (value-curve)/sigma.
 * Values with s > c1 are replaced with a value that yields
    a modified s' = c1+(c2-c1)*tanh((s-c1)/(c2-c1)).
 * c1 is the threshold value of s for a 'spike' [default c1=2.5].
 * c2 is the upper range of the allowed deviation from the curve:
    s=[c1..infinity) is mapped to s'=[c1..c2)   [default c2=4].

Options:
 -ignore I  = Ignore the first I points in the time series:
               these values will just be copied to the
               output dataset [default I=0].
 -corder L  = Set the curve fit order to L:
               the curve that is fit to voxel data v(t) is

                       k=L [        (2*PI*k*t)          (2*PI*k*t) ]
 f(t) = a+b*t+c*t*t + SUM  [ d * sin(--------) + e * cos(--------) ]
                       k=1 [  k     (    T   )    k     (    T   ) ]

               where T = duration of time series;
               the a,b,c,d,e parameters are chosen to minimize
               the sum over t of |v(t)-f(t)| (L1 regression);
               this type of fitting is is insensitive to large
               spikes in the data.  The default value of L is
               NT/30, where NT = number of time points.

 -cut c1 c2 = Alter default values for the spike cut values
               [default c1=2.5, c2=4.0].
 -prefix pp = Save de-spiked dataset with prefix 'pp'
               [default pp='despike']
 -ssave ttt = Save 'spikiness' measure s for each voxel into a
               3D+time dataset with prefix 'ttt' [default=no save]
 -nomask    = Process all voxels
               [default=use a mask of high-intensity voxels, ]
               [as created via '3dAutomask -dilate 4 dataset'].
 -dilate nd = Dilate 'nd' times (as in 3dAutomask).  The default
               value of 'nd' is 4.
 -q[uiet]   = Don't print '++' informational messages.

 -localedit = Change the editing process to the following:
                If a voxel |s| value is >= c2, then replace
                the voxel value with the average of the two
                nearest non-spike (|s| < c2) values; the first
                one previous and the first one after.
                Note that the c1 cut value is not used here.

 -NEW       = Use the 'new' method for computing the fit, which
              should be faster than the L1 method for long time
              series (200+ time points); however, the results
              are similar but NOT identical. [29 Nov 2013]
              * You can also make the program use the 'new'
                method by setting the environment variable
                  AFNI_3dDespike_NEW
                to the value YES; as in
                  setenv AFNI_3dDespike_NEW YES  (csh)
                  export AFNI_3dDespike_NEW=YES  (bash)
              * If this variable is set to YES, you can turn off
                the '-NEW' processing by using the '-OLD' option.
          -->>* For time series more than 500 points long, the
                '-OLD' algorithm is tremendously slow.  You should
                use the '-NEW' algorith in such cases.
             ** At some indeterminate point in the future, the '-NEW'
                method will become the default!

Caveats:
* Despiking may interfere with image registration, since head
   movement may produce 'spikes' at the edge of the brain, and
   this information would be used in the registration process.
   This possibility has not been explored or calibrated.
* [LATER] Actually, it seems like the registration problem
   does NOT happen, and in fact, despiking seems to help!
* Check your data visually before and after despiking and
   registration!
   [Hint: open 2 AFNI controllers, and turn Time Lock on.]

 =========================================================================
* This binary version of 3dDespike is compiled using OpenMP, a semi-
   automatic parallelizer software toolkit, which splits the work across
   multiple CPUs/cores on the same shared memory computer.
* OpenMP is NOT like MPI -- it does not work with CPUs connected only
   by a network (e.g., OpenMP doesn't work with 'cluster' setups).
* For implementation and compilation details, please see
   https://afni.nimh.nih.gov/pub/dist/doc/misc/OpenMP.html
* The number of CPU threads used will default to the maximum number on
   your system.  You can control this value by setting environment variable
   OMP_NUM_THREADS to some smaller value (including 1).
* Un-setting OMP_NUM_THREADS resets OpenMP back to its default state of
   using all CPUs available.
   ++ However, on some systems (such as the NIH Biowulf), it seems to be
      necessary to set OMP_NUM_THREADS explicitly, or you only get one CPU.
   ++ On other systems with many CPUS, you probably want to limit the CPU
      count, since using more than (say) 16 threads is probably useless.
* You must set OMP_NUM_THREADS in the shell BEFORE running the program,
   since OpenMP queries this variable BEFORE the program actually starts.
   ++ You can't usefully set this variable in your ~/.afnirc file or on the
      command line with the '-D' option.
* How many threads are useful?  That varies with the program, and how well
   it was coded.  You'll have to experiment on your own systems!
* The number of CPUs on this particular computer system is ...... 6.
* The maximum number of CPUs that will be used is now set to .... 6.
 =========================================================================

++ Compile date = Aug 28 2016

All this does is prove that the program exists in your PATH, and so should be accessible to nipype. Then, when you run your nipype workflow, if you still get the same error, we’ll know there’s some other kind of issue.

If that happens, along with your output, if you could copy the output of the bash command env (run from the same terminal as the failure), that would be useful.

1 Like

So its definitely in the path. I tested it in the Terminal. Then I try to run the workflow from within Python and got the following message.

/disk2/home/granit/anaconda2/bin/python /disk2/home/granit/Scripts/spmPreprocessing.py
170524-03:28:08,493 workflow INFO:
 Converting dotfile: /storage/Output/Preprocess/graph.dot to png format
170524-03:28:08,723 workflow INFO:
 ['check', 'execution', 'logging']
170524-03:28:08,731 workflow INFO:
 Running serially.
170524-03:28:08,732 workflow INFO:
 Executing node despike in dir: /tmp/tmpKo_shM/preproc/despike
170524-03:28:08,739 workflow INFO:
 Executing node _despike0 in dir: /tmp/tmpKo_shM/preproc/despike/mapflow/_despike0
170524-03:28:08,743 workflow INFO:
 Running: 3dDespike -prefix /disk2/home/granit/Output/Preprocess /tmp/tmpKo_shM/preproc/despike/mapflow/_despike0/prox1.nii
 170524-03:28:08,745 workflow ERROR:
 ['Node despike failed to run on host psy-slab4.']
170524-03:28:08,745 workflow INFO:
 Saving crash info to /storage/Scripts/crash-20170524-032808-granit-despike-11bc471d-72ad-44a1-b265-7ddacce82097.pklz
170524-03:28:08,745 workflow INFO:
 Traceback (most recent call last):
File "/disk2/home/granit/anaconda2/lib/python2.7/site-packages/nipype-0.12.1-py2.7.egg/nipype/pipeline/plugins/linear.py", line 39, in run
node.run(updatehash=updatehash)
File "/disk2/home/granit/anaconda2/lib/python2.7/site-packages/nipype-0.12.1-py2.7.egg/nipype/pipeline/engine/nodes.py", line 394, in run
self._run_interface()
File "/disk2/home/granit/anaconda2/lib/python2.7/site-packages/nipype-0.12.1-py2.7.egg/nipype/pipeline/engine/nodes.py", line 1306, in _run_interface
updatehash=updatehash))
File "/disk2/home/granit/anaconda2/lib/python2.7/site-packages/nipype-0.12.1-py2.7.egg/nipype/pipeline/engine/nodes.py", line 1211, in _collate_results
(self.name, '\n'.join(msg)))
Exception: Subnodes of node: despike failed:
Subnode 0 failed
Error:
 command '3dDespike' could not be found on host psy-slab4.psyk.uu.se
 Interface Despike failed to run. 

 170524-03:28:08,753 workflow INFO:
 Executing node fssource in dir: /tmp/tmpMKFzuK/preproc/fssource
 170524-03:28:08,766 workflow INFO:
 Runtime memory and threads stats unavailable
 170524-03:28:08,766 workflow INFO:
 ***********************************
 170524-03:28:08,767 workflow ERROR:
 could not run node: preproc.despike
 170524-03:28:08,767 workflow INFO:
 crashfile: /storage/Scripts/crash-20170524-032808-granit-despike-11bc471d-72ad-44a1-b265-7ddacce82097.pklz
 170524-03:28:08,767 workflow INFO:
 ***********************************
Traceback (most recent call last):

File"/disk2/home/granit/Scripts/spmPreprocessing.py", line 152, in
preproc.run()
File “/disk2/home/granit/anaconda2/lib/python2.7/site-packages/nipype-0.12.1-py2.7.egg/nipype/pipeline/engine/workflows.py”, line 597, in run
runner.run(execgraph, updatehash=updatehash, config=self.config)
File “/disk2/home/granit/anaconda2/lib/python2.7/site-packages/nipype-0.12.1-py2.7.egg/nipype/pipeline/plugins/linear.py”, line 57, in run
report_nodes_not_run(notrun)
File “/disk2/home/granit/anaconda2/lib/python2.7/site-packages/nipype-0.12.1-py2.7.egg/nipype/pipeline/plugins/base.py”, line 95, in report_nodes_not_run
raise RuntimeError(('Workflow did not execute cleanly. ’
RuntimeError: Workflow did not execute cleanly. Check log for details

 Process finished with exit code 1

In the Terminal, I ran the following:

3dDespike -ignore 0 -nomask -prefix despiked /..../functional.nii

It produced output with .brik and something else. I then used:

3dAFNItoNIFTI ..../despike+orig.BRIK

and it created a despiked.nii file.

Then I tried the following in Python and it failed:

from nipype.interfaces import afni as afni
despike = afni.Despike()
despike.inputs.in_file = '/../functional.nii'
despike.cmdline
'3dDespike -prefix functional_despike functional.nii
res = despike.run( )

Result:

IOError: command '3dDespike' could not be found on host py-g.psyk.se
Interface Despike failed to run.

Alright, it seems that the problem with not finding 3dDespike is solved now.
I went out and sat in the sun for an hour then came back.
Instead of having the path in /etc/profile (maybe I did it wrong) I added the path in .bashrc.
At the end of the document I added: export PATH=$PATH:/usr/lib/afni/bin:$PATH

Then, my code looks different from the one that I posted before. I am not using a Workflow.
Instead it looks like this, specificaly, for despiking:

import os
from nipype.interfaces.afni import Despike

 despike = Despike( )
 despike.inputs.in_file = 'functional.nii'
 despike.inputs.out_file = 'despiked' + 'functional.nii'  # name of output file
 despike.inputs.ignore_exception = True 
 despike.inputs.outputtype = 'NIFTI'
 res = despike.run( ) 

I will try with the preprocessing steps, building on the output from despike.run( )