A preprocessing workflow for Siemens resting state data.
This workflow makes use of:
For example:
python rsfmri_preprocessing.py -d /data/12345-34-1.dcm -f /data/Resting.nii
-s subj001 -n 2 --despike -o output
-p PBS --plugin_args "dict(qsub_args='-q many')"
This workflow takes resting timeseries and a Siemens dicom file corresponding to it and preprocesses it to produce timeseries coordinates or grayordinates.
This workflow also requires 2mm subcortical atlas and templates that are available from:
http://mindboggle.info/data.html
specifically the 2mm versions of:
The 2mm version was generated with:
>>> from nipype import freesurfer as fs
>>> rs = fs.Resample()
>>> rs.inputs.in_file = 'OASIS-TRT-20_jointfusion_DKT31_CMA_labels_in_MNI152.nii.gz'
>>> rs.inputs.resampled_file = 'OASIS-TRT-20_jointfusion_DKT31_CMA_labels_in_MNI152_2mm.nii.gz'
>>> rs.inputs.voxel_size = (2., 2., 2.)
>>> rs.inputs.args = '-rt nearest -ns 1'
>>> res = rs.run()
import os
from nipype.interfaces.base import CommandLine
CommandLine.set_default_terminal_output('file')
from nipype import config
config.enable_provenance()
from nipype import (ants, afni, fsl, freesurfer, nipy, Function, DataSink)
from nipype import Workflow, Node, MapNode
from nipype.algorithms.rapidart import ArtifactDetect
from nipype.algorithms.misc import TSNR
from nipype.interfaces.fsl import EPIDeWarp
from nipype.interfaces.io import FreeSurferSource
from nipype.interfaces.c3 import C3dAffineTool
from nipype.interfaces.utility import Merge, IdentityInterface
from nipype.utils.filemanip import filename_to_list
import numpy as np
import scipy as sp
import nibabel as nb
from dcmstack.extract import default_extractor
from dicom import read_file
imports = ['import os',
'import nibabel as nb',
'import numpy as np',
'import scipy as sp',
'from nipype.utils.filemanip import filename_to_list'
]
def get_info(dicom_files):
"""Given a Siemens dicom file return metadata
Returns
-------
RepetitionTime
Slice Acquisition Times
Spacing between slices
"""
meta = default_extractor(read_file(filename_to_list(dicom_files)[0],
stop_before_pixels=True,
force=True))
return (meta['RepetitionTime']/1000., meta['CsaImage.MosaicRefAcqTimes'],
meta['SpacingBetweenSlices'])
def median(in_files):
"""Computes an average of the median of each realigned timeseries
Parameters
----------
in_files: one or more realigned Nifti 4D time series
Returns
-------
out_file: a 3D Nifti file
"""
average = None
for idx, filename in enumerate(filename_to_list(in_files)):
img = nb.load(filename)
data = np.median(img.get_data(), axis=3)
if not average:
average = data
else:
average = average + data
median_img = nb.Nifti1Image(average/float(idx + 1),
img.get_affine(), img.get_header())
filename = os.path.join(os.getcwd(), 'median.nii.gz')
median_img.to_filename(filename)
return filename
def get_aparc_aseg(files):
"""Return the aparc+aseg.mgz file"""
for name in files:
if 'aparc+aseg.mgz' in name:
return name
raise ValueError('aparc+aseg.mgz not found')
def motion_regressors(motion_params, order=2, derivatives=2):
"""Compute motion regressors upto given order and derivative
motion + d(motion)/dt + d2(motion)/dt2 (linear + quadratic)
"""
out_files = []
for idx, filename in enumerate(filename_to_list(motion_params)):
params = np.genfromtxt(filename)
out_params = params
for d in range(1, derivatives + 1):
cparams = np.vstack((np.repeat(params[0, :][None, :], d, axis=0),
params))
out_params = np.hstack((out_params, np.diff(cparams, d, axis=0)))
out_params2 = out_params
for i in range(2, order + 1):
out_params2 = np.hstack((out_params2, np.power(out_params, i)))
filename = os.path.join(os.getcwd(), "motion_regressor%02d.txt" % idx)
np.savetxt(filename, out_params2, fmt="%.10f")
out_files.append(filename)
return out_files
def build_filter1(motion_params, comp_norm, outliers):
"""Builds a regressor set comprisong motion parameters, composite norm and
outliers
The outliers are added as a single time point column for each outlier
Parameters
----------
motion_params: a text file containing motion parameters and its derivatives
comp_norm: a text file containing the composite norm
outliers: a text file containing 0-based outlier indices
Returns
-------
components_file: a text file containing all the regressors
"""
out_files = []
for idx, filename in enumerate(filename_to_list(motion_params)):
params = np.genfromtxt(filename)
norm_val = np.genfromtxt(filename_to_list(comp_norm)[idx])
out_params = np.hstack((params, norm_val[:, None]))
try:
outlier_val = np.genfromtxt(filename_to_list(outliers)[idx])
except IOError:
outlier_val = np.empty((0))
for index in np.atleast_1d(outlier_val):
outlier_vector = np.zeros((out_params.shape[0], 1))
outlier_vector[index] = 1
out_params = np.hstack((out_params, outlier_vector))
filename = os.path.join(os.getcwd(), "filter_regressor%02d.txt" % idx)
np.savetxt(filename, out_params, fmt="%.10f")
out_files.append(filename)
return out_files
def extract_noise_components(realigned_file, mask_file, num_components=6):
"""Derive components most reflective of physiological noise
Parameters
----------
realigned_file: a 4D Nifti file containing realigned volumes
mask_file: a 3D Nifti file containing white matter + ventricular masks
num_components: number of components to use for noise decomposition
Returns
-------
components_file: a text file containing the noise components
"""
imgseries = nb.load(realigned_file)
noise_mask = nb.load(mask_file)
voxel_timecourses = imgseries.get_data()[np.nonzero(noise_mask.get_data())]
voxel_timecourses = voxel_timecourses.byteswap().newbyteorder()
voxel_timecourses[np.isnan(np.sum(voxel_timecourses, axis=1)), :] = 0
_, _, v = sp.linalg.svd(voxel_timecourses, full_matrices=False)
components_file = os.path.join(os.getcwd(), 'noise_components.txt')
np.savetxt(components_file, v[:num_components, :].T)
return components_file
def extract_subrois(timeseries_file, label_file, indices):
"""Extract voxel time courses for each subcortical roi index
Parameters
----------
timeseries_file: a 4D Nifti file
label_file: a 3D file containing rois in the same space/size of the 4D file
indices: a list of indices for ROIs to extract.
Returns
-------
out_file: a text file containing time courses for each voxel of each roi
The first four columns are: freesurfer index, i, j, k positions in the
label file
"""
img = nb.load(timeseries_file)
data = img.get_data()
roiimg = nb.load(label_file)
rois = roiimg.get_data()
out_ts_file = os.path.join(os.getcwd(), 'subcortical_timeseries.txt')
with open(out_ts_file, 'wt') as fp:
for fsindex in indices:
ijk = np.nonzero(rois == fsindex)
ts = data[ijk]
for i0, row in enumerate(ts):
fp.write('%d,%d,%d,%d,' % (fsindex, ijk[0][i0],
ijk[1][i0], ijk[2][i0]) +
','.join(['%.10f' % val for val in row]) + '\n')
return out_ts_file
def combine_hemi(left, right):
"""Combine left and right hemisphere time series into a single text file
"""
lh_data = nb.load(left).get_data()
rh_data = nb.load(right).get_data()
indices = np.vstack((1000000 + np.arange(0, lh_data.shape[0])[:, None],
2000000 + np.arange(0, rh_data.shape[0])[:, None]))
all_data = np.hstack((indices, np.vstack((lh_data.squeeze(),
rh_data.squeeze()))))
filename = 'combined_surf.txt'
np.savetxt(filename, all_data,
fmt=','.join(['%d'] + ['%.10f'] * (all_data.shape[1] - 1)))
return os.path.abspath(filename)
Creates the main preprocessing workflow
def create_workflow(files,
subject_id,
n_vol=0,
despike=True,
TR=None,
slice_times=None,
slice_thickness=None,
fieldmap_images=[],
norm_threshold=1,
num_components=6,
vol_fwhm=None,
surf_fwhm=None,
lowpass_freq=-1,
highpass_freq=-1,
sink_directory=os.getcwd(),
FM_TEdiff=2.46,
FM_sigma=2,
FM_echo_spacing=.7,
target_subject=['fsaverage3', 'fsaverage4'],
name='resting'):
wf = Workflow(name=name)
# Skip starting volumes
remove_vol = MapNode(fsl.ExtractROI(t_min=n_vol, t_size=-1),
iterfield=['in_file'],
name="remove_volumes")
remove_vol.inputs.in_file = files
# Run AFNI's despike. This is always run, however, whether this is fed to
# realign depends on the input configuration
despiker = MapNode(afni.Despike(outputtype='NIFTI_GZ'),
iterfield=['in_file'],
name='despike')
#despiker.plugin_args = {'qsub_args': '-l nodes=1:ppn='}
wf.connect(remove_vol, 'roi_file', despiker, 'in_file')
# Run Nipy joint slice timing and realignment algorithm
realign = Node(nipy.SpaceTimeRealigner(), name='realign')
realign.inputs.tr = TR
realign.inputs.slice_times = slice_times
realign.inputs.slice_info = 2
if despike:
wf.connect(despiker, 'out_file', realign, 'in_file')
else:
wf.connect(remove_vol, 'roi_file', realign, 'in_file')
# Comute TSNR on realigned data regressing polynomials upto order 2
tsnr = MapNode(TSNR(regress_poly=2), iterfield=['in_file'], name='tsnr')
wf.connect(realign, 'out_file', tsnr, 'in_file')
# Compute the median image across runs
calc_median = Node(Function(input_names=['in_files'],
output_names=['median_file'],
function=median,
imports=imports),
name='median')
wf.connect(tsnr, 'detrended_file', calc_median, 'in_files')
# Coregister the median to the surface
register = Node(freesurfer.BBRegister(),
name='bbregister')
register.inputs.subject_id = subject_id
register.inputs.init = 'fsl'
register.inputs.contrast_type = 't2'
register.inputs.out_fsl_file = True
register.inputs.epi_mask = True
# Compute fieldmaps and unwarp using them
if fieldmap_images:
fieldmap = Node(interface=EPIDeWarp(), name='fieldmap_unwarp')
fieldmap.inputs.tediff = FM_TEdiff
fieldmap.inputs.esp = FM_echo_spacing
fieldmap.inputs.sigma = FM_sigma
fieldmap.inputs.mag_file = fieldmap_images[0]
fieldmap.inputs.dph_file = fieldmap_images[1]
wf.connect(calc_median, 'median_file', fieldmap, 'exf_file')
dewarper = MapNode(interface=fsl.FUGUE(), iterfield=['in_file'],
name='dewarper')
wf.connect(tsnr, 'detrended_file', dewarper, 'in_file')
wf.connect(fieldmap, 'exf_mask', dewarper, 'mask_file')
wf.connect(fieldmap, 'vsm_file', dewarper, 'shift_in_file')
wf.connect(fieldmap, 'exfdw', register, 'source_file')
else:
wf.connect(calc_median, 'median_file', register, 'source_file')
# Get the subject's freesurfer source directory
fssource = Node(FreeSurferSource(),
name='fssource')
fssource.inputs.subject_id = subject_id
fssource.inputs.subjects_dir = os.environ['SUBJECTS_DIR']
# Extract wm+csf, brain masks by eroding freesurfer lables and then
# transform the masks into the space of the median
wmcsf = Node(freesurfer.Binarize(), name='wmcsfmask')
mask = wmcsf.clone('anatmask')
wmcsftransform = Node(freesurfer.ApplyVolTransform(inverse=True,
interp='nearest'),
name='wmcsftransform')
wmcsftransform.inputs.subjects_dir = os.environ['SUBJECTS_DIR']
wmcsf.inputs.wm_ven_csf = True
wmcsf.inputs.match = [4, 5, 14, 15, 24, 31, 43, 44, 63]
wmcsf.inputs.binary_file = 'wmcsf.nii.gz'
wmcsf.inputs.erode = int(np.ceil(slice_thickness))
wf.connect(fssource, ('aparc_aseg', get_aparc_aseg), wmcsf, 'in_file')
if fieldmap_images:
wf.connect(fieldmap, 'exf_mask', wmcsftransform, 'source_file')
else:
wf.connect(calc_median, 'median_file', wmcsftransform, 'source_file')
wf.connect(register, 'out_reg_file', wmcsftransform, 'reg_file')
wf.connect(wmcsf, 'binary_file', wmcsftransform, 'target_file')
mask.inputs.binary_file = 'mask.nii.gz'
mask.inputs.dilate = int(np.ceil(slice_thickness)) + 1
mask.inputs.erode = int(np.ceil(slice_thickness))
mask.inputs.min = 0.5
wf.connect(fssource, ('aparc_aseg', get_aparc_aseg), mask, 'in_file')
masktransform = wmcsftransform.clone("masktransform")
if fieldmap_images:
wf.connect(fieldmap, 'exf_mask', masktransform, 'source_file')
else:
wf.connect(calc_median, 'median_file', masktransform, 'source_file')
wf.connect(register, 'out_reg_file', masktransform, 'reg_file')
wf.connect(mask, 'binary_file', masktransform, 'target_file')
# Compute Art outliers
art = Node(interface=ArtifactDetect(use_differences=[True, False],
use_norm=True,
norm_threshold=norm_threshold,
zintensity_threshold=3,
parameter_source='NiPy',
bound_by_brainmask=True,
save_plot=False,
mask_type='file'),
name="art")
if fieldmap_images:
wf.connect(dewarper, 'unwarped_file', art, 'realigned_files')
else:
wf.connect(tsnr, 'detrended_file', art, 'realigned_files')
wf.connect(realign, 'par_file',
art, 'realignment_parameters')
wf.connect(masktransform, 'transformed_file', art, 'mask_file')
# Compute motion regressors
motreg = Node(Function(input_names=['motion_params', 'order',
'derivatives'],
output_names=['out_files'],
function=motion_regressors,
imports=imports),
name='getmotionregress')
wf.connect(realign, 'par_file', motreg, 'motion_params')
# Create a filter to remove motion and art confounds
createfilter1 = Node(Function(input_names=['motion_params', 'comp_norm',
'outliers'],
output_names=['out_files'],
function=build_filter1,
imports=imports),
name='makemotionbasedfilter')
wf.connect(motreg, 'out_files', createfilter1, 'motion_params')
wf.connect(art, 'norm_files', createfilter1, 'comp_norm')
wf.connect(art, 'outlier_files', createfilter1, 'outliers')
# Filter the motion and art confounds
filter1 = MapNode(fsl.GLM(out_res_name='timeseries.nii.gz',
demean=True),
iterfield=['in_file', 'design'],
name='filtermotion')
if fieldmap_images:
wf.connect(dewarper, 'unwarped_file', filter1, 'in_file')
else:
wf.connect(tsnr, 'detrended_file', filter1, 'in_file')
wf.connect(createfilter1, 'out_files', filter1, 'design')
wf.connect(masktransform, 'transformed_file', filter1, 'mask')
# Create a filter to remove noise components based on white matter and CSF
createfilter2 = MapNode(Function(input_names=['realigned_file', 'mask_file',
'num_components'],
output_names=['out_files'],
function=extract_noise_components,
imports=imports),
iterfield=['realigned_file'],
name='makecompcorrfilter')
createfilter2.inputs.num_components = num_components
wf.connect(filter1, 'out_res', createfilter2, 'realigned_file')
wf.connect(masktransform, 'transformed_file', createfilter2, 'mask_file')
# Filter noise components
filter2 = MapNode(fsl.GLM(out_res_name='timeseries_cleaned.nii.gz',
demean=True),
iterfield=['in_file', 'design'],
name='filtercompcorr')
wf.connect(filter1, 'out_res', filter2, 'in_file')
wf.connect(createfilter2, 'out_files', filter2, 'design')
wf.connect(masktransform, 'transformed_file', filter2, 'mask')
# Smoothing using surface and volume smoothing
smooth = MapNode(freesurfer.Smooth(),
iterfield=['in_file'],
name='smooth')
smooth.inputs.proj_frac_avg = (0.1, 0.9, 0.1)
if surf_fwhm is None:
surf_fwhm = 5 * slice_thickness
smooth.inputs.surface_fwhm = surf_fwhm
if vol_fwhm is None:
vol_fwhm = 2 * slice_thickness
smooth.inputs.vol_fwhm = vol_fwhm
wf.connect(filter2, 'out_res', smooth, 'in_file')
wf.connect(register, 'out_reg_file', smooth, 'reg_file')
# Bandpass filter the data
bandpass = MapNode(fsl.TemporalFilter(),
iterfield=['in_file'],
name='bandpassfilter')
if highpass_freq < 0:
bandpass.inputs.highpass_sigma = -1
else:
bandpass.inputs.highpass_sigma = 1. / (2 * TR * highpass_freq)
if lowpass_freq < 0:
bandpass.inputs.lowpass_sigma = -1
else:
bandpass.inputs.lowpass_sigma = 1. / (2 * TR * lowpass_freq)
wf.connect(smooth, 'smoothed_file', bandpass, 'in_file')
# Convert aparc to subject functional space
aparctransform = wmcsftransform.clone("aparctransform")
if fieldmap_images:
wf.connect(fieldmap, 'exf_mask', aparctransform, 'source_file')
else:
wf.connect(calc_median, 'median_file', aparctransform, 'source_file')
wf.connect(register, 'out_reg_file', aparctransform, 'reg_file')
wf.connect(fssource, ('aparc_aseg', get_aparc_aseg),
aparctransform, 'target_file')
# Sample the average time series in aparc ROIs
sampleaparc = MapNode(freesurfer.SegStats(avgwf_txt_file=True,
default_color_table=True),
iterfield=['in_file'],
name='aparc_ts')
sampleaparc.inputs.segment_id = ([8] + range(10, 14) + [17, 18, 26, 47] +
range(49, 55) + [58] + range(1001, 1036) +
range(2001, 2036))
wf.connect(aparctransform, 'transformed_file',
sampleaparc, 'segmentation_file')
wf.connect(bandpass, 'out_file', sampleaparc, 'in_file')
# Sample the time series onto the surface of the target surface. Performs
# sampling into left and right hemisphere
target = Node(IdentityInterface(fields=['target_subject']), name='target')
target.iterables = ('target_subject', filename_to_list(target_subject))
samplerlh = MapNode(freesurfer.SampleToSurface(),
iterfield=['source_file'],
name='sampler_lh')
samplerlh.inputs.sampling_method = "average"
samplerlh.inputs.sampling_range = (0.1, 0.9, 0.1)
samplerlh.inputs.sampling_units = "frac"
samplerlh.inputs.interp_method = "trilinear"
#samplerlh.inputs.cortex_mask = True
samplerlh.inputs.out_type = 'niigz'
samplerlh.inputs.subjects_dir = os.environ['SUBJECTS_DIR']
samplerrh = samplerlh.clone('sampler_rh')
samplerlh.inputs.hemi = 'lh'
wf.connect(bandpass, 'out_file', samplerlh, 'source_file')
wf.connect(register, 'out_reg_file', samplerlh, 'reg_file')
wf.connect(target, 'target_subject', samplerlh, 'target_subject')
samplerrh.set_input('hemi', 'rh')
wf.connect(bandpass, 'out_file', samplerrh, 'source_file')
wf.connect(register, 'out_reg_file', samplerrh, 'reg_file')
wf.connect(target, 'target_subject', samplerrh, 'target_subject')
# Combine left and right hemisphere to text file
combiner = MapNode(Function(input_names=['left', 'right'],
output_names=['out_file'],
function=combine_hemi,
imports=imports),
iterfield=['left', 'right'],
name="combiner")
wf.connect(samplerlh, 'out_file', combiner, 'left')
wf.connect(samplerrh, 'out_file', combiner, 'right')
# Compute registration between the subject's structural and MNI template
# This is currently set to perform a very quick registration. However, the
# registration can be made significantly more accurate for cortical
# structures by increasing the number of iterations
# All parameters are set using the example from:
# https://github.com/stnava/ANTs/blob/master/Scripts/newAntsExample.sh
reg = Node(ants.Registration(), name='antsRegister')
reg.inputs.output_transform_prefix = "output_"
reg.inputs.transforms = ['Translation', 'Rigid', 'Affine', 'SyN']
reg.inputs.transform_parameters = [(0.1,), (0.1,), (0.1,), (0.2, 3.0, 0.0)]
# reg.inputs.number_of_iterations = ([[10000, 111110, 11110]]*3 +
# [[100, 50, 30]])
reg.inputs.number_of_iterations = [[100, 100, 100]] * 3 + [[100, 20, 10]]
reg.inputs.dimension = 3
reg.inputs.write_composite_transform = True
reg.inputs.collapse_output_transforms = False
reg.inputs.metric = ['Mattes'] * 3 + [['Mattes', 'CC']]
reg.inputs.metric_weight = [1] * 3 + [[0.5, 0.5]]
reg.inputs.radius_or_number_of_bins = [32] * 3 + [[32, 4]]
reg.inputs.sampling_strategy = ['Regular'] * 3 + [[None, None]]
reg.inputs.sampling_percentage = [0.3] * 3 + [[None, None]]
reg.inputs.convergence_threshold = [1.e-8] * 3 + [-0.01]
reg.inputs.convergence_window_size = [20] * 3 + [5]
reg.inputs.smoothing_sigmas = [[4, 2, 1]] * 3 + [[1, 0.5, 0]]
reg.inputs.sigma_units = ['vox'] * 4
reg.inputs.shrink_factors = [[6, 4, 2]] + [[3, 2, 1]]*2 + [[4, 2, 1]]
reg.inputs.use_estimate_learning_rate_once = [True] * 4
reg.inputs.use_histogram_matching = [False] * 3 + [True]
reg.inputs.output_warped_image = 'output_warped_image.nii.gz'
reg.inputs.fixed_image = \
os.path.abspath('OASIS-30_Atropos_template_in_MNI152_2mm.nii.gz')
reg.inputs.num_threads = 4
reg.plugin_args = {'qsub_args': '-l nodes=1:ppn=4'}
# Convert T1.mgz to nifti for using with ANTS
convert = Node(freesurfer.MRIConvert(out_type='niigz'), name='convert2nii')
wf.connect(fssource, 'T1', convert, 'in_file')
# Mask the T1.mgz file with the brain mask computed earlier
maskT1 = Node(fsl.BinaryMaths(operation='mul'), name='maskT1')
wf.connect(mask, 'binary_file', maskT1, 'operand_file')
wf.connect(convert, 'out_file', maskT1, 'in_file')
wf.connect(maskT1, 'out_file', reg, 'moving_image')
# Convert the BBRegister transformation to ANTS ITK format
convert2itk = MapNode(C3dAffineTool(),
iterfield=['transform_file', 'source_file'],
name='convert2itk')
convert2itk.inputs.fsl2ras = True
convert2itk.inputs.itk_transform = True
wf.connect(register, 'out_fsl_file', convert2itk, 'transform_file')
if fieldmap_images:
wf.connect(fieldmap, 'exf_mask', convert2itk, 'source_file')
else:
wf.connect(calc_median, 'median_file', convert2itk, 'source_file')
wf.connect(convert, 'out_file', convert2itk, 'reference_file')
# Concatenate the affine and ants transforms into a list
pickfirst = lambda x: x[0]
merge = MapNode(Merge(2), iterfield=['in2'], name='mergexfm')
wf.connect(convert2itk, 'itk_transform', merge, 'in2')
wf.connect(reg, ('composite_transform', pickfirst), merge, 'in1')
# Apply the combined transform to the time series file
sample2mni = MapNode(ants.ApplyTransforms(),
iterfield=['input_image', 'transforms'],
name='sample2mni')
sample2mni.inputs.input_image_type = 3
sample2mni.inputs.interpolation = 'BSpline'
sample2mni.inputs.invert_transform_flags = [False, False]
sample2mni.inputs.reference_image = \
os.path.abspath('OASIS-30_Atropos_template_in_MNI152_2mm.nii.gz')
sample2mni.inputs.terminal_output = 'file'
wf.connect(bandpass, 'out_file', sample2mni, 'input_image')
wf.connect(merge, 'out', sample2mni, 'transforms')
# Sample the time series file for each subcortical roi
ts2txt = MapNode(Function(input_names=['timeseries_file', 'label_file',
'indices'],
output_names=['out_file'],
function=extract_subrois,
imports=imports),
iterfield=['timeseries_file'],
name='getsubcortts')
ts2txt.inputs.indices = [8] + range(10, 14) + [17, 18, 26, 47] +\
range(49, 55) + [58]
ts2txt.inputs.label_file = \
os.path.abspath(('OASIS-TRT-20_jointfusion_DKT31_CMA_labels_in_MNI152_'
'2mm.nii.gz'))
wf.connect(sample2mni, 'output_image', ts2txt, 'timeseries_file')
# Save the relevant data into an output directory
datasink = Node(interface=DataSink(), name="datasink")
datasink.inputs.base_directory = sink_directory
datasink.inputs.container = subject_id
datasink.inputs.substitutions = [('_target_subject_', '')]
datasink.inputs.regexp_substitutions = (r'(/_.*(\d+/))', r'/run\2')
wf.connect(despiker, 'out_file', datasink, 'resting.qa.despike')
wf.connect(realign, 'par_file', datasink, 'resting.qa.motion')
wf.connect(tsnr, 'tsnr_file', datasink, 'resting.qa.tsnr')
wf.connect(tsnr, 'mean_file', datasink, 'resting.qa.tsnr.@mean')
wf.connect(tsnr, 'stddev_file', datasink, 'resting.qa.@tsnr_stddev')
if fieldmap_images:
wf.connect(fieldmap, 'exf_mask', datasink, 'resting.reference')
else:
wf.connect(calc_median, 'median_file', datasink, 'resting.reference')
wf.connect(art, 'norm_files', datasink, 'resting.qa.art.@norm')
wf.connect(art, 'intensity_files', datasink, 'resting.qa.art.@intensity')
wf.connect(art, 'outlier_files', datasink, 'resting.qa.art.@outlier_files')
wf.connect(mask, 'binary_file', datasink, 'resting.mask')
wf.connect(masktransform, 'transformed_file',
datasink, 'resting.mask.@transformed_file')
wf.connect(register, 'out_reg_file', datasink, 'resting.registration.bbreg')
wf.connect(reg, ('composite_transform', pickfirst),
datasink, 'resting.registration.ants')
wf.connect(register, 'min_cost_file',
datasink, 'resting.qa.bbreg.@mincost')
wf.connect(smooth, 'smoothed_file', datasink, 'resting.timeseries.fullpass')
wf.connect(bandpass, 'out_file', datasink, 'resting.timeseries.bandpassed')
wf.connect(sample2mni, 'output_image', datasink, 'resting.timeseries.mni')
wf.connect(createfilter1, 'out_files',
datasink, 'resting.regress.@regressors')
wf.connect(createfilter2, 'out_files',
datasink, 'resting.regress.@compcorr')
wf.connect(sampleaparc, 'summary_file',
datasink, 'resting.parcellations.aparc')
wf.connect(sampleaparc, 'avgwf_txt_file',
datasink, 'resting.parcellations.aparc.@avgwf')
wf.connect(ts2txt, 'out_file',
datasink, 'resting.parcellations.grayo.@subcortical')
datasink2 = Node(interface=DataSink(), name="datasink2")
datasink2.inputs.base_directory = sink_directory
datasink2.inputs.container = subject_id
datasink2.inputs.substitutions = [('_target_subject_', '')]
datasink2.inputs.regexp_substitutions = (r'(/_.*(\d+/))', r'/run\2')
wf.connect(combiner, 'out_file',
datasink2, 'resting.parcellations.grayo.@surface')
return wf
Creates the full workflow including getting information from dicom files
def create_resting_workflow(args, name='resting'):
TR = args.TR
slice_times = args.slice_times
slice_thickness = None
if args.dicom_file:
TR, slice_times, slice_thickness = get_info(args.dicom_file)
slice_times = (np.array(slice_times)/1000.).tolist()
if slice_thickness is None:
from nibabel import load
img = load(args.files[0])
slice_thickness = max(img.get_header().get_zooms()[:3])
kwargs = dict(files=[os.path.abspath(filename) for
filename in args.files],
subject_id=args.subject_id,
n_vol=args.n_vol,
despike=args.despike,
TR=TR,
slice_times=slice_times,
slice_thickness=slice_thickness,
lowpass_freq=args.lowpass_freq,
highpass_freq=args.highpass_freq,
sink_directory=os.path.abspath(args.sink),
name=name)
if args.field_maps:
kwargs.update(**dict(fieldmap_images=args.field_maps,
FM_TEdiff=args.TE_diff,
FM_echo_spacing=args.echo_spacing,
FM_sigma=args.sigma))
wf = create_workflow(**kwargs)
return wf
if __name__ == "__main__":
from argparse import ArgumentParser
parser = ArgumentParser(description=__doc__)
parser.add_argument("-d", "--dicom_file", dest="dicom_file",
help="an example dicom file from the resting series")
parser.add_argument("-f", "--files", dest="files", nargs="+",
help="4d nifti files for resting state",
required=True)
parser.add_argument("-s", "--subject_id", dest="subject_id",
help="FreeSurfer subject id", required=True)
parser.add_argument("-n", "--n_vol", dest="n_vol", default=0, type=int,
help="Volumes to skip at the beginning")
parser.add_argument("--despike", dest="despike", default=False,
action="store_true", help="Use despiked data")
parser.add_argument("--TR", dest="TR", default=None,
help="TR if dicom not provided in seconds")
parser.add_argument("--slice_times", dest="slice_times", nargs="+",
type=float, help="Slice times in seconds")
parser.add_argument("-l", "--lowpass_freq", dest="lowpass_freq",
default=-1, help="Low pass frequency (Hz)")
parser.add_argument("-u", "--highpass_freq", dest="highpass_freq",
default=-1, help="High pass frequency (Hz)")
parser.add_argument("-o", "--output_dir", dest="sink",
help="Output directory base")
parser.add_argument("-w", "--work_dir", dest="work_dir",
help="Output directory base")
parser.add_argument("-p", "--plugin", dest="plugin",
default='Linear',
help="Plugin to use")
parser.add_argument("--plugin_args", dest="plugin_args",
help="Plugin arguments")
parser.add_argument("--field_maps", dest="field_maps", nargs="+",
help="field map niftis")
parser.add_argument("--fm_echospacing", dest="echo_spacing", type=float,
help="field map echo spacing")
parser.add_argument("--fm_TE_diff", dest='TE_diff', type=float,
help="field map echo time difference")
parser.add_argument("--fm_sigma", dest='sigma', type=float,
help="field map sigma value")
args = parser.parse_args()
wf = create_resting_workflow(args)
if args.work_dir:
work_dir = os.path.abspath(args.work_dir)
else:
work_dir = os.getcwd()
wf.base_dir = work_dir
if args.plugin_args:
wf.run(args.plugin, plugin_args=eval(args.plugin_args))
else:
wf.run(args.plugin)
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.