Comment

Command

load library

library(ANTsR)

file to find (raw)

task<-"MR"

file to save (output)

taskout<-"MR"

regressor formula

lmformula<-"boldmatp[Censor==1,]~regc$EasyC+regc$MedC+regc$HardC+regc$EasyE+regc$MedE+regc$HardE"

# This is where you decide what regressors you care about

study location

study<-"/mnt/nfs2/users/alans/Config/ANTS"

move to study

setwd(study)

make subjectlist

subs<-Sys.glob("Config_???")

walk through list

for (sub in subs){

if done skip

    if(!(file.exists(paste(study,"/",sub,"/",sub,taskout,"_ANTS.nii.gz",sep="")))){

if raw exists continue

      if(file.exists(paste(study,"/",sub,"/",sub,task,".nii.gz",sep=""))){

move to subject directory

        setwd(paste(study,sub,sep="/"))

load generic regressor if available

        if(file.exists(paste(task,"Rwav",sep=""))){reg<-read.delim(paste(task,"Rwav",sep=""),header=TRUE)}

load personal regressor if available

        if(file.exists(paste(sub,task,"Rwav",sep=""))){reg<-read.delim(paste(sub,task,"Rwav",sep=""),header=TRUE)}

make OL file

        system(paste("3dToutcount -automask ",sub,task,".nii.gz > ",sub,task,"_OL.txt",sep=""))

# for now I'm optioning out data quality check to AFNI

make censor file

        system(paste("PkCensorSD ",sub,task,"_OL.txt 1.5 ",sub,task,"Censor",sep=""))

# for now I'm creating Censors based on the AFNI OL file

load censor file

        Censor<-read.delim(paste(sub,task,"Censor",sep=""),header=FALSE)

load T1

        t1<-antsImageRead(paste(sub,".nii.gz",sep=""),3)

clean T1 (field correct)

        t1c<-abpN4(t1)

save clean T1

        antsImageWrite(t1c,paste(sub,"_c.nii.gz",sep=""))

load target brain

        tem<-antsImageRead('/mnt/nfs2/users/alans/niiTemplates/mni_icbm152_nlin_sym_09c/mni_icbm152_t1_tal_nlin_sym_09c.nii.gz')

load target brain mask

        temmask<-antsImageRead('/mnt/nfs2/users/alans/niiTemplates/mni_icbm152_nlin_sym_09c/mni_icbm152_t1_tal_nlin_sym_09c_mask.nii.gz')

strip brain

        bm<-abpBrainExtraction(t1c,tem,temmask)

move T1 to target

        bmbrain <- antsApplyTransforms(fixed=tem,moving=bm$brain,transformlist=bm$fwdtransforms )

save MNI T1

        antsImageWrite(bmbrain,paste(sub,"_t1.nii.gz",sep=""))

save morph nifti

        file.copy(bm$fwdtransforms[1],paste(sub,"_ft.nii.gz",sep=""))

save affine mat

        file.copy(bm$fwdtransforms[2],paste(sub,"_ft.mat",sep=""))

move segmented T1 to target

        bmseg <- antsApplyTransforms(fixed=tem,moving=bm$kmeansseg,transformlist=bm$fwdtransforms )

save MNI segment

        antsImageWrite(bmseg,paste(sub,"_seg.nii.gz",sep=""))

functional despike if needed

        if(!(file.exists(paste(study,"/",sub,"/",sub,task,"_ds.nii.gz",sep="")))){system(paste("3dDespike -ignore 4 -prefix ",sub,task,"_ds.nii.gz ",sub,task,".nii.gz",sep=""))}

# while R offers really nice methods for despiking it takes about 20 minutes versus a few seconds in AFNI

load despiked function

        bold<-antsImageRead(paste(sub,task,"_ds.nii.gz",sep=""),4)

calculate bold brain

        avg<-getAverageOfTimeSeries(bold)

# this creates the basis of you bold mask and alignment. If you have excessive movement before this it can cause some problems

field correct avg bold

        avg<-n3BiasFieldCorrection(avg,1)

# this correction allows for better alignment

make mask from avg bold

        mask<-getMask(avg,mean(avg),Inf,2)

# this mask will allow transfer from multidimension ANTs files to 2 dimensional matrix

correct for slice timing

        boldst<-

# adjust for slice timing

convert to matrix

        boldstm<-timeseries2matrix(boldst,mask)

# some functions work with matricies so you can flip them

whiten bold

        boldstmw<-temporalwhiten(boldstm, myord = 2)

# this whitens the data (aka fsl or afni:3dREML)

convert to ANTs file

        boldstw<-matrix2timeseries(bold,mask,boldstmw)

# this converts the other way

preprocess brick

        boldpre<-preprocessfMRI( boldstw, mask, useMotionCorrectedImage = TRUE,spatialSmoothingType="perona-malik",spatialSmoothingParameters=c(0.25, 5))

# The basic principals of this preprocessing follows: Powers et al.
# does: compcor/motion correction, nuisance regression, band-pass filtering, and spatial smoothing. As requested.
# compcor is similar to retroicor but using ica for noise determination
# spacial smoothing can be important to reduce noise you can use an anisotropic approach: related example Nam et al.

calculate bold to raw T1

        boldtx<-antsRegistration(fixed=avg, moving=bm$brain, typeofTransform="SyNBold")

build full transform list

        fulltx<-c(bm$fwdtransforms[1],bm$fwdtransforms[2],boldtx$invtransforms[1],boldtx$invtransforms[2])

resample MNI for intended save

        temrs<-resampleImage(tem,c(4,4,4),interpType=4)

# the betas will be moved to the target matrix so we change it here

convert bold to matrix

        boldmat<-timeseries2matrix(boldpre$cleanBoldImage,boldpre$maskImage)

# move this to matrix

norm bold data

        boldmatp<-(100*boldmat)/colMeans(boldstm)

# norm bold signal

trim regressor to bold max

        regclip<-reg[1:dim(bold)[4],]

# just in case the regressor is longer than the bold

remove censored trials

        regc<-regclip[Censor==1,]

# censor the regressor for regression

run regressor

        mylm<-lm(lmformula)

# run regression — you can do a robust regression but it takes 30+ minutes and you can't use same steps below

sum lm for saving

        biglm<-bigLMStats(mylm)

# makes the lm results easy to save

walk through beta values

        for (i in 1:(dim(biglm$beta)[1])) {

 make beta ANTs file

         beta<-makeImage(boldpre$maskImage,biglm$beta[i,])

 move beta to MNI

         beta_al<-antsApplyTransforms(fixed=temrs, moving=beta, transformlist=fulltx, whichtoinvert=c(F,F,T,F) )

 save MNI beta

         antsImageWrite(beta_al,paste(sub,taskout,"_b",i,".nii.gz",sep=""))}

# this creates nifti files from betas

walk through T values

        for (i in 1:(dim(biglm$beta.t)[1])) {

 make T value ANTs file

         beta<-makeImage(boldpre$maskImage,biglm$beta.t[i,])

 move T value to MNI

         beta_al<-antsApplyTransforms(fixed=temrs, moving=beta, transformlist=fulltx, whichtoinvert=c(F,F,T,F) )

 save MNI beta

         antsImageWrite(beta_al,paste(sub,taskout,"_b",i,"t.nii.gz",sep=""))}

# this creates nifti files from tvalues

make bucket of beta/Tvalues

        system(paste("3dbucket -prefix ",sub,taskout,"_ANTS.nii.gz ",sub,taskout,"_b*.nii.gz",sep=""))

# combine them into one file

delete sub files

        system(paste("rm ",sub,taskout,"_b*.nii.gz",sep=""))

clear notes list

        refit<-""

walk through beta values

        for (i in 1:(dim(biglm$beta)[1])) {

add name of beta to notes list

          refit_b <-paste(" -sublabel",2*(i-1),names(mylm$coefficients[,1])[i+1])

add name of T value to notes list

          refit_t <-paste(" -sublabel ",2*(i-1)+1," t_",names(mylm$coefficients[,1])[i+1],sep="")

combine notes list

          refit<-paste(refit,refit_b,refit_t)}

# this reads the lm for the name of the saved files

clear out reserved words

        refit<-gsub("sub","-sub",gsub("regc","",gsub("[[:punct:]]","",refit)))

label bucket file and specificy space

        system(paste("3drefit -space MNI -view tlrc ",refit," ",sub,taskout,"_ANTS.nii.gz",sep=""))

        system(paste("3drefit -space MNI -view tlrc ",sub,"_t1.nii.gz",sep=""))

end (if exists)

    }

end (if no raw)

  }

end (do next sub)

}