ANTsR Notes: 11/8/2019
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.
To allow for flexibility in processing we have now split scripts into
- Preprocessing: R_preProcessMR290.R (txt)
- Postprocessing: R_fmri_robust_MRttl.R (txt)
Preprocessing walkthrough:
| 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"))) |
| }} |