| Comment | Command |
| Load library | library(ANTsR) |
| specify task | tasks=c("MR290","MR") |
| setup study | study="/data/asimmons/data/CDA/ANTs/" |
| pull in brain in mni space | tem<-antsImageRead('/data/asimmons/data/niiTemplates/mni_icbm152_nlin_sym_09c/mni_icbm152_t1_tal_nlin_sym_09c.nii.gz') |
| pull in a matched mask in mni space | temmask<-antsImageRead('/data/asimmons/data/niiTemplates/mni_icbm152_nlin_sym_09c/mni_icbm152_t1_tal_nlin_sym_09c_mask.nii.gz') |
| move to study | setwd(study) |
| find all subjects | subs<-Sys.glob("TC_???_?") |
| loop through all subjects in directory | for (sub in subs){ |
| print to command line | print(paste("Checking and processing ",sub)) |
| if t1 doesn't exist make it | if(!(file.exists(paste(study,"/",sub,"/",sub,"_t1.nii.gz",sep="")))){ |
| move to subject directory | setwd(paste(study,sub,sep="/")) |
| load T1 | t1<-antsImageRead(paste(sub,".nii.gz",sep=""),3) |
| do Field Bias correction | t1c<-abpN4(t1) |
| save it | antsImageWrite(t1c,paste(sub,"_c.nii.gz",sep="")) |
| print to command line | print("Skull Strip") |
| do the skull stripping | bm<-abpBrainExtraction(t1c,tem,temmask) |
| move to mni space | bmbrain <- antsApplyTransforms(fixed=tem,moving=bm$brain,transformlist=bm$fwdtransforms ) |
| save it in orig space | antsImageWrite(bm$brain,paste(sub,"_t1orig.nii.gz",sep="")) |
| save it in mni space | antsImageWrite(bmbrain,paste(sub,"_t1.nii.gz",sep="")) |
| save in transforms | file.copy(bm$fwdtransforms[1],paste(sub,"_ft.nii.gz",sep="")) |
| save in transforms | file.copy(bm$fwdtransforms[2],paste(sub,"_ft.mat",sep="")) |
| save in transforms | file.copy(bm$invtransforms[1],paste(sub,"_invt.nii.gz",sep="")) |
| save in transforms | file.copy(bm$invtransforms[2],paste(sub,"_invt.mat",sep="")) |
| move seg file to mni | bmseg <- antsApplyTransforms(fixed=tem,moving=bm$kmeansseg,transformlist=bm$fwdtransforms ) |
| save seg file | antsImageWrite(bmseg,paste(sub,"_seg.nii.gz",sep=""))} |
| loop through all tasks | for (task in tasks) { |
| save an output name | taskout<-paste(task,"Prep",sep="") |
| if output doesn't exist continue | if(!(file.exists(paste(study,"/",sub,"/",sub,taskout,".nii.gz",sep="")))){ |
| if input exists continue | if(file.exists(paste(study,"/",sub,"/",sub,task,".nii.gz",sep=""))){ |
| move to subject directory | setwd(paste(study,sub,sep="/")) |
| run task | print(paste("running data for ",task," in ",sub,sep="")) |
| if not despiked then despike | if(!(file.exists(paste(sub,task,"_ds.nii.gz",sep="")))) |
| {system(paste("3dDespike -ignore 4 -prefix ",sub,task,"_ds.nii.gz ",sub,task,".nii.gz",sep=""))} |
| if not time shifted then time shift | if(!(file.exists(paste(sub,task,"_dst.nii.gz",sep="")))) |
| {system(paste("3dTshift -ignore 2 -tzero 0 -tpattern seqplus -prefix ",sub,task,"_dst.nii.gz ",sub,task,"_ds.nii.gz",sep=""))} |
| read in bold | bold=antsImageRead(paste(sub,task,"_dst.nii.gz",sep=""),4) |
| print("n3 of BOLD") |
| make average for alignment | avg <- getAverageOfTimeSeries(bold) |
| do Field Bias correction | avgn4 <- abpN4(avg) |
| make mask | mask<-getMask(avgn4,mean(avgn4),Inf,2) |
| do Field Bias correction | bold_n3<-timeseriesN3(bold,mask) |
| make vector mask for matrix conversions | avgm<-as.vector(avgn4[mask==1]) |
| convert to matrix | boldstm<-timeseries2matrix(bold_n3,mask) |
| convert to PSC | boldstmp<-100*t(t(boldstm)/avgm) |
| do temporal whitening | boldstmw<-temporalwhiten(boldstmp, myord = 2) |
| convert to timeseries | boldstw<-matrix2timeseries(bold,mask,boldstmw) |
| print("preprocess fMRI") |
| pull in t1 | bmbrainorig=antsImageRead(paste(sub,"_t1orig.nii.gz",sep="")) |
| pull in transforms | bmfwdtransforms1=paste(sub,"_ft.nii.gz",sep="") |
| pull in transforms | bmfwdtransforms2=paste(sub,"_ft.mat",sep="") |
| Do data cleaning | boldpre<-preprocessfMRI( boldstw, mask, useMotionCorrectedImage = TRUE) |
| # 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. |
| register data | boldtx<-antsRegistration(fixed=avg, moving=bmbrainorig, typeofTransform="SyNBold", spatialSmoothingType="perona-malik",spatialSmoothingParameters=c(0.25, 5)) |
| concat transforms | fulltx<-c(bmfwdtransforms1,bmfwdtransforms2,boldtx$invtransforms[1],boldtx$invtransforms[2]) |
| resample for mni outputs | temrs<-resampleImage(tem,c(4,4,4),interpType=4) |
| save output in orig space | antsImageWrite(boldpre$cleanBoldImage,paste(sub,taskout,".nii.gz",sep="")) |
| move to mni space | boldpreM<-antsApplyTransforms(fixed=temrs,moving=boldpre$cleanBoldImage,transformlist=fulltx,whichtoinvert=c(F,F,T,F),imagetype=3) |
| save output in mni space | antsImageWrite(boldpreM,paste(sub,taskout,"MNI.nii.gz",sep="")) |
| move mask to mni space | maskM<-antsApplyTransforms(fixed=temrs, moving=mask, transformlist=fulltx, whichtoinvert=c(F,F,T,F)) |
| save mask | antsImageWrite(maskM,paste(sub,taskout,"mask.nii.gz",sep="")) |
| pretty up for pdf | boldpreMa=getAverageOfTimeSeries(boldpreM) |
| boldpreMa=boldpreMa*boldpreMa |
| save nuisance variables | write.csv(boldpre$nuisanceVariables,paste("NV",sub,taskout,"1D",sep=".")) |
| make pdf | pdf(paste("Prep",sub,taskout,"pdf",sep=".")) |
| output orig space alignment | plot(t1,avgn4,window.overlay=c(mean(avgn4),Inf)) |
| output DVARS | plot(boldpre$DVARS,type='l',ylab='DVARS',xlab='TR') |
| output mni space alignment | plot(bmbrain,boldpreMa,axis=3) |
| save pdf | dev.off() |
| }} |
| clean up for loop | try(rm(list=c("bold","boldpre","boldstm","boldstmw","boldpreM"))) |
| }} |