Automatically label peak voxels from mincLM LMER with an atlas

This page builds on previous documentation from MICe.

What you need:

1. a mincLM or mincLMER object
2. a set of labels (such as from DSURQE) on your average. Pro tip: these can be obtained by the same methods as you register a mask onto your average, described here. Essentially, you're registering a label file onto your average, instead of a mask. Make sure your template and model have the same resolution.
3. a csv of definitions for the labels. The label values are numeric (such as 1), so you need the csv to link the numeric value to a region name (such as 1 is "left lateral ventricle")

Here we will a) automatically detect peak voxels b) connect them to labels c) produce plots for these peak voxels.

Automatically detecting peak voxels.

In an R script, first load in the csv and set of labels.

#definitions of labels
defs = read.csv("../../make_mask/DSURQE_40micron_R_mapping.csv")
#Labels.
atlasLabel <- mincArray(mincGetVolume("../../make_mask/labels_on_average.mnc"))

As in when you visualize your model results (see here), we must extract a matrix of the thresholds of statistical significance for the terms in our model:

thresholds = attr(mincFDR(model), "thresholds")

Then, from thresholds, you can extract the t-value corresponding with a specific FDR level. For example, here we select the t-value corresponding with 1% FDR or 0.01.

lowerthreshold = round(thresholds["0.01",'tvalue-sexmale'],digits=2)

Now, we can call the function to identifiy peak voxels. Here, inputStats points to our model, column identifies which term to use, direction ensures it will find positive and negative, threshold indicates only peaks more significant than 1% FDR will be identified, and minDistance specifies that peaks should not be identified between 1mm of each other (in order to ensure only one voxel from clusters of peak voxels are highlighted).

peaks <- mincFindPeaks(inputStats = mlm2,column = "tvalue-sexmale",direction="both", threshold = lowerthreshold, minDistance = 1)

The variable peaks is now a dataframe indicating the coordinates of the voxel and the value.

Labelling peak voxels

Now that we identified peak voxels, we need to connect them to a label. We will define a function that will do this.

The input to this function is the csv we called "defs"

label_peaks <- function(defs){
    #Create a variable in peaks called "label" that takes the coordinates from peaks (in columns 1:3) and extracts the label from the atlas corresponding with that voxel coordinate.
    peaks$label <- as.integer(atlasLabel[as.matrix(peaks[,1:3])])
    #Pipe the first three columns of defs to a function (gather) that will combine left.label and right.label into a single column called "variable" and the value of the label into a single column called value
    mdefs <- defs[,1:3] %>% gather("variable","value", c("left.label", "right.label"))
    #remove the words "label" from left and right
    mdefs$variable <- sub('.label','', mdefs$variable, fixed=T)
    #use paste to add left and right to the structure name. 
    mdefs$Structure <- paste(mdefs[,2], mdefs[,1])
    #In case any of your peak voxels are in an ulabelled region, add a 0 label for unlabelled regions
    mdefs <- rbind(mdefs, c("unlabelled","both",0))
    #turn value (the label) into a numeric
    mdefs$value <- as.numeric(mdefs$value)
    #go through peaks, add a column labelling the peak. 
    for (i in 1:nrow(peaks)){
      if(any(peaks$label[i]==mdefs$value)){
        ars<-which(peaks$label[i]==mdefs$value)
        peaks$label[i]<-mdefs$Structure[ars]
      }
     }
      return(peaks)
}

Now, you can call the function.

peaks <- label_peaks(defs)

The output, should be a dataframe with the label as one column:

> head(peaks)
  d1  d2 d3           x          y   z value                                           label
1 61  81 79 -0.26999989 -0.1899995 3.6 7.813                  left Cingulate cortex: area 30
2 62  91 78 -0.16999989  0.8100006 3.5 7.598                     left Secondary motor cortex
3 64 119 74  0.03000011  3.6100006 3.1 7.308                  left Cingulate cortex: area 32
4 66 103 76  0.23000012  2.0100006 3.3 6.400                right Cingulate cortex: area 24b
5 71  77 80  0.73000012 -0.5899995 3.7 6.350                                      unlabelled
6 52  61 77 -1.16999990 -2.1899995 3.4 6.236 left Secondary visual cortex: mediolateral area