In Mouse-Imaging-Centre/RMINC: Statistical Tools for Medical Imaging NetCDF (MINC) Files

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)
library(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)
mask <- "/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)

And we take a look

# first step - load the background anatomy
anatVol <- mincGetVolume("/Users/jason/data/CREB/21aug15_est_conserv-nlin-3.mnc")
# 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)

Sure 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)

Mouse-Imaging-Centre/RMINC documentation built on Nov. 12, 2022, 1:50 p.m.