1#!/usr/local/bin/python3.8 2__version__="v2.5 beta11" 3welcome_block=""" 4# Multi-Echo ICA, Version %s 5# 6# Kundu, P., Brenowitz, N.D., Voon, V., Worbe, Y., Vertes, P.E., Inati, S.J., Saad, Z.S., 7# Bandettini, P.A. & Bullmore, E.T. Integrated strategy for improving functional 8# connectivity mapping using multiecho fMRI. PNAS (2013). 9# 10# Kundu, P., Inati, S.J., Evans, J.W., Luh, W.M. & Bandettini, P.A. Differentiating 11# BOLD and non-BOLD signals in fMRI time series using multi-echo EPI. NeuroImage (2011). 12# http://dx.doi.org/10.1016/j.neuroimage.2011.12.028 13# 14# meica.py version 2.5 (c) 2014 Prantik Kundu 15# PROCEDURE 1 : Preprocess multi-echo datasets and apply multi-echo ICA based on spatial concatenation 16# -Check arguments, input filenames, and filesystem for dependencies 17# -Calculation of motion parameters based on images with highest constrast 18# -Application of motion correction and coregistration parameters 19# -Misc. EPI preprocessing (temporal alignment, smoothing, etc) in appropriate order 20# -Compute PCA and ICA in conjuction with TE-dependence analysis 21""" % (__version__) 22 23import sys 24import commands 25from re import split as resplit 26import re 27from os import system,getcwd,mkdir,chdir,popen 28import os.path 29from string import rstrip,split 30from optparse import OptionParser,OptionGroup,SUPPRESS_HELP 31 32#Filename parser for NIFTI and AFNI files 33def dsprefix(idn): 34 def prefix(datasetname): 35 return split(datasetname,'+')[0] 36 if len(split(idn,'.'))!=0: 37 if split(idn,'.')[-1]=='HEAD' or split(idn,'.')[-1]=='BRIK' or split(idn,'.')[-2:]==['BRIK','gz']: 38 return prefix(idn) 39 elif split(idn,'.')[-1]=='nii' and not split(idn,'.')[-1]=='nii.gz': 40 return '.'.join(split(idn,'.')[:-1]) 41 elif split(idn,'.')[-2:]==['nii','gz']: 42 return '.'.join(split(idn,'.')[:-2]) 43 else: 44 return prefix(idn) 45 else: 46 return prefix(idn) 47 48def dssuffix(idna): 49 suffix = idna.split(dsprefix(idna))[-1] 50 #print suffix 51 spl_suffix=suffix.split('.') 52 #print spl_suffix 53 if len(spl_suffix[0])!=0 and spl_suffix[0][0] == '+': return spl_suffix[0] 54 else: return suffix 55 56def version_checker(cur_ver,ref_ver): 57 """ 58 Checks version in major/minor format of a current version versus a reference version. 59 Supports 2 or more level versioning, only 2 levels checked) 60 61 Input: 62 cur_ver: float or string of version to check 63 ref_ver: float or string of reference version 64 65 Returns: 66 bool for pass or fail 67 """ 68 cur_ver = str(cur_ver) 69 ref_ver = str(ref_ver) 70 cur_Varr = [int(vvv) for vvv in cur_ver.split('.')[0:2]] 71 cur_V = cur_Varr[0]*10**2 + cur_Varr[1]%100 72 ref_Varr = [int(vvv) for vvv in ref_ver.split('.')[0:2]] 73 ref_V = ref_Varr[0]*10**2 + ref_Varr[1]%100 74 if cur_V > ref_V: return True 75 return False 76 77#Run dependency check 78def dep_check(): 79 print '++ Checking system for dependencies...' 80 fails=0 81 numpy_installed = 0 82 scipy_installed = 0 83 python_version_ok = 0 84 global grayweight_ok 85 grayweight_ok = 0 86 print " + Python version: %s" % ('.'.join([str(v) for v in sys.version_info[0:3]])) 87 if sys.version_info < (2,6) or sys.version_info > (3,0): 88 print "*+ Python 2.x is <2.6, please upgrade to Python 2.x >= 2.6 & Numpy >= 1.5.x." 89 fails+=1 90 else: 91 python_version_ok = 1 92 try: 93 import numpy 94 numpy_installed = 1 95 except: 96 print "*+ Can't import Numpy! Please check Numpy installation for this Python version." 97 fails+=1 98 try: 99 import scipy 100 scipy_installed = 1 101 except: 102 print "*+ Can't import Scipy! Please check Scipy installation for this Python version." 103 fails+=1 104 if numpy_installed: 105 print " + Numpy version: %s" % (numpy.__version__) 106 if version_checker(numpy.__version__,1.5)==False: 107 fails+=1 108 print "*+ Numpy version is too old! Please upgrade to Numpy >=1.5.x!" 109 import numpy.__config__ as nc 110 if nc.blas_opt_info == {}: 111 fails+=1 112 print "*+ Numpy is not linked to BLAS! Please check Numpy installation." 113 if scipy_installed: 114 print " + Scipy version: %s" % (scipy.__version__) 115 if version_checker(scipy.__version__,0.11)==False: 116 fails+=1 117 print "*+ Scipy version is too old! Please upgrade to Scipy >=0.11.x!" 118 afnicheck = commands.getstatusoutput("3dinfo") 119 afnisegcheck = commands.getstatusoutput("3dSeg -help") 120 if afnicheck[0]!=0: 121 print "*+ Can't run AFNI binaries. Make sure AFNI is on the path!" 122 fails+=1 123 elif not afnicheck[1].__contains__('Alternate Alternative Usage') and len(afnisegcheck[1]) < 1000: 124 print "*+ This seems like an old version of AFNI. Please upgrade to latest version of AFNI." 125 fails+=1 126 if afnisegcheck[0]==0 and len(afnisegcheck[1]) >= 1000: 127 print " + Using AFNI 3dSeg for gray matter weighted anatomical-functional coregistration" 128 grayweight_ok = 1 129 if grayweight_ok==0: 130 print "*+ WARNING: AFNI 3dSeg not available for gray matter weighted coregistration. See README." 131 if fails==0: 132 print " + Dependencies OK." 133 else: 134 print "*+ EXITING. Please see error messages." 135 sys.exit() 136 137def getdsname(e_ii,prefixonly=False): 138 if shorthand_dsin: dsname = '%s%s%s%s' % (prefix,datasets[e_ii],trailing,isf) 139 else: dsname = datasets_in[e_ii] 140 if prefixonly: return dsprefix(dsname) 141 else: return dsname 142 143def logcomment(comment,level=3): 144 majmark='\n' 145 leading='--------' 146 if level==3: majmark='' 147 if level==1: 148 leading="+* " 149 sl.append("""echo "\n++++++++++++++++++++++++" """) 150 majmark='' 151 sl.append("""%secho %s"%s" """ % (majmark,leading,comment)) 152 153#Configure options and help dialog 154parser=OptionParser() 155parser.add_option('-e',"",dest='tes',help="ex: -e 14.5,38.5,62.5 Echo times (in ms)",default='') 156parser.add_option('-d',"",dest='dsinputs',help="ex: -d RESTe1.nii.gz,RESTe2.nii.gz,RESTe3.nii.gz",default='') 157parser.add_option('-a',"",dest='anat',help="ex: -a mprage.nii.gz Anatomical dataset (optional)",default='') 158parser.add_option('-b',"",dest='basetime',help="ex: -b 10s OR -b 10v Time to steady-state equilibration in seconds(s) or volumes(v). Default 0. ",default='0') 159parser.add_option('',"--MNI",dest='mni',action='store_true',help="Warp to MNI space using high-resolution template",default=False) 160extopts=OptionGroup(parser,"Additional processing options") 161extopts.add_option('',"--qwarp",dest='qwarp',action='store_true',help="Nonlinear anatomical normalization to MNI (or --space template) using 3dQWarp, after affine",default=False) 162extopts.add_option('',"--fres",dest='fres',help="Specify functional voxel dim. in mm (iso.) for resampling during preprocessing. Default none. ex: --fres=2.5", default=False) 163extopts.add_option('',"--space",dest='space',help="Path to specific standard space template for affine anatomical normalization",default=False) 164extopts.add_option('',"--no_skullstrip",action="store_true",dest='no_skullstrip',help="Anatomical is already intensity-normalized and skull-stripped (if -a provided)",default=False) 165extopts.add_option('',"--no_despike",action="store_true",dest='no_despike',help="Do not de-spike functional data. Default is to de-spike, recommended.",default=False) 166extopts.add_option('',"--no_axialize",action="store_true",dest='no_axialize',help="Do not re-write dataset in axial-first order. Default is to axialize, recommended.",default=False) 167extopts.add_option('',"--mask_mode",dest='mask_mode',help="Mask functional with help from anatomical or standard space images: use 'anat' or 'template' ",default='func') 168extopts.add_option('',"--coreg_mode",dest='coreg_mode',help="Coregistration with Local Pearson and T2* weights (default), or use align_epi_anat.py (edge method): use 'lp-t2s' or 'aea'",default='lp-t2s') 169extopts.add_option('',"--smooth",dest='FWHM',help="Data FWHM smoothing (3dBlurInMask). Default off. ex: --smooth 3mm ",default='0mm') 170extopts.add_option('',"--align_base",dest='align_base',help="Explicitly specify base dataset for volume registration",default='') 171extopts.add_option('',"--TR",dest='TR',help="The TR. Default read from input dataset header",default='') 172extopts.add_option('',"--tpattern",dest='tpattern',help="Slice timing (i.e. alt+z, see 3dTshift -help). Default from header. (N.B. This is important!)",default='') 173extopts.add_option('',"--align_args",dest='align_args',help="Additional arguments to anatomical-functional co-registration routine",default='') 174extopts.add_option('',"--ted_args",dest='ted_args',help="Additional arguments to TE-dependence analysis routine",default='') 175extopts.add_option('',"--daw",dest='daw',help=SUPPRESS_HELP,default='10') 176extopts.add_option('',"--tlrc",dest='space',help=SUPPRESS_HELP,default=False) #For backwards compat. with existing scripts 177extopts.add_option('',"--highpass",dest='highpass',help=SUPPRESS_HELP,default=0.0) 178extopts.add_option('',"--detrend",dest='detrend',help=SUPPRESS_HELP,default=0.) 179extopts.add_option('',"--initcost",dest='initcost',help=SUPPRESS_HELP,default='tanh') 180extopts.add_option('',"--finalcost",dest='finalcost',help=SUPPRESS_HELP,default='tanh') 181extopts.add_option('',"--sourceTEs",dest='sourceTEs',help=SUPPRESS_HELP,default='-1') 182parser.add_option_group(extopts) 183runopts=OptionGroup(parser,"Run optipns") 184runopts.add_option('',"--prefix",dest='prefix',help="Prefix for final ME-ICA output datasets.",default='') 185runopts.add_option('',"--cpus",dest='cpus',help="Maximum number of CPUs (OpenMP threads) to use. Default 2.",default='2') 186runopts.add_option('',"--label",dest='label',help="Label to tag ME-ICA analysis folder.",default='') 187runopts.add_option('',"--test_proc",action="store_true",dest='test_proc',help="Align and preprocess 1 dataset then exit, for testing",default=False) 188runopts.add_option('',"--script_only",action="store_true",dest='script_only',help="Generate script only, then exit",default=0) 189runopts.add_option('',"--pp_only",action="store_true",dest='pp_only',help="Preprocess only, then exit.",default=False) 190runopts.add_option('',"--keep_int",action="store_true",dest='keep_int',help="Keep preprocessing intermediates. Default delete.",default=False) 191runopts.add_option('',"--skip_check",action="store_true",dest='skip_check',help="Skip dependency checks during initialization.",default=False) 192runopts.add_option('',"--OVERWRITE",dest='overwrite',action="store_true",help="If meica.xyz directory exists, overwrite. ",default=False) 193parser.add_option_group(runopts) 194(options,args) = parser.parse_args() 195 196#Welcome line 197print """\n-- Multi-Echo Independent Components Analysis (ME-ICA) %s -- 198 199Please cite: 200Kundu, P., Brenowitz, N.D., Voon, V., Worbe, Y., Vertes, P.E., Inati, S.J., Saad, Z.S., 201Bandettini, P.A. & Bullmore, E.T. Integrated strategy for improving functional 202connectivity mapping using multiecho fMRI. PNAS (2013). 203 204Kundu, P., Inati, S.J., Evans, J.W., Luh, W.M. & Bandettini, P.A. Differentiating 205BOLD and non-BOLD signals in fMRI time series using multi-echo EPI. NeuroImage (2011). 206""" % (__version__) 207 208#Parse dataset input names 209if options.dsinputs=='' or options.TR==0: 210 dep_check() 211 print "*+ Need at least dataset inputs and TE. Try meica.py -h" 212 sys.exit() 213if os.path.abspath(os.path.curdir).__contains__('meica.'): 214 print "*+ You are inside a ME-ICA directory! Please leave this directory and rerun." 215 sys.exit() 216 217#Parse shorthand input file specification and TEs 218tes=split(options.tes,',') 219outprefix=options.prefix 220if '[' in options.dsinputs: 221 shorthand_dsin = True 222 dsinputs=dsprefix(options.dsinputs) 223 prefix=resplit(r'[\[\],]',dsinputs)[0] 224 datasets=resplit(r'[\[\],]',dsinputs)[1:-1] 225 trailing=resplit(r'[\]+]',dsinputs)[-1] 226 isf= dssuffix(options.dsinputs) 227 setname=prefix+''.join(datasets)+trailing+options.label 228else: 229 #Parse longhand input file specificiation 230 shorthand_dsin = False 231 datasets_in = options.dsinputs.split(',') 232 datasets = [str(vv+1) for vv in range(len(tes))] 233 prefix = dsprefix(datasets_in[0]) 234 isf = dssuffix(datasets_in[0]) 235 if '.nii' in isf: isf='.nii' 236 trailing='' 237 setname=prefix+options.label 238 239if not shorthand_dsin and len(datasets)!=len(datasets_in): 240 print "*+ Can't understand dataset specification. Try double quotes around -d argument." 241 sys.exit() 242 243if len(options.tes.split(','))!=len(datasets): 244 print "*+ Number of TEs and input datasets must be equal and matched in order. Or try double quotes around -d argument." 245 sys.exit() 246 247#Prepare script 248startdir=rstrip(popen('pwd').readlines()[0]) 249meicadir=os.path.dirname(os.path.abspath(os.path.expanduser(sys.argv[0]))) 250sl = [] #Script command list 251sl.append('#'+" ".join(sys.argv).replace('"',r"\"")) 252sl.append(welcome_block) 253osf='.nii.gz' #Using NIFTI outputs 254 255#Check if input files exist 256notfound=0 257for ds_ii in range(len(datasets)): 258 if commands.getstatusoutput('3dinfo %s' % (getdsname(ds_ii)))[0]!=0: 259 print "*+ Can't find/load dataset %s !" % (getdsname(ds_ii)) 260 notfound+=1 261if options.anat!='' and commands.getstatusoutput('3dinfo %s' % (options.anat))[0]!=0: 262 print "*+ Can't find/load anatomical dataset %s !" % (options.anat) 263 notfound+=1 264if notfound!=0: 265 print "++ EXITING. Check dataset names." 266 sys.exit() 267 268#Check dependencies 269grayweight_ok = 0 270if not options.skip_check: 271 dep_check() 272 print "++ Continuing with preprocessing." 273else: 274 print "*+ Skipping dependency checks." 275 grayweight_ok = 1 276 277#Parse timing arguments 278if options.TR!='':tr=float(options.TR) 279else: 280 tr=float(os.popen('3dinfo -tr %s' % (getdsname(0))).readlines()[0].strip()) 281 options.TR=str(tr) 282if 'v' in str(options.basetime): 283 basebrik = int(options.basetime.strip('v')) 284else: 285 timetoclip=0 286 timetoclip = float(options.basetime.strip('s')) 287 basebrik=int(round(timetoclip/tr)) 288 289#Misc. command parsing 290if options.mni: options.space='MNI_caez_N27+tlrc' 291if options.qwarp and (options.anat=='' or not options.space): 292 print "*+ Can't specify Qwarp nonlinear coregistration without anatomical and SPACE template!" 293 sys.exit() 294 295if not options.mask_mode in ['func','anat','template']: 296 print "*+ Mask mode option '%s' is not recognized!" % options.mask_mode 297 sys.exit() 298 299#Parse alignment options 300if options.coreg_mode == 'aea': options.t2salign=False 301elif 'lp' in options.coreg_mode : options.t2salign=True 302align_base = basebrik 303align_interp='cubic' 304align_interp_final='wsinc5' 305oblique_epi_read = 0 306oblique_anat_read = 0 307zeropad_opts = " -I %s -S %s -A %s -P %s -L %s -R %s " % (tuple([1]*6)) 308if options.anat!='': 309 oblique_anat_read = int(os.popen('3dinfo -is_oblique %s' % (options.anat)).readlines()[0].strip()) 310 epicm = [float(coord) for coord in os.popen("3dCM %s" % (getdsname(0))).readlines()[0].strip().split()] 311 anatcm = [float(coord) for coord in os.popen("3dCM %s" % (options.anat)).readlines()[0].strip().split()] 312 maxvoxsz = float(os.popen("3dinfo -dk %s" % (getdsname(0))).readlines()[0].strip()) 313 deltas = [abs(epicm[0]-anatcm[0]),abs(epicm[1]-anatcm[1]),abs(epicm[2]-anatcm[2])] 314 cmdist = 20+sum([dd**2. for dd in deltas])**.5 315 cmdif = max(abs(epicm[0]-anatcm[0]),abs(epicm[1]-anatcm[1]),abs(epicm[2]-anatcm[2])) 316 addslabs = abs(int(cmdif/maxvoxsz))+10 317 zeropad_opts=" -I %s -S %s -A %s -P %s -L %s -R %s " % (tuple([addslabs]*6)) 318oblique_epi_read = int(os.popen('3dinfo -is_oblique %s' % (getdsname(0))).readlines()[0].strip()) 319if oblique_epi_read or oblique_anat_read: 320 oblique_mode = True 321 sl.append("echo Oblique data detected.") 322else: oblique_mode = False 323if options.fres: 324 if options.qwarp: qwfres="-dxyz %s" % options.fres 325 else: alfres = "-mast_dxyz %s" % options.fres 326else: 327 if options.qwarp: qwfres="-dxyz ${voxsize}" #See section called "Preparing functional masking for this ME-EPI run" 328 else: alfres="-mast_dxyz ${voxsize}" 329if options.anat=='' and options.mask_mode!='func': 330 print "*+ Can't do anatomical-based functional masking without an anatomical!" 331 sys.exit() 332 333#Detect if current AFNI has old 3dNwarpApply 334if " -affter aaa = *** THIS OPTION IS NO LONGER AVAILABLE" in commands.getstatusoutput("3dNwarpApply -help")[1]: old_qwarp = False 335else: old_warp = True 336 337#Detect AFNI direcotry 338afnidir = os.path.dirname(os.popen('which 3dSkullStrip').readlines()[0]) 339 340#Prepare script and enter MEICA directory 341logcomment("Set up script run environment",level=1) 342sl.append('set -e') 343sl.append('export OMP_NUM_THREADS=%s' % (options.cpus)) 344sl.append('export MKL_NUM_THREADS=%s' % (options.cpus)) 345sl.append('export DYLD_FALLBACK_LIBRARY_PATH=%s' % (afnidir)) 346sl.append('export AFNI_3dDespike_NEW=YES') 347if options.overwrite: 348 sl.append('rm -rf meica.%s' % (setname)) 349else: 350 sl.append("if [[ -e meica.%s ]]; then echo ME-ICA directory exists, exiting; exit; fi" % (setname)) 351sl.append('mkdir -p meica.%s' % (setname)) 352sl.append("cp _meica_%s.sh meica.%s/" % (setname,setname)) 353sl.append("cd meica.%s" % setname) 354thecwd= "%s/meica.%s" % (getcwd(),setname) 355 356ica_datasets = sorted(datasets) 357 358#Parse anatomical processing options, process anatomical 359if options.anat != '': 360 logcomment("Deoblique, unifize, skullstrip, and/or autobox anatomical, in starting directory (may take a little while)", level=1) 361 nsmprage = options.anat 362 anatprefix=dsprefix(nsmprage) 363 pathanatprefix="%s/%s" % (startdir,anatprefix) 364 if oblique_mode: 365 sl.append("if [ ! -e %s_do.nii.gz ]; then 3dWarp -overwrite -prefix %s_do.nii.gz -deoblique %s/%s; fi" % (pathanatprefix,pathanatprefix,startdir,nsmprage)) 366 nsmprage="%s_do.nii.gz" % (anatprefix) 367 if not options.no_skullstrip: 368 sl.append("if [ ! -e %s_ns.nii.gz ]; then 3dUnifize -overwrite -prefix %s_u.nii.gz %s/%s; 3dSkullStrip -shrink_fac_bot_lim 0.3 -orig_vol -overwrite -prefix %s_ns.nii.gz -input %s_u.nii.gz; 3dAutobox -overwrite -prefix %s_ns.nii.gz %s_ns.nii.gz; fi" % (pathanatprefix,pathanatprefix,startdir,nsmprage,pathanatprefix,pathanatprefix,pathanatprefix,pathanatprefix)) 369 nsmprage="%s_ns.nii.gz" % (anatprefix) 370 371# Copy in functional datasets as NIFTI (if not in NIFTI already), calculate rigid body alignment 372vrbase=getdsname(0,True) 373logcomment("Copy in functional datasets, reset NIFTI tags as needed", level=1) 374for e_ii in range(len(datasets)): 375 ds = datasets[e_ii] 376 sl.append("3dcalc -a %s/%s -expr 'a' -prefix ./%s.nii" % (startdir,getdsname(e_ii),getdsname(e_ii,True) ) ) 377 if '.nii' in isf: 378 sl.append("nifti_tool -mod_hdr -mod_field sform_code 1 -mod_field qform_code 1 -infiles ./%s.nii -overwrite" % ( getdsname(e_ii,True) )) 379isf = '.nii' 380 381logcomment("Calculate and save motion and obliquity parameters, despiking first if not disabled, and separately save and mask the base volume",level=1) 382#Determine input to volume registration 383vrAinput = "./%s%s" % (vrbase,isf) 384#Compute obliquity matrix 385if oblique_mode: 386 if options.anat!='': sl.append("3dWarp -verb -card2oblique %s[0] -overwrite -newgrid 1.000000 -prefix ./%s_ob.nii.gz %s/%s | \grep -A 4 '# mat44 Obliquity Transformation ::' > %s_obla2e_mat.1D" % (vrAinput,anatprefix,startdir,nsmprage,prefix)) 387 else: sl.append("3dWarp -overwrite -prefix %s -deoblique %s" % (vrAinput,vrAinput)) 388#Despike and axialize 389if not options.no_despike: 390 sl.append("3dDespike -overwrite -prefix ./%s_vrA%s %s " % (vrbase,osf,vrAinput)) 391 vrAinput = "./%s_vrA%s" % (vrbase,osf) 392if not options.no_axialize: 393 sl.append("3daxialize -overwrite -prefix ./%s_vrA%s %s" % (vrbase,osf,vrAinput)) 394 vrAinput = "./%s_vrA%s" % (vrbase,osf) 395#Set eBbase 396external_eBbase=False 397if options.align_base!='': 398 if options.align_base.isdigit(): 399 basevol = '%s[%s]' % (vrAinput,options.align_base) 400 else: 401 basevol = options.align_base 402 external_eBbase=True 403else: 404 basevol = '%s[%s]' % (vrAinput,basebrik) 405sl.append("3dcalc -a %s -expr 'a' -prefix eBbase.nii.gz " % (basevol) ) 406if external_eBbase: 407 if oblique_mode: sl.append("3dWarp -overwrite -deoblique eBbase.nii.gz eBbase.nii.gz") 408 if not options.no_axialize: sl.append("3daxialize -overwrite -prefix eBbase.nii.gz eBbase.nii.gz") 409#Compute motion parameters 410sl.append("3dvolreg -overwrite -tshift -quintic -prefix ./%s_vrA%s -base eBbase.nii.gz -dfile ./%s_vrA.1D -1Dmatrix_save ./%s_vrmat.aff12.1D %s" % \ 411 (vrbase,osf,vrbase,prefix,vrAinput)) 412vrAinput = "./%s_vrA%s" % (vrbase,osf) 413sl.append("1dcat './%s_vrA.1D[1..6]{%s..$}' > motion.1D " % (vrbase,basebrik)) 414e2dsin = prefix+datasets[0]+trailing 415 416logcomment("Preliminary preprocessing of functional datasets: despike, tshift, deoblique, and/or axialize",level=1) 417#Do preliminary preproc for this run 418if shorthand_dsin: datasets.sort() 419for echo_ii in range(len(datasets)): 420 #Determine dataset name 421 echo = datasets[echo_ii] 422 indata = getdsname(echo_ii) 423 dsin = 'e'+echo 424 if echo_ii==0: e1_dsin = dsin 425 logcomment("Preliminary preprocessing dataset %s of TE=%sms to produce %s_ts+orig" % (indata,str(tes[echo_ii]),dsin) ) 426 #Pre-treat datasets: De-spike, RETROICOR in the future? 427 intsname = '%s%s' % (dsprefix(indata),isf) 428 if not options.no_despike: 429 intsname = "./%s_pt.nii.gz" % dsprefix(indata) 430 sl.append("3dDespike -overwrite -prefix %s %s%s" % (intsname,dsprefix(indata),isf)) 431 #Time shift datasets 432 if options.tpattern!='': 433 tpat_opt = ' -tpattern %s ' % options.tpattern 434 else: 435 tpat_opt = '' 436 sl.append("3dTshift -heptic %s -prefix ./%s_ts+orig %s" % (tpat_opt,dsin,intsname) ) 437 #Force +orig label on dataset 438 sl.append("3drefit -view orig %s_ts*HEAD" % (dsin)) 439 if oblique_mode and options.anat=="": 440 sl.append("3dWarp -overwrite -deoblique -prefix ./%s_ts+orig ./%s_ts+orig" % (dsin,dsin)) 441 #Axialize functional dataset 442 if not options.no_axialize: 443 sl.append("3daxialize -overwrite -prefix ./%s_ts+orig ./%s_ts+orig" % (dsin,dsin)) 444 if oblique_mode: sl.append("3drefit -deoblique -TR %s %s_ts+orig" % (options.TR,dsin)) 445 else: sl.append("3drefit -TR %s %s_ts+orig" % (options.TR,dsin)) 446 447#Compute T2*, S0, and OC volumes from raw data 448logcomment("Prepare T2* and S0 volumes for use in functional masking and (optionally) anatomical-functional coregistration (takes a little while).",level=1) 449dss = datasets 450dss.sort() 451stackline="" 452for echo_ii in range(len(dss)): 453 echo = datasets[echo_ii] 454 dsin = 'e'+echo 455 sl.append("3dAllineate -overwrite -final NN -NN -float -1Dmatrix_apply %s_vrmat.aff12.1D'{%i..%i}' -base eBbase.nii.gz -input %s_ts+orig'[%i..%i]' -prefix %s_vrA.nii.gz" % \ 456 (prefix,int(basebrik),int(basebrik)+5,dsin,int(basebrik),int(basebrik)+5,dsin)) 457 stackline+=" %s_vrA.nii.gz" % (dsin) 458sl.append("3dZcat -prefix basestack.nii.gz %s" % (stackline)) 459sl.append("%s %s -d basestack.nii.gz -e %s" % (sys.executable, '/'.join([meicadir,'meica.libs','t2smap.py']),options.tes)) 460sl.append("3dUnifize -prefix ./ocv_uni+orig ocv.nii") 461sl.append("3dSkullStrip -prefix ./ocv_ss.nii.gz -overwrite -input ocv_uni+orig") 462sl.append("3dcalc -overwrite -a t2svm.nii -b ocv_ss.nii.gz -expr 'a*ispositive(a)*step(b)' -prefix t2svm_ss.nii.gz" ) 463sl.append("3dcalc -overwrite -a s0v.nii -b ocv_ss.nii.gz -expr 'a*ispositive(a)*step(b)' -prefix s0v_ss.nii.gz" ) 464if not options.no_axialize: 465 sl.append("3daxialize -overwrite -prefix t2svm_ss.nii.gz t2svm_ss.nii.gz") 466 sl.append("3daxialize -overwrite -prefix ocv_ss.nii.gz ocv_ss.nii.gz") 467 sl.append("3daxialize -overwrite -prefix s0v_ss.nii.gz s0v_ss.nii.gz") 468 469# Calculate affine anatomical warp if anatomical provided, then combine motion correction and coregistration parameters 470if options.anat!='': 471 #Copy in anatomical and make sure its in +orig space 472 logcomment("Copy anatomical into ME-ICA directory and process warps",level=1) 473 sl.append("cp %s/%s* ." % (startdir,nsmprage)) 474 abmprage = nsmprage 475 refanat = nsmprage 476 if options.space: 477 sl.append("afnibinloc=`which 3dSkullStrip`") 478 if '/' in options.space: 479 sl.append("ll=\"%s\"; templateloc=${ll%%/*}/" % options.space) 480 options.space=options.space.split('/')[-1] 481 else: 482 sl.append("templateloc=${afnibinloc%/*}") 483 atnsmprage = "%s_at.nii.gz" % (dsprefix(nsmprage)) 484 if not dssuffix(nsmprage).__contains__('nii'): sl.append("3dcalc -float -a %s -expr 'a' -prefix %s.nii.gz" % (nsmprage,dsprefix(nsmprage))) 485 logcomment("If can't find affine-warped anatomical, copy native anatomical here, compute warps (takes a while) and save in start dir. ; otherwise link in existing files") 486 sl.append("if [ ! -e %s/%s ]; then \@auto_tlrc -no_ss -init_xform AUTO_CENTER -base ${templateloc}/%s -input %s.nii.gz -suffix _at" % (startdir,atnsmprage,options.space,dsprefix(nsmprage))) 487 sl.append("cp %s.nii %s" % (dsprefix(atnsmprage),startdir)) 488 sl.append("gzip -f %s/%s.nii" % (startdir,dsprefix(atnsmprage))) 489 sl.append("else ln -s %s/%s ." % (startdir,atnsmprage)) 490 refanat = '%s/%s' % (startdir,atnsmprage) 491 sl.append("fi") 492 sl.append("3dcopy %s/%s.nii.gz %s" % (startdir,dsprefix(atnsmprage),dsprefix(atnsmprage))) 493 sl.append("3drefit -view orig %s+tlrc " % dsprefix(atnsmprage) ) 494 sl.append("3dAutobox -prefix ./abtemplate.nii.gz ${templateloc}/%s" % options.space) 495 abmprage = 'abtemplate.nii.gz' 496 if options.qwarp: 497 logcomment("If can't find non-linearly warped anatomical, compute, save back; otherwise link") 498 nlatnsmprage="%s_atnl.nii.gz" % (dsprefix(nsmprage)) 499 sl.append("if [ ! -e %s/%s ]; then " % (startdir,nlatnsmprage)) 500 logcomment("Compute non-linear warp to standard space using 3dQwarp (get lunch, takes a while) ") 501 sl.append("3dUnifize -overwrite -GM -prefix ./%su.nii.gz %s/%s" % (dsprefix(atnsmprage),startdir,atnsmprage)) 502 sl.append("3dQwarp -iwarp -overwrite -resample -useweight -blur 2 2 -duplo -workhard -base ${templateloc}/%s -prefix %s/%snl.nii.gz -source ./%su.nii.gz" % (options.space,startdir,dsprefix(atnsmprage),dsprefix(atnsmprage))) 503 sl.append("fi") 504 sl.append("ln -s %s/%s ." % (startdir,nlatnsmprage)) 505 refanat = '%s/%snl.nii.gz' % (startdir,dsprefix(atnsmprage)) 506 507 #Set anatomical reference for anatomical-functional co-registration 508 if oblique_mode: alnsmprage = "./%s_ob.nii.gz" % (anatprefix) 509 else: alnsmprage = "%s/%s" % (startdir,nsmprage) 510 if options.coreg_mode=='lp-t2s': 511 logcomment("Using alignp_mepi_anat.py to drive T2*-map weighted anatomical-functional coregistration") 512 ama_alnsmprage = alnsmprage 513 if not options.no_axialize: 514 ama_alnsmprage = os.path.basename(alnsmprage) 515 sl.append("3daxialize -overwrite -prefix ./%s %s" % (ama_alnsmprage,alnsmprage)) 516 t2salignpath = 'meica.libs/alignp_mepi_anat.py' 517 sl.append("%s %s -t t2svm_ss.nii.gz -a %s -p mepi %s" % \ 518 (sys.executable, '/'.join([meicadir,t2salignpath]),ama_alnsmprage,options.align_args)) 519 sl.append("cp alignp.mepi/mepi_al_mat.aff12.1D ./%s_al_mat.aff12.1D" % anatprefix) 520 elif options.coreg_mode=='aea': 521 logcomment("Using AFNI align_epi_anat.py to drive anatomical-functional coregistration ") 522 sl.append("3dcopy %s ./ANAT_ns+orig " % alnsmprage) 523 sl.append("align_epi_anat.py -anat2epi -volreg off -tshift off -deoblique off -anat_has_skull no -save_script aea_anat_to_ocv.tcsh -anat ANAT_ns+orig -epi ocv_uni+orig -epi_base 0 %s" % (options.align_args) ) 524 sl.append("cp ANAT_ns_al_mat.aff12.1D %s_al_mat.aff12.1D" % (anatprefix)) 525 if options.space: 526 tlrc_opt = "%s/%s::WARP_DATA -I" % (startdir,atnsmprage) 527 sl.append("cat_matvec -ONELINE %s > %s/%s_ns2at.aff12.1D" % (tlrc_opt,startdir,anatprefix)) 528 else: tlrc_opt = "" 529 if oblique_mode: oblique_opt = "%s_obla2e_mat.1D" % prefix 530 else: oblique_opt = "" 531 sl.append("cat_matvec -ONELINE %s %s %s_al_mat.aff12.1D -I > %s_wmat.aff12.1D" % (tlrc_opt,oblique_opt,anatprefix,prefix)) 532 sl.append("cat_matvec -ONELINE %s %s %s_al_mat.aff12.1D -I %s_vrmat.aff12.1D > %s_vrwmat.aff12.1D" % (tlrc_opt,oblique_opt,anatprefix,prefix,prefix)) 533 534else: sl.append("cp %s_vrmat.aff12.1D %s_vrwmat.aff12.1D" % (prefix,prefix)) 535 536#Preprocess datasets 537if shorthand_dsin: datasets.sort() 538logcomment("Extended preprocessing of functional datasets",level=1) 539 540#Compute grand mean scaling factor 541sl.append("3dBrickStat -mask eBbase.nii.gz -percentile 50 1 50 %s_ts+orig[%i] > gms.1D" % (e1_dsin,basebrik)) 542sl.append("gms=`cat gms.1D`; gmsa=($gms); p50=${gmsa[1]}") 543 544for echo_ii in range(len(datasets)): 545 546 #Determine dataset name 547 echo = datasets[echo_ii] 548 indata = getdsname(echo_ii) 549 dsin = 'e'+echo 550 551 if echo_ii == 0: 552 logcomment("Preparing functional masking for this ME-EPI run",2 ) 553 if options.anat: almaster="-master %s" % abmprage 554 else: almaster="" 555 sl.append("3dZeropad %s -prefix eBvrmask.nii.gz ocv_ss.nii.gz[0]" % (zeropad_opts)) 556 sl.append("voxsize=`ccalc $(3dinfo -voxvol eBvrmask.nii.gz)**.33`") #Set voxel size 557 #Create base mask 558 if options.anat and options.space and options.qwarp: 559 if old_qwarp: affter_string = "'-affter '%s_wmat.aff12.1D'" % prefix 560 else: affter_string = "" 561 sl.append("3dNwarpApply -overwrite -nwarp '%s/%s_WARP.nii.gz' %s %s -source eBvrmask.nii.gz -interp %s -prefix ./eBvrmask.nii.gz " % \ 562 (startdir,dsprefix(nlatnsmprage),prefix,almaster,qwfres,'NN')) 563 if options.t2salign or options.mask_mode!='func': 564 sl.append("3dNwarpApply -overwrite -nwarp '%s/%s_WARP.nii.gz' %s %s %s -source t2svm_ss.nii.gz -interp %s -prefix ./t2svm_ss_vr.nii.gz " % \ 565 (startdir,dsprefix(nlatnsmprage),affter_string,almaster,qwfres,'NN')) 566 sl.append("3dNwarpApply -overwrite -nwarp '%s/%s_WARP.nii.gz' %s %s %s -source ocv_uni+orig -interp %s -prefix ./ocv_uni_vr.nii.gz " % \ 567 (startdir,dsprefix(nlatnsmprage),affter_string,almaster,qwfres,'NN')) 568 sl.append("3dNwarpApply -overwrite -nwarp '%s/%s_WARP.nii.gz' %s %s %s -source s0v_ss.nii.gz -interp %s -prefix ./s0v_ss_vr.nii.gz " % \ 569 (startdir,dsprefix(nlatnsmprage),affter_string,almaster,qwfres,'NN')) 570 elif options.anat: 571 sl.append("3dAllineate -overwrite -final %s -%s -float -1Dmatrix_apply %s_wmat.aff12.1D -base eBvrmask.nii.gz -input eBvrmask.nii.gz -prefix ./eBvrmask.nii.gz %s %s" % \ 572 ('NN','NN',prefix,almaster,alfres)) 573 if options.t2salign or options.mask_mode!='func': 574 sl.append("3dAllineate -overwrite -final %s -%s -float -1Dmatrix_apply %s_wmat.aff12.1D -base eBvrmask.nii.gz -input t2svm_ss.nii.gz -prefix ./t2svm_ss_vr.nii.gz %s %s" % \ 575 ('NN','NN',prefix,almaster,alfres)) 576 sl.append("3dAllineate -overwrite -final %s -%s -float -1Dmatrix_apply %s_wmat.aff12.1D -base eBvrmask.nii.gz -input ocv_uni+orig -prefix ./ocv_uni_vr.nii.gz %s %s" % \ 577 ('NN','NN',prefix,almaster,alfres)) 578 sl.append("3dAllineate -overwrite -final %s -%s -float -1Dmatrix_apply %s_wmat.aff12.1D -base eBvrmask.nii.gz -input s0v_ss.nii.gz -prefix ./s0v_ss_vr.nii.gz %s %s" % \ 579 ('NN','NN',prefix,almaster,alfres)) 580 581 if options.anat and options.mask_mode != 'func': 582 if options.space and options.mask_mode == 'template': 583 sl.append("3dfractionize -template eBvrmask.nii.gz -input abtemplate.nii.gz -prefix ./anatmask_epi.nii.gz -clip 1") 584 logcomment("Preparing functional mask using information from standard space template (takes a little while)") 585 if options.mask_mode == 'anat': 586 sl.append("3dfractionize -template eBvrmask.nii.gz -input %s -prefix ./anatmask_epi.nii.gz -clip 0.5" % (refanat) ) 587 logcomment("Preparing functional mask using information from anatomical (takes a little while)") 588 sl.append("3dBrickStat -mask eBvrmask.nii.gz -percentile 50 1 50 t2svm_ss_vr.nii.gz > t2s_med.1D") 589 sl.append("3dBrickStat -mask eBvrmask.nii.gz -percentile 50 1 50 s0v_ss_vr.nii.gz > s0v_med.1D") 590 sl.append("t2sm=`cat t2s_med.1D`; t2sma=($t2sm); t2sm=${t2sma[1]}") 591 sl.append("s0vm=`cat s0v_med.1D`; s0vma=($s0vm); s0vm=${s0vma[1]}") 592 sl.append("3dcalc -a ocv_uni_vr.nii.gz -b anatmask_epi.nii.gz -c t2svm_ss_vr.nii.gz -d s0v_ss_vr.nii.gz -expr \"a-a*equals(equals(b,0)+isnegative(c-${t2sm})+ispositive(d-${s0vm}),3)\" -overwrite -prefix ocv_uni_vr.nii.gz ") 593 sl.append("3dSkullStrip -overwrite -input ocv_uni_vr.nii.gz -prefix eBvrmask.nii.gz ") 594 if options.fres: resstring = "-dxyz %s %s %s" % (options.fres,options.fres,options.fres) 595 else: resstring = "-dxyz ${voxsize} ${voxsize} ${voxsize}" 596 sl.append("3dresample -overwrite -master %s %s -input eBvrmask.nii.gz -prefix eBvrmask.nii.gz" % (abmprage,resstring)) 597 598 if options.anat=='': 599 logcomment("Trim empty space off of mask dataset and/or resample") 600 sl.append("3dAutobox -overwrite -prefix eBvrmask%s eBvrmask%s" % (osf,osf) ) 601 if options.fres: 602 resstring = "-dxyz %s %s %s" % (options.fres,options.fres,options.fres) 603 sl.append("3dresample -overwrite -master eBvrmask.nii.gz %s -input eBvrmask.nii.gz -prefix eBvrmask.nii.gz" % (resstring)) 604 605 sl.append("3dcalc -float -a eBvrmask.nii.gz -expr 'notzero(a)' -overwrite -prefix eBvrmask.nii.gz") 606 607 #logcomment("Extended preprocessing dataset %s of TE=%sms to produce %s_in.nii.gz" % (indata,str(tes[echo_ii]),dsin),level=2 ) 608 logcomment("Apply combined normalization/co-registration/motion correction parameter set to %s_ts+orig" % dsin) 609 if options.qwarp: 610 if old_qwarp: affter_string = "-affter %s_vrwmat.aff12.1D" % prefix 611 else: affter_string = "" 612 sl.append("3dNwarpApply -nwarp '%s/%s_WARP.nii.gz' %s -master eBvrmask.nii.gz -source %s_ts+orig -interp %s -prefix ./%s_vr%s " % \ 613 (startdir,dsprefix(nlatnsmprage),affter_string,dsin,align_interp,dsin,osf)) 614 else: sl.append("3dAllineate -final %s -%s -float -1Dmatrix_apply %s_vrwmat.aff12.1D -base eBvrmask%s -input %s_ts+orig -prefix ./%s_vr%s" % \ 615 (align_interp_final,align_interp,prefix,osf,dsin,dsin,osf)) 616 if echo_ii == 0: 617 sl.append("3dTstat -min -prefix ./%s_vr_min%s ./%s_vr%s" % (dsin,osf,dsin,osf) ) 618 sl.append("3dcalc -a eBvrmask.nii.gz -b %s_vr_min%s -expr 'step(a)*step(b)' -overwrite -prefix eBvrmask.nii.gz " % (dsin,osf)) 619 if options.FWHM=='0mm': 620 sl.append("3dcalc -float -overwrite -a eBvrmask.nii.gz -b ./%s_vr%s[%i..$] -expr 'step(a)*b' -prefix ./%s_sm%s " % (dsin,osf,basebrik,dsin,osf)) 621 else: 622 sl.append("3dBlurInMask -fwhm %s -mask eBvrmask%s -prefix ./%s_sm%s ./%s_vr%s[%i..$]" % (options.FWHM,osf,dsin,osf,dsin,osf,basebrik)) 623 sl.append("3dcalc -float -overwrite -a ./%s_sm%s -expr \"a*10000/${p50}\" -prefix ./%s_sm%s" % (dsin,osf,dsin,osf)) 624 sl.append("3dTstat -prefix ./%s_mean%s ./%s_sm%s" % (dsin,osf,dsin,osf)) 625 if options.detrend: sl.append("3dDetrend -polort %s -overwrite -prefix ./%s_sm%s ./%s_sm%s " % (options.detrend,dsin,osf,dsin,osf) ) 626 if options.highpass: sl.append("3dBandpass -prefix ./%s_in%s %f 99 ./%s_sm%s " % (dsin,osf,float(options.highpass),dsin,osf) ) 627 else: sl.append("mv %s_sm%s %s_in%s" % (dsin,osf,dsin,osf)) 628 sl.append("3dcalc -float -overwrite -a ./%s_in%s -b ./%s_mean%s -expr 'a+b' -prefix ./%s_in%s" % (dsin,osf,dsin,osf,dsin,osf)) 629 sl.append("3dTstat -stdev -prefix ./%s_std%s ./%s_in%s" % (dsin,osf,dsin,osf)) 630 if options.test_proc: sl.append("exit") 631 if not (options.test_proc or options.keep_int): sl.append("rm -f %s_ts+orig* %s_vr%s %s_sm%s" % (dsin,dsin,osf,dsin,osf)) 632 633#Concatenate for ICA 634if len(ica_datasets)==1: 635 dsin = ''.join(ica_datasets)+trailing 636 ica_prefix = dsin 637 ica_input="./%s_in%s" % (prefix,''.join(ica_datasets),trailing) 638 ica_mask="eBvrmask.nii.gz" 639else: 640 ica_input = "zcat_ffd.nii.gz" 641 ica_mask = "zcat_mask.nii.gz" 642 zcatstring="" 643 for echo in ica_datasets: 644 dsin ='e'+echo 645 zcatstring = "%s ./%s_in%s" % (zcatstring,dsin,osf) 646 sl.append("3dZcat -overwrite -prefix %s %s" % (ica_input,zcatstring) ) 647 sl.append("3dcalc -float -overwrite -a %s[0] -expr 'notzero(a)' -prefix %s" % (ica_input,ica_mask)) 648 649if options.pp_only: tedflag='#' 650else: tedflag = '' 651 652if os.path.exists('%s/meica.libs' % (meicadir)): tedanapath = 'meica.libs/tedana.py' 653else: tedanapath = 'tedana.py' 654logcomment("Perform TE-dependence analysis (takes a good while)",level=1) 655sl.append("%s%s %s -e %s -d %s --sourceTEs=%s --kdaw=%s --rdaw=1 --initcost=%s --finalcost=%s --conv=2.5e-5 %s" % (tedflag,sys.executable, '/'.join([meicadir,tedanapath]),options.tes,ica_input,options.sourceTEs,options.daw,options.initcost,options.finalcost,options.ted_args)) 656sl.append("#") 657if outprefix=='': outprefix=setname 658 659logcomment("Copying results to start directory",level=1) 660 661sl.append("%scp TED/ts_OC.nii TED/%s_tsoc.nii" % (tedflag,outprefix)) 662sl.append("%scp TED/dn_ts_OC.nii TED/%s_medn.nii" % (tedflag,outprefix)) 663sl.append("%scp TED/betas_hik_OC.nii TED/%s_mefc.nii" % (tedflag,outprefix)) 664sl.append("%scp TED/betas_OC.nii TED/%s_mefl.nii" % (tedflag,outprefix)) 665sl.append("%scp TED/comp_table.txt %s/%s_ctab.txt" % (tedflag,startdir,outprefix)) 666hist_line = "%s" % (" ".join(sys.argv).replace('"',r"\"")) 667note_line = "Denoised timeseries (including thermal noise), produced by ME-ICA v2.5" 668sl.append("%s3dNotes -h \'%s (%s)\' TED/%s_medn.nii" % (tedflag,hist_line,note_line,outprefix)) 669note_line = "Denoised ICA coeff. set for ME-ICR seed-based FC analysis, produced by ME-ICA v2.5" 670sl.append("%s3dNotes -h \'%s (%s)\' TED/%s_mefc.nii" % (tedflag,hist_line,note_line,outprefix)) 671note_line = "Full ICA coeff. set for component assessment, produced by ME-ICA v2.5" 672sl.append("%s3dNotes -h \'%s (%s)\' TED/%s_mefc.nii" % (tedflag,hist_line,note_line,outprefix)) 673note_line = "T2* weighted average of ME time series, produced by ME-ICA v2.5" 674sl.append("%s3dNotes -h \'%s (%s)\' TED/%s_tsoc.nii" % (tedflag,hist_line,note_line,outprefix)) 675if options.anat!='' and options.space!=False: 676 sl.append("%snifti_tool -mod_hdr -mod_field sform_code 2 -mod_field qform_code 2 -infiles TED/%s_tsoc.nii -overwrite" % (tedflag,outprefix)) 677 sl.append("%snifti_tool -mod_hdr -mod_field sform_code 2 -mod_field qform_code 2 -infiles TED/%s_medn.nii -overwrite" % (tedflag,outprefix)) 678 sl.append("%snifti_tool -mod_hdr -mod_field sform_code 2 -mod_field qform_code 2 -infiles TED/%s_mefc.nii -overwrite" % (tedflag,outprefix)) 679 sl.append("%snifti_tool -mod_hdr -mod_field sform_code 2 -mod_field qform_code 2 -infiles TED/%s_mefl.nii -overwrite" % (tedflag,outprefix)) 680sl.append("%sgzip -f TED/%s_medn.nii TED/%s_mefc.nii TED/%s_tsoc.nii TED/%s_mefl.nii" % (tedflag,outprefix,outprefix,outprefix,outprefix)) 681sl.append("%smv TED/%s_medn.nii.gz TED/%s_mefc.nii.gz TED/%s_tsoc.nii.gz TED/%s_mefl.nii.gz %s" % (tedflag,outprefix,outprefix,outprefix,outprefix,startdir)) 682 683 684#Write the preproc script and execute it 685ofh = open('_meica_%s.sh' % setname ,'w') 686print "++ Writing script file: _meica_%s.sh" % (setname) 687ofh.write("\n".join(sl)+"\n") 688ofh.close() 689if not options.script_only: 690 print "++ Executing script file: _meica_%s.sh" % (setname) 691 system('bash _meica_%s.sh' % setname) 692