macjob.py

#!/usr/bin/python

import os, shutil, errno
import nipype
from nipype.interfaces import fsl, ants, c3
from nipype.caching import Memory
from nipype.interfaces.utility import IdentityInterface, Merge

#set up variables
datadir = "/scratch/users/macshine/data/nki/0102157"
bolddir = os.path.join(datadir,'session_1/RfMRI_mx_645/')
boldfile = os.path.join(bolddir,'rest.nii.gz')
t1dir = os.path.join(datadir,'anat')
t1_head = os.path.join(t1dir,'mprage.nii.gz')
t1_brain = os.path.join(t1dir,'mprage_brain2.nii.gz')
mem = Memory(base_dir='/scratch/users/chrisgor/macjob')
#results_dir = os.path.abspath("results")

#if not os.path.exists(results_dir):
#  os.mkdir(results_dir)

#def force_symlink(file1, file2):
#  os.symlink(file1, file2)


#mcflirt
mcflirt = mem.cache(fsl.MCFLIRT)
mcflirt_results = mcflirt(in_file = boldfile,
                          save_plots = True,
                          mean_vol = True)
#mcflirt_results.outputs

#save mean BOLD image
#force_symlink(mcflirt_results.outputs.mean_img, os.path.join(results_dir, "meanbold.nii.gz"))
#meanbold=os.path.join(results_dir,"meanbold.nii.gz")

#brain mask
bet = mem.cache(fsl.BET)
bet_results = bet(functional = True,
                  in_file = mcflirt_results.outputs.mean_img,
                  mask = True)
#force_symlink(bet_results.outputs.mask_file, os.path.join(results_dir, "mask.nii.gz"))
#force_symlink(bet_results.outputs.out_file, os.path.join(results_dir, "rest_brain.nii.gz"))
#force_symlink(bet_results.outputs.out_file, os.path.join(results_dir, "rest_brain_trans.nii.gz"))
#bet_results.outputs
#bet_brain=os.path.join(results_dir,"rest_brain.nii.gz")


#load motion data
#hm=numpy.loadtxt(mcflirt_results.outputs.par_file)
#force_symlink(mcflirt_results.outputs.par_file, os.path.join(results_dir, "hm.txt"))

#registration
epi_reg = mem.cache(fsl.EpiReg)
epi_reg_results = epi_reg(epi=mcflirt_results.outputs.mean_img,
                          t1_head=t1_head,
                          t1_brain=t1_brain,
                          out_base="epi2struct")
#force_symlink(epi_reg_results.outputs.epi2str_mat, os.path.join(results_dir, "epi2struct.mat"))
#epi_reg_results.outputs

#normalization using ANTs
ants_reg = mem.cache(ants.Registration)
ants_results = ants_reg(fixed_image = os.path.join(os.getenv('FSLDIR'),
                                                   'data/standard/MNI152_T1_2mm_brain.nii.gz'),
                        moving_image=t1_brain,
                        collapse_output_transforms=True,
                        initial_moving_transform_com=True,
                        num_threads=1,
                        output_warped_image=True,
                        sigma_units=['vox']*3,
                        transforms=['Rigid', 'Affine', 'SyN'],
                        winsorize_lower_quantile=0.005,
                        winsorize_upper_quantile=0.995,
                        convergence_threshold=[1e-08, 1e-08, -0.01],
                        convergence_window_size=[20, 20, 5],
                        metric=['Mattes', 'Mattes', ['Mattes', 'CC']],
                        metric_weight=[1.0, 1.0, [0.5, 0.5]],
                        number_of_iterations=[[10000, 11110, 11110], [10000, 11110, 11110], [100, 30, 20]],
                        radius_or_number_of_bins=[32, 32, [32, 4]],
                        sampling_percentage=[0.3, 0.3, [None, None]],
                        sampling_strategy=['Regular', 'Regular', [None, None]],
                        shrink_factors=[[3, 2, 1], [3, 2, 1], [4, 2, 1]],
                        smoothing_sigmas=[[4.0, 2.0, 1.0], [4.0, 2.0, 1.0], [1.0, 0.5, 0.0]],
                        transform_parameters=[(0.1,), (0.1,), (0.2, 3.0, 0.0)],
                        use_estimate_learning_rate_once=[True]*3,
                        use_histogram_matching=[False, False, True],
                        write_composite_transform=True)
#force_symlink(ants_results.outputs.composite_transform[0], os.path.join(results_dir, "t1_to_mni_warp.h5"))
#ants_results.outputs
#warp_file=ants_results.outputs.composite_transform


#Convert transforms
convert2itk=mem.cache(c3.C3dAffineTool)
convert2itk_results = convert2itk(source_file = bet_results.outputs.out_file,
                                  transform_file = epi_reg_results.outputs.epi2str_mat,
                                  reference_file = t1_brain,
                                  itk_transform = True,
                                  fsl2ras = True)
#force_symlink(c3_results.outputs.itk_transform, os.path.join(results_dir, "affine.txt"))
#c3_results.outputs

# Apply Transforms
at = mem.cache(ants.ApplyTransforms)
at_results = at(dimension=3,
                input_image = mcflirt_results.outputs.mean_img,
                reference_image = os.path.join(os.getenv('FSLDIR'),
                                               'data/standard/MNI152_T1_2mm_brain.nii.gz'),
                interpolation = 'BSpline',
                default_value = 0,
                transforms = [convert2itk_results.outputs.itk_transform, ants_results.outputs.composite_transform],
                invert_transform_flags = [False, False])

#force_symlink(at_results.outputs.output_image, os.path.join(results_dir, "rest_pp_new.nii.gz"))
#at_results.outputs