ANTsR Notes: 8/4/2016

 General Notes: You will want to make .nii.gz files for use before running this. This underlined steps need to be modified for each task other commands are generalizable. If you copy the table below into excel you should be able to modify the second column and directly paste into R. Later I may source this file but for now I prefer being aware of all substeps listed.

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=""))}

  # slice time has been pulled form ANTsR so use system to run AFNI
functional slice time if needed
        if(!(file.exists(paste(study,"/",sub,"/",sub,task,"_dst.nii.gz",sep=""))){system(paste("3dTshift -ignore 2 -tzero 0 -prefix ",sub,task,"_dst.nii.gz ",sub,task,"._ds.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)

}