RMINC visualization
Jason Lerch
February 20, 2016
Here we will step through some of the new visualization tools built into RMINC 1.3. So, let’s start with getting a dataset ready - we’ll look at sex differences in the brain.
# load libraries
library(RMINC)## Warning in rgl.init(initValue, onlyNULL): RGL: unable to open X11 display## Warning: 'rgl_init' failed, running with rgl.useNULL = TRUElibrary(plotrix) # for legends
# some weirdness potentially particular to my mac. Images are flipped, so fix that:
options(RMINC_flip_image=TRUE)
# load an existing dataset
load("~/data/CREB/CREB-analyses/wtonly.RData")
table(gfwt$Sex)##
## F M
## 38 42mask <- "/Users/jason/data/CREB/mask.mnc"Let’s start by running a linear model testing the effect of sex on local brain volume.
vs <- mincLm(relJac02 ~ Sex, gfwt, mask=mask)## Method: lm
## Number of volumes: 80
## Volume sizes: 201 439 298
## N: 80 P: 2
## In slice
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200
## DoneAnd we take a look
# first step - load the background anatomy
anatVol <- mincGetVolume("/Users/jason/data/CREB/21aug15_est_conserv-nlin-3.mnc")## Start: 0 0 0
## Count: 201 439 298# and a series of slices through the brain
mincPlotSliceSeries(mincArray(anatVol), # the anatomical volume
mincArray(vs, "tvalue-SexM"), # pull out one column of the stats
anatLow=700, anatHigh=1400, # set anatomy thresholds
low=2.5, high=10, # set stats thresholds
symmetric=T) # show separate upper and lower
This introduces the first, and one of the most useful, of the visualization functions: mincPlotSliceSeries. It does what the name claims - shows a series of slices through the brain, in this case overlaying the mincLm results on the average anatomy background. Before proceeding, one note: mincArray is called for both the anatomical volume and the stats; for the moment it’s a shortcut to add dimension information to each volume. Future versions of RMINC will likely get rid of that and store the necessary information right after outputs are created by the likes of mincLm.
Anyway, let’s make the plot a bit prettier: there are too many slices at the beginning and end that show no useful info, so let’s get rid of them.
mincPlotSliceSeries(mincArray(anatVol), # the anatomical volume
mincArray(vs, "tvalue-SexM"), # pull out one column of the stats
anatLow=700, anatHigh=1400, # set anatomy thresholds
low=2.5, high=10, # set stats thresholds
symmetric=T, # show separate upper and lower
begin=50, end=-90) # remove slices from both sides 
Looks prettier already. But it would help to have an indication of where the slices are and what the colours represent. So let’s add that:
mincPlotSliceSeries(mincArray(anatVol), # the anatomical volume
mincArray(vs, "tvalue-SexM"), # pull out one column of the stats
anatLow=700, anatHigh=1400, # set anatomy thresholds
low=2.5, high=10, # set stats thresholds
symmetric=T, # show separate upper and lower
begin=50, end=-90, # remove slices from both sides
legend="t-statistics") 
And let’s show the other dimensions:
mincPlotSliceSeries(mincArray(anatVol), # the anatomical volume
mincArray(vs, "tvalue-SexM"), # pull out one column of the stats
anatLow=700, anatHigh=1400, # set anatomy thresholds
low=2.5, high=10, # set stats thresholds
symmetric=T, # show separate upper and lower
begin=35, end=-35 , # remove slices from both sides
legend="t-statistics",
dimension = 1) 
mincPlotSliceSeries(mincArray(anatVol), # the anatomical volume
mincArray(vs, "tvalue-SexM"), # pull out one column of the stats
anatLow=700, anatHigh=1400, # set anatomy thresholds
low=2.5, high=10, # set stats thresholds
symmetric=T, # show separate upper and lower
begin=25, end=-25 , # remove slices from both sides
legend="t-statistics",
dimension = 3) 
Let’s do something fancier - or at least requiring more customization. Let’s create a slice showing the estimated effect (i.e. the beta) with the significant regions contoured.
The core function underlying it all is mincImage. So we begin by putting together the plot with a simple rendition of the underlying anatomy.
mincImage(mincArray(anatVol), slice=250)
A good start, but the axes are quite unnecessary, so let’s get rid of them.
mincImage(mincArray(anatVol), slice=250, axes=F)
And let’s return to a cleaner setting of thresholds:
mincImage(mincArray(anatVol), slice=250, axes=F, low=700, high=1400)
Now we add the beta coefficients on top.
mincImage(mincArray(anatVol), slice=250, axes=F, low=700, high=1400)
mincImage(mincArray(vs, "beta-SexM"), slice=250, low=0.02, high=0.1, add=T, col=heat.colors(255), underTransparent = T)
Starting to look good, but we’ll want slighty different colours and both the positive and the negative effects shown:
poscolours = colorRampPalette(c("red", "yellow"))(255)
negcolours = colorRampPalette(c("blue", "turquoise1"))(255)
mincImage(mincArray(anatVol), slice=250, axes=F, low=700, high=1400)
mincImage(mincArray(vs, "beta-SexM"), slice=250, low=0.02, high=0.1, add=T, col=poscolours, underTransparent = T)
mincImage(mincArray(vs, "beta-SexM"), slice=250, low=-0.02, high=-0.1, add=T, col=negcolours, underTransparent = T)
This slice shows any effects above (or below) a log jacobians of 0.02. Is any of that significant? Let’s find out via FDR.
mincFDR(vs, mask=mask)##
## Computing FDR threshold for all columns
## Start: 0 0 0
## Count: 201 439 298
## Computing threshold for F-statistic
## Computing threshold for tvalue-(Intercept)
## Computing threshold for tvalue-SexM## Multidimensional MINC volume
## Columns: F-statistic tvalue-(Intercept) tvalue-SexM
## NULL
## Degrees of Freedom: c(1, 78) 78 78
## FDR Thresholds:
## F-statistic tvalue-(Intercept) tvalue-SexM
## 0.01 12.720632 2.830987 3.569582
## 0.05 7.627336 2.159808 2.763280
## 0.1 5.579595 1.821244 2.363340
## 0.15 4.411686 1.601513 2.101459
## 0.2 3.597821 1.432507 1.897762Sure enough, there is. So let’s add contours to where we have significance.
poscolours = colorRampPalette(c("red", "yellow"))(255)
negcolours = colorRampPalette(c("blue", "turquoise1"))(255)
mincImage(mincArray(anatVol), slice=250, axes=F, low=700, high=1400)
mincImage(mincArray(vs, "beta-SexM"), slice=250, low=0.02, high=0.1, add=T, col=poscolours, underTransparent = T)
mincImage(mincArray(vs, "beta-SexM"), slice=250, low=-0.02, high=-0.1, add=T, col=negcolours, underTransparent = T)
mincContour(abs(mincArray(vs, "tvalue-SexM")), slice=250, levels=c(2.36, 3.57), add=T)
We’ve added lines for 10% and 1%; let’s label them and make them easier to distinguish:
poscolours = colorRampPalette(c("red", "yellow"))(255)
negcolours = colorRampPalette(c("blue", "turquoise1"))(255)
mincImage(mincArray(anatVol), slice=250, axes=F, low=700, high=1400)
mincImage(mincArray(vs, "beta-SexM"), slice=250, low=0.02, high=0.1, add=T, col=poscolours, underTransparent = T)
mincImage(mincArray(vs, "beta-SexM"), slice=250, low=-0.02, high=-0.1, add=T, col=negcolours, underTransparent = T)
mincContour(abs(mincArray(vs, "tvalue-SexM")), slice=250, levels=c(2.36, 3.57),
labels=c("10%", "1%"), lty=c(2,1), add=T)
Fairly pretty, but I’d rather have the labels in a legend rather than on the image itself. Doable:
poscolours = colorRampPalette(c("red", "yellow"))(255)
negcolours = colorRampPalette(c("blue", "turquoise1"))(255)
mincImage(mincArray(anatVol), slice=250, axes=F, low=700, high=1400)
mincImage(mincArray(vs, "beta-SexM"), slice=250, low=0.02, high=0.1, add=T, col=poscolours, underTransparent = T)
mincImage(mincArray(vs, "beta-SexM"), slice=250, low=-0.02, high=-0.1, add=T, col=negcolours, underTransparent = T)
mincContour(abs(mincArray(vs, "tvalue-SexM")), slice=250, levels=c(2.36, 3.57),
drawlabels=F, lty=c(2,1), col="white", add=T)
legend(200, 200, c("10% FDR", "1% FDR"), lty=c(2,1), col="white",
text.col="white", bg="black", bty="n")
Still missing the colour bars - let’s add those to the bottom. This will also need some mucking about with the par plot settings to get the whole background to be black and pretty.
poscolours = colorRampPalette(c("red", "yellow"))(255)
negcolours = colorRampPalette(c("blue", "turquoise1"))(255)
graycolours = gray.colors(255)
opar <- par(bg=graycolours[1]) # set the background to be the same colour as the under colour
mincImage(mincArray(anatVol), slice=250, axes=F, low=700, high=1400, col=graycolours)
mincImage(mincArray(vs, "beta-SexM"), slice=250, low=0.02, high=0.1, add=T, col=poscolours, underTransparent = T)
mincImage(mincArray(vs, "beta-SexM"), slice=250, low=-0.02, high=-0.1, add=T, col=negcolours, underTransparent = T)
mincContour(abs(mincArray(vs, "tvalue-SexM")), slice=250, levels=c(2.36, 3.57),
drawlabels=F, lty=c(2,1), col="white", add=T)
legend(200, 200, c("10% FDR", "1% FDR"), lty=c(2,1), col="white",
text.col="white", bg="black", bty="n")
color.legend(20, -20, 130, 0, c("-0.1", "-0.02"), rev(negcolours), col="white", align="rb")
color.legend(170, -20, 280, 0, c("0.02", "0.1"), poscolours, col="white", align="rb")
# restore the old settings
par(opar)