R function to write CIFTI files

Update: This post is outdated.  Use our ciftiTools R package instead!  ciftiTools is available to install from CRAN in R:

install.packages(“ciftiTools”)

or from GitHub for the latest release:

library(devtools)

install_github(‘mandymejia/ciftiTools’)

For more information about ciftiTools, check out our paper in NeuroImage:

Pham, D., Muschelli, J., & Mejia, A. (2022). ciftiTools: A package for reading, writing, visualizing, and manipulating CIFTI files in R. NeuroImage, 118877. https://doi.org/10.1016/j.neuroimage.2022.118877

As I wrote about previously, CIFTI files are a new and useful file format for storing, analyzing and visualizing imaging data on the cortical surface AND subcortical areas (gray matter structures and cerebellum).  Much of the imaging data available through the Human Connectome Project is in CIFTI format, so anyone who wants to be able to work with that data needs to be able to work with CIFTI files.  Some tools have been developed to do this, but they are not yet available for R (to my knowledge at the time of writing).  My solution until recently was to just use MATLAB for all of my analyses, but I’m working on a new project that uses INLA, which has been programmed in R.  Cue the tedious back-and-forth between R and MATLAB.  Today I finally had enough, so I wrote up a quick function to write CIFTI file(s) from a vector or matrix of values.  (For an intro on how to use MATLAB from within R, see this previous post.)  For visualization purposes, the Connectome Workbench is still needed, but this has drastically streamlined my workflow.  Later I may wrap this into a package, but for now here’s the function code.

You’ll need to replace the default values for the template and toolboxloc arguments, but I leave them in as an example.  The CIFTI file you point to should be of the lower resolution variety.  As mentioned in the function comments, this assumes you have MATLAB installed and can run it from the command line (I have “module load matlab” in my .bashrc file).  You’ll also need to download the cifti-matlab or fieldtrip toolbox, both of which include the ft_read_cifti and ft_write_cifti functions.  I use cifti-matlab because fieldtrip has a bunch of other stuff I don’t need.  You also need a fair amount of memory to read and write CIFTI files, so I usually do this on the cluster.  My laptop has 8GB of RAM, which is not enough, but my desktop has 32G, and this seems to work.  Once I have the CIFTI files written, I just transfer them to my machine to visualize them with the Connectome Workbench.

Comments and suggestions welcome!

makeCIFTIs <- function(table, names, hemisphere='both', 
 template = '~/HCP/data/groupICA/groupICA_3T_Q1-Q6related468_MSMsulc_d25.ica/melodic_IC_ftb.dlabel.nii', 
 toolboxloc = '~/matlab_toolboxes/cifti-matlab/',
 scratchloc = '~/'){

 # table = VxQ matrix, each column will be written as a CIFTI file
 # table can be a vector if only a single CIFTI is needed
 # V = number of voxels in both or one hemisphere
 # Q = number of CIFTI files to be written
 # names = vector length Q, containing file names to be written (no extension)
 # hemisphere = 'both', 'left', 'right', or 'all' (includes subcortical voxels)
 # if hemisphere='both', the first half of table should be the left hemisphere
 # template = name of existing CIFTI file that can be used as a template
 # toolboxloc = location and name of folder containing cifti-matlab or fieldtrip toolbox
 # scratchloc = location where CSV files will be written then deleted

 # this function writes a MATLAB script and executes it
 # must have MATLAB installed
 # must have the cifti-matlab or fieldtrip toolbox installed
 # currently can only write CIFTIs with a vector of values (e.g. not a timeseries or a set of label maps)

if(is.vector(table)) table <- matrix(table, nrow=length(table, ncol=1))
 Q <- ncol(table)
 V <- nrow(table)

 if(hemisphere=='all') if(V != 96854) error('Number of rows must equal 96854 with hemisphere="all"')
 if(hemisphere=='both') if(V != 64984) error('Number of rows must equal 64984 with hemisphere="both"')
 if(hemisphere=='left') if(V != 32492) error('Number of rows must equal 32492 with hemisphere="left"')
 if(hemisphere=='right') if(V != 32492) error('Number of rows must equal 32492 with hemisphere="right"')

# Write table to CSV in scratch location
 tmp <- table
 tmp[is.na(tmp)] <- 0 #set NAs to zero
 fname <- file.path(scratchloc, 'tmp.csv')
 write.table(tmp, file=fname, row.names=FALSE, col.names=FALSE, sep=',')

 # Write MATLAB script

 line1 <- paste0("addpath '", toolboxloc, "'")
 line2 <- paste0("cifti = ft_read_cifti('",template,"');") #reads in structure with cell 'indexmax'
 line3 <- paste0("indsL = (cifti.brainstructure == 1);", #L cortex only
 "indsR = (cifti.brainstructure == 2);", #R cortex only
 "indsLR = (cifti.brainstructure <= 2);", #L and R cortices only
 "inds_nan = isnan(cifti.indexmax);") #location of NaNs

 line5 <- paste0("cd '", getwd(), "'")
 line6 <- paste0("fname = '",fname,"';",
 "coefs = csvread(fname);")

 line7 <- paste0("names = {'",paste(names_beta, collapse="';'"),"'};")

 line8 <- paste0("for k=1:",Q)
 line9 <- "k"
 line10 <- "vals = coefs(:,k);"
 line11 <- "cifti.indexmax = cifti.indexmax*0;"

 if(hemisphere=='all') line12 <- "cifti.indexmax = vals;"
 if(hemisphere=='left') line12 <- "cifti.indexmax(indsL) = vals;"
 if(hemisphere=='right') line12 <- "cifti.indexmax(indsR) = vals;"
 if(hemisphere=='both') line12 <- "cifti.indexmax(indsLR) = vals;"

 line13 <- "cifti.indexmax(inds_nan) = 0/0;" #put in NaNs
 line14 <- "ft_write_cifti(names{k}, cifti, 'parameter', 'indexmax');"
 line15 <- "end"

 matlab_lines <- c(line1, line2, line3, line5, line6, line7, line8,
 line9, line10, line11, line12, line13, line14, line15)
 writeLines(matlab_lines, con=file.path(scratchloc, 'myscript.m'))

 system("matlab -nodisplay -r \"run('~/myscript.m'); exit\"")

 file.remove(fname)
 file.remove(file.path(scratchloc, 'myscript.m'))
}

And since a pretty picture always helps, here’s one of the products of my current analysis based on a CIFTI file I wrote this way, visualized using the Connectome Workbench.