Permalink
Switch branches/tags
Nothing to show
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
378 lines (307 sloc) 13.7 KB
#!/usr/local/bin/python
#############################################################################
# This script create and runs a preprocessing pipeline using Nipype and SPM
# See http://bit.ly/1tfaWom
#
# Copyright (C) 2014 Luke Chang
#
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#############################################################################
import matplotlib
matplotlib.use('Agg')
from nipype.interfaces import spm
import nipype.interfaces.io as nio # Data i/o
import nipype.interfaces.utility as util # utility
from nipype.pipeline.engine import Node, Workflow
from nipype.interfaces.base import BaseInterface, TraitedSpec, File, traits
import nipype.algorithms.rapidart as ra # artifact detection
from nipype.interfaces.nipy.preprocess import ComputeMask
import nipype.interfaces.matlab as mlab
import os
import nibabel as nib
from IPython.display import Image
import glob
# Specify various inputs files for pipeline
spm_path = '/Users/lukechang/Documents/MATLAB/spm8/'
data_dir='/Users/lukechang/Dropbox/PTSD/Data/Imaging'
canonical_file = spm_path + 'canonical/single_subj_T1.nii'
t1_template_file = spm_path + 'templates/T1.nii'
# Set the way matlab should be called
mlab.MatlabCommand.set_default_matlab_cmd("matlab -nodesktop -nosplash")
mlab.MatlabCommand.set_default_paths(spm_path)
###############################
## Functions
###############################
def get_n_slices(volume):
import nibabel as nib
nii = nib.load(volume)
return nii.get_shape()[2]
def get_ta(tr, n_slices):
return tr - tr/float(n_slices)
def get_slice_order(volume):
import nibabel as nib
nii = nib.load(volume)
n_slices = nii.get_shape()[2]
return range(1,n_slices+1)
def get_vox_dims(volume):
import nibabel as nib
if isinstance(volume, list):
volume = volume[0]
nii = nib.load(volume)
hdr = nii.get_header()
voxdims = hdr.get_zooms()
return [float(voxdims[0]), float(voxdims[1]), float(voxdims[2])]
class Plot_Coregistration_Montage_InputSpec(TraitedSpec):
wra_img = File(exists=True, mandatory=True)
canonical_img = File(exists=True, mandatory=True)
title = traits.Str("Normalized Functional Check", usedefault=True)
class Plot_Coregistration_Montage_OutputSpec(TraitedSpec):
plot = File(exists=True)
class Plot_Coregistration_Montage(BaseInterface):
# This function creates an axial montage of the average normalized functional data
# and overlays outline of the normalized single subject overlay.
# Could probably pick a better overlay later.
input_spec = Plot_Coregistration_Montage_InputSpec
output_spec = Plot_Coregistration_Montage_OutputSpec
def _run_interface(self, runtime):
import matplotlib
matplotlib.use('Agg')
import nibabel as nib
from nilearn import plotting, datasets, image
from nipype.interfaces.base import isdefined
import numpy as np
import pylab as plt
import os
wra_img = nib.load(self.inputs.wra_img)
canonical_img = nib.load(self.inputs.canonical_img)
title = self.inputs.title
mean_wraimg = image.mean_img(wra_img)
if title != "":
filename = title.replace(" ", "_")+".pdf"
else:
filename = "plot.pdf"
fig = plotting.plot_anat(mean_wraimg, title="wrafunc & canonical single subject", cut_coords=range(-40, 40, 10), display_mode='z')
fig.add_edges(canonical_img)
fig.savefig(filename)
fig.close()
self._plot = filename
runtime.returncode=0
return runtime
def _list_outputs(self):
outputs = self._outputs().get()
outputs["plot"] = os.path.abspath(self._plot)
return outputs
class PlotRealignmentParametersInputSpec(TraitedSpec):
realignment_parameters = File(exists=True, mandatory=True)
outlier_files = File(exists=True)
title = traits.Str("Realignment parameters", usedefault=True)
dpi = traits.Int(300, usedefault = True)
class PlotRealignmentParametersOutputSpec(TraitedSpec):
plot = File(exists=True)
class PlotRealignmentParameters(BaseInterface):
#This function is adapted from Chris Gorgolewski and creates a figure of the realignment parameters
input_spec = PlotRealignmentParametersInputSpec
output_spec = PlotRealignmentParametersOutputSpec
def _run_interface(self, runtime):
import matplotlib
matplotlib.use('Agg')
from nipype.interfaces.base import isdefined
import numpy as np
import pylab as plt
import os
realignment_parameters = np.loadtxt(self.inputs.realignment_parameters)
title = self.inputs.title
F = plt.figure(figsize=(8.3,11.7))
F.text(0.5, 0.96, self.inputs.title, horizontalalignment='center')
ax1 = plt.subplot2grid((2,2),(0,0), colspan=2)
handles =ax1.plot(realignment_parameters[:,0:3])
ax1.legend(handles, ["x translation", "y translation", "z translation"], loc=0)
ax1.set_xlabel("image #")
ax1.set_ylabel("mm")
ax1.set_xlim((0,realignment_parameters.shape[0]-1))
ax1.set_ylim(bottom = realignment_parameters[:,0:3].min(), top = realignment_parameters[:,0:3].max())
ax2 = plt.subplot2grid((2,2),(1,0), colspan=2)
handles= ax2.plot(realignment_parameters[:,3:6]*180.0/np.pi)
ax2.legend(handles, ["pitch", "roll", "yaw"], loc=0)
ax2.set_xlabel("image #")
ax2.set_ylabel("degrees")
ax2.set_xlim((0,realignment_parameters.shape[0]-1))
ax2.set_ylim(bottom=(realignment_parameters[:,3:6]*180.0/np.pi).min(), top= (realignment_parameters[:,3:6]*180.0/np.pi).max())
if isdefined(self.inputs.outlier_files):
try:
outliers = np.loadtxt(self.inputs.outlier_files)
except IOError as e:
if e.args[0] == "End-of-file reached before encountering data.":
pass
else:
raise
else:
if outliers.size > 0:
ax1.vlines(outliers, ax1.get_ylim()[0], ax1.get_ylim()[1])
ax2.vlines(outliers, ax2.get_ylim()[0], ax2.get_ylim()[1])
if title != "":
filename = title.replace(" ", "_")+".pdf"
else:
filename = "plot.pdf"
F.savefig(filename,papertype="a4",dpi=self.inputs.dpi)
plt.clf()
plt.close()
del F
self._plot = filename
runtime.returncode=0
return runtime
def _list_outputs(self):
outputs = self._outputs().get()
outputs["plot"] = os.path.abspath(self._plot)
return outputs
class Create_Covariates_InputSpec(TraitedSpec):
realignment_parameters = File(exists=True, mandatory=True)
spike_id = File(exists=True, mandatory=True)
class Create_Covariates_OutputSpec(TraitedSpec):
covariates = File(exists=True)
class Create_Covariates(BaseInterface):
input_spec = Create_Covariates_InputSpec
output_spec = Create_Covariates_OutputSpec
def _run_interface(self, runtime):
import pandas as pd
import numpy as np
ra = pd.read_table(self.inputs.realignment_parameters, header=None, sep=r"\s*", names=['ra' + str(x) for x in range(1,7)])
spike = pd.read_table(self.inputs.spike_id, header=None,names=['Spikes'])
ra = ra-ra.mean() #mean center
ra[['rasq' + str(x) for x in range(1,7)]] = ra**2 #add squared
ra[['radiff' + str(x) for x in range(1,7)]] = pd.DataFrame(ra[ra.columns[0:6]].diff()) #derivative
ra[['radiffsq' + str(x) for x in range(1,7)]] = pd.DataFrame(ra[ra.columns[0:6]].diff())**2 #derivatives squared
#build spike regressors
for i,loc in enumerate(spike['Spikes']):
ra['spike' + str(i+1)] = 0
ra['spike' + str(i+1)].iloc[loc] = 1
filename = 'covariates.csv'
ra.to_csv(filename, index=False) #write out to file
self._covariates = filename
runtime.returncode=0
return runtime
def _list_outputs(self):
outputs = self._outputs().get()
outputs["covariates"] = os.path.abspath(self._covariates)
return outputs
def create_preproc_func_pipeline(data_dir=None, subject_id=None, task_list=None):
###############################
## Set up Nodes
###############################
ds = Node(nio.DataGrabber(infields=['subject_id', 'task_id'], outfields=['func', 'struc']),name='datasource')
ds.inputs.base_directory = os.path.abspath(data_dir + '/' + subject_id)
ds.inputs.template = '*'
ds.inputs.sort_filelist = True
ds.inputs.template_args = {'func': [['task_id']], 'struc':[]}
ds.inputs.field_template = {'func': 'Functional/Raw/%s/func.nii','struc': 'Structural/SPGR/spgr.nii'}
ds.inputs.subject_id = subject_id
ds.inputs.task_id = task_list
ds.iterables = ('task_id',task_list)
# ds.run().outputs #show datafiles
# #Setup Data Sinker for writing output files
# datasink = Node(nio.DataSink(), name='sinker')
# datasink.inputs.base_directory = '/path/to/output'
# workflow.connect(realigner, 'realignment_parameters', datasink, 'motion.@par')
# datasink.inputs.substitutions = [('_variable', 'variable'),('file_subject_', '')]
#Get Timing Acquisition for slice timing
tr = 2
ta = Node(interface=util.Function(input_names=['tr', 'n_slices'], output_names=['ta'], function = get_ta), name="ta")
ta.inputs.tr=tr
#Slice Timing: sequential ascending
slice_timing = Node(interface=spm.SliceTiming(), name="slice_timing")
slice_timing.inputs.time_repetition = tr
slice_timing.inputs.ref_slice = 1
#Realignment - 6 parameters - realign to first image of very first series.
realign = Node(interface=spm.Realign(), name="realign")
realign.inputs.register_to_mean = True
#Plot Realignment
plot_realign = Node(interface=PlotRealignmentParameters(), name="plot_realign")
#Artifact Detection
art = Node(interface=ra.ArtifactDetect(), name="art")
art.inputs.use_differences = [True,False]
art.inputs.use_norm = True
art.inputs.norm_threshold = 1
art.inputs.zintensity_threshold = 3
art.inputs.mask_type = 'file'
art.inputs.parameter_source = 'SPM'
#Coregister - 12 parameters, cost function = 'nmi', fwhm 7, interpolate, don't mask
#anatomical to functional mean across all available data.
coregister = Node(interface=spm.Coregister(), name="coregister")
coregister.inputs.jobtype = 'estimate'
# Segment structural, gray/white/csf,mni,
segment = Node(interface=spm.Segment(), name="segment")
segment.inputs.save_bias_corrected = True
#Normalize - structural to MNI - then apply this to the coregistered functionals
normalize = Node(interface=spm.Normalize(), name = "normalize")
normalize.inputs.template = os.path.abspath(t1_template_file)
#Plot normalization Check
plot_normalization_check = Node(interface=Plot_Coregistration_Montage(), name="plot_normalization_check")
plot_normalization_check.inputs.canonical_img = canonical_file
#Create Mask
compute_mask = Node(interface=ComputeMask(), name="compute_mask")
#remove lower 5% of histogram of mean image
compute_mask.inputs.m = .05
#Smooth
#implicit masking (.im) = 0, dtype = 0
smooth = Node(interface=spm.Smooth(), name = "smooth")
fwhmlist = [8]
smooth.iterables = ('fwhm',fwhmlist)
#Create Covariate matrix
make_covariates = Node(interface=Create_Covariates(), name="make_covariates")
###############################
## Create Pipeline
###############################
Preprocessed = Workflow(name="Preprocessed")
Preprocessed.base_dir = os.path.abspath(data_dir + '/' + subject_id + '/Functional')
Preprocessed.connect([(ds, ta, [(('func', get_n_slices), "n_slices")]),
(ta, slice_timing, [("ta", "time_acquisition")]),
(ds, slice_timing, [('func', 'in_files'),
(('func', get_n_slices), "num_slices"),
(('func', get_slice_order), "slice_order"),
]),
(slice_timing, realign, [('timecorrected_files', 'in_files')]),
(realign, compute_mask, [('mean_image','mean_volume')]),
(realign,coregister, [('mean_image', 'target')]),
(ds,coregister, [('struc', 'source')]),
(coregister,segment, [('coregistered_source', 'data')]),
(segment, normalize, [('transformation_mat','parameter_file'),
('bias_corrected_image', 'source'),]),
(realign,normalize, [('realigned_files', 'apply_to_files'),
(('realigned_files', get_vox_dims), 'write_voxel_sizes')]),
(normalize, smooth, [('normalized_files', 'in_files')]),
(compute_mask,art,[('brain_mask','mask_file')]),
(realign,art,[('realignment_parameters','realignment_parameters')]),
(realign,art,[('realigned_files','realigned_files')]),
(realign,plot_realign, [('realignment_parameters', 'realignment_parameters')]),
(normalize, plot_normalization_check, [('normalized_files', 'wra_img')]),
(realign, make_covariates, [('realignment_parameters', 'realignment_parameters')]),
(art, make_covariates, [('outlier_files', 'spike_id')]),
])
return Preprocessed
###############################
## Create Pipeline and Run
###############################
# Create Pipeline for subject
sublist = sorted([x.split('/')[-1] for x in glob.glob(data_dir + '/subj*')])
#Run first half of subjects
#for sub in reversed(sublist[0:35]):
for sub in sublist[55:-1]:
#Glob Subject runs as they vary
runlist = [x.split('/')[-1] for x in glob.glob(data_dir + '/' + sub + '/Functional/Raw/*')]
Preprocessed = create_preproc_func_pipeline(data_dir=data_dir, subject_id = sub, task_list=runlist)
print data_dir + '/' + sub
# Write out pipeline as a DAG
Preprocessed.write_graph(dotfilename=data_dir + '/' + sub + "/Functional/Preprocessed_Workflow.dot.svg",format='svg')
# Image(filename=data_dir + '/' + sub + "/Functional/Preprocessed_Workflow.dot.png")
Preprocessed.run(plugin='MultiProc', plugin_args={'n_procs' : 12}) #Run on workstation in parallel