Running arbitrary voxelwise functions in parallel (qMincApply)

Defining and running arbitrary functions across all voxels in parallel

There may come a time when you want to run a function across all the voxels in your data and there is no built-in function in RMINC to accomplish this. Some functions, such as mean and standard deviation are already built-in statistics that can be applied (see here: https://mouse-imaging-centre.github.io/RMINC/reference/mincSummary.html), but what if you want to calculate the mean, add one, and then raise it to the power of 15 or something silly like that?

The RMINC function qMincApply (https://mouse-imaging-centre.github.io/RMINC/reference/qMincApply.html) is the perfect function for you! It allows you to pass in a list of filenames (which point to minc files), perform a function on each voxel in parallel, and output the results. This output can be used in R or saved as a minc file.

For this function, you will need:

library(RMINC)
library(batchtools)

Set up your csv as you would for running voxelwise linear models or linear mixed effects models in R with RMINC. Specifically you need, at the very least a csv with one column that includes paths to each minc file you want to use. For this example, the dataframe is called data and the column with the mincs is called file. You'll also need to define a variable that points to a mask in minc format.

# load in data
data <- read.csv("mncl_full_smooth_jacs.csv")
averagemask = "manual_mask.mnc"

The next thing you have to do is define your function. the function below, called my_fun takes as input a list of voxels, calculates the mean, raises it to the 15th power, then returns the output.

my_fun = function(vox_dat){
    mew = mean(vox_dat)
    out = (mew+1)^15
    return(out)
}

After defining the function, you need to call it in qMincApply. Here we pass qMincApply several important variables: data$file : this is the data frame specifying the column that has the paths to the minc files. fun = my_fun : This points to the function we previously defines mask = averagemask : This specifies the mask within which we would like to include the voxels in our analysis. wait = FALSE : This is an important thing to include when running the function on Trillium/scinet. If you exclude it, jobs will be submitted to the cluster, then qMincApply will try to reduce the outputs of these jobs to a collated object in R, even if the jobs haven't finished running yet! by setting wait = FALSE, we allow the jobs to be queued and run. By default, this command will create 4 jobs of 15 minutes each. For bigger jobs, see the section below on Timing and Scheduling qMincApply Jobs.

qMincApply(data$file, fun = my_fun, mask = averagemask, wait = FALSE)

If you are on Trillium/scinet, after you submit the jobs with qMincApply, you can check their status in a terminal with sq. Once the jobs have finished, run the following code to automatically reduce and collate the output into an object:

qminc_test <- qMincReduce(batchtools::getDefaultRegistry(), collate = simplify2minc)

This code creates an object called qminc_test. If you look at its structure, you should see something like the following:

str(qminc_test)
 'mincSingleDim' Named num [1:5429424] 0 0 0 0 0 0 0 0 0 0 ...
 - attr(*, "names")= chr [1:5429424] "inds" "vals" NA NA ...
 - attr(*, "filenames")= chr [1:165] "resample_dbm_full_smooth_mncl/MCH_NEO_THC_01_01_10_1-2-1_fix_lsq6_manual_extracted_fwhm_4vox.mnc" "resample_dbm_full_smooth_mncl/MCH_NEO_THC_01_01_03_1-2-1_fix_lsq6_manual_extracted_fwhm_4vox.mnc" "resample_dbm_full_smooth_mncl/MCH_NEO_THC_01_01_05_1-2-1_fix_lsq6_manual_extracted_fwhm_4vox.mnc" "resample_dbm_full_smooth_mncl/MCH_NEO_THC_01_01_07_1-2-1_fix_lsq6_manual_extracted_fwhm_4vox.mnc" ...
 - attr(*, "likeVolume")= chr "resample_dbm_full_smooth_mncl/MCH_NEO_THC_01_01_10_1-2-1_fix_lsq6_manual_extracted_fwhm_4vox.mnc"
 - attr(*, "mask")= chr "/scratch/lanicupo/neonate_power/mncl2_investigation/voxelwise/data/pMincMask429ee7588d27c.mnc"

It is a minc object and can be written out into a minc like this:

mincWriteVolume(qminc_test, "my_fun_minc.mnc", like.filename = averagemask)

Now, we can visualize it in tools like Display or register from the terminal:

Timing and Scheduling qMincApply Jobs on Trillium