afni: 3dTstat with not-censored timepoints only

In a few recent posts I've shown images of the mean and standard deviation (calculated across time for each voxel), for QC tests. These are easy to calculate in afni (example here), but the 3dTstat command I used includes all timepoints (TRs), unless you specify otherwise. As described previously, we've been using a threshold of FD > 0.9 for censoring high-motion frames before doing GLMs. Thus, I wanted to calculate the mean and standard deviation images only including frames that were not marked for censoring (i.e., restrict the frames used by 3dTstat). This was a bit of a headache to code up, so R and afni code are after the jump, in the hopes it will be useful for others.


I mostly use R these days, and sometimes use R to generate the afni commands and execute them on the fly, as in this example. I ran this on a linux machine (with neurodebian); the system2 call in particular may need changing for other OS.

I include the full path to the afni commands (e.g., /usr/local/pkg/linux_openmp_64/1d_tool.py) to avoid path issues; replace as needed for your installation. This code takes two input files: the 4d NIfTI (_functional.nii.gz) and _FD.txt, a single column text file with the FD for each frame in the rows (so the number of rows in _FD.txt is the same as the 4th dimension in the NIfTI). The code then thresholds at FD > 0.9 and writes out a temporary text file (temp.1D, also one column, one row per frame), with 1s for volumes to keep and 0s for frames to censor (not include in the 3dTstat calculation). Of course, you could generate this censor list using whatever metrics or code you wish.

 in.path <- "/scratch2/input/";  # path to input files  
 out.path <- "/scratch2/output/";  # where to write output files  
 sub.ids <- c("subject1", "subject2"); # subject IDs, as in input file names  
   
 for (sid in 1:length(sub.ids)) {   # sid <- 1; rid <- 1;  
  fname.fd <- paste0(in.path, sub.ids[sid], "_FD.txt"); # 1D text file of framewise FD values (or whatever for censoring)  
  fname.4d <- paste0(in.path, sub.ids[sid], "_functional.nii.gz"); # 3d+time brain image; average each voxel over time  
  fname.out <- paste0(out.path, sub.ids[sid], "_sdcens.nii.gz"); # output filename (will be written by afni)  
  if (file.exists(fname.fd) & file.exists(fname.4d) & !file.exists(fname.out)) {  # check required files present (or not)  
   fd.vec <- read.table(fname.fd)[,1]; # read the FD column in as a vector  
     
   # make a single-column text file with 0 for volumes to OMIT (censor), 1 for volumes to KEEP.  
   # this makes a temporary 1D file from the FD values; you can make the file however you need, as long as it's in this format  
   out.vec <- rep(1, length(fd.vec));   
   out.vec[which(fd.vec > 0.9)] <- 0;   # set volumes above threshold to 0  
   write.table(out.vec, paste0(out.path, "temp.1D"), row.names=FALSE, col.names=FALSE); # write out temporary file for 1d_tool.py  
     
   # these next two commands actually do the work. The first gets the censored volumes (TRs with 0 in temp.1D) in tmp,   
   # in the afni-required text format (e.g., [5..8,15]).  
   tmp <- system2("/usr/local/pkg/linux_openmp_64/1d_tool.py", args=paste0("-infile ", out.path, "temp.1D -show_trs_uncensored encoded"), stdout=TRUE);  
     
   # now, 3dTstat is called, feeding it the frames in the tmp string.  
   system2("/usr/local/pkg/linux_openmp_64/3dTstat", args=paste0("-stdev -prefix ", fname.out, " ", fname.4d, "'[", tmp, "]'"), stdout=TRUE);  
  }  
 }  

The key part is calling 1d_tool.py to generate the list of volumes in the afni-required format (I put it into the tmp variable), then feeding this string back to afni in the 3dTstat command. These strings can be VERY long, which is why I found it easier to build the commands in R rather than directly.

Here are the key afni calls, without the R code:

 1d_tool.py -infile /scratch2/output/temp.1D -show_trs_uncensored encoded  
   
 3dTstat -stdev -prefix /scratch2/output/subject1_sdcens.nii.gz /scratch2/input/subject1_functional.nii.gz'[0..26,28..138,140..148]'  

Here, I copy-pasted the string returned by 1d_tool.py into the '[]' required by 3dTstat.