Phantom Morphometry Study with ANTsR

Brian B. Avants et al.

2017-07-30

# set this for your own compilation
set.seed( 11 )
library( ANTsR )
## Loading required package: ANTsRCore
## 
## Attaching package: 'ANTsRCore'
## The following object is masked from 'package:stats':
## 
##     var
## The following objects are masked from 'package:base':
## 
##     all, any, apply, max, min, prod, range, sum
library( pander )
library( visreg )
bd="/Users/bavants/data/antsExamples/ANTsTutorial/"
templatefn=paste(bd,"phantomData/phantomtemplate.jpg",sep='')
if ( ! file.exists(templatefn) )
  stop( paste( "no image", templatefn ) )
allfns=Sys.glob( paste(bd,"phantomData/ph*wmgm.jpg",sep='')  )
demog=read.csv(paste(bd,'phantomData/phantpredictors.csv',sep=''))
demog$vol[ demog$vol == 1 ] = "Smaller"
demog$vol[ demog$vol == 0 ] = "Bigger"
ilist=imageFileNames2ImageList( allfns )
sccansig=0

Phantom Morphometry Study

Phantom Morphometry Study

Look at population: S1

plot( ilist[[1]], window.img=c(0,255) )

NULL

Look at population: S2

plot( ilist[[2]], window.img=c(0,255) )

NULL

Look at population: S5

plot( ilist[[5]], window.img=c(0,255) )

NULL

Look at population: S6

plot( ilist[[6]], window.img=c(0,255) )

NULL

Segmentation

First: segment the images

template=antsImageRead( templatefn )
ilist=imageFileNames2ImageList( allfns )
plist=imageFileNames2ImageList( allfns )
for ( x in 1:length(ilist) )
  {
  mask = getMask( ilist[[ x ]] )  %>% iMath("MD",1)
  img  = antsImageClone( ilist[[ x ]] )
  km   = kmeansSegmentation(  img, k = 3, kmask = mask ) # segment
  gm   = km$probabilityimages[[2]]
  plist[[ x ]] = gm
  }

Quick look at segmentation results

plot( ilist[[1]] )

## NULL

Quick look at segmentation results

mask = getMask( ilist[[1]] )  %>% iMath("MD",1)
img  = ilist[[ 1 ]]
km   = kmeansSegmentation(  img, 3, mask ) # segment
plot( km$segmentation, window.img=c(0,3) )

## NULL

Registration

Second: register the images

template = antsImageRead( templatefn )
jlist = imageFileNames2ImageList( allfns )
txList = list( )
jList2 = list( )
myvolsj = rep( NA, length( ilist ) )
myvolsi = rep( NA, length( ilist ) )
for ( x in 1:length(ilist) )
  {
  gm = plist[[x]]
  tx = antsRegistration( template, ilist[[x]], "SyNCC", 
                         gradStep = 0.1, synSampling = 2,
                         flowSigma = 3, totalSigma = 0.0,
                         regIterations = c(100,50,40) )
  gmw = antsApplyTransforms( template, gm, tx$fwdtransforms)
  plist[[ x ]] = gmw
  jlist[[ x ]] = createJacobianDeterminantImage( template,
    tx$fwdtransforms[1], 1, 0 ) %>% smoothImage(1.5)
  txList[[ x ]] = composeTransformsToField( template, tx$fwdtransforms )
  jList2[[ x ]] = createJacobianDeterminantImage(  template, txList[[ x ]] , FALSE, TRUE )
  lo = 120
  hi = 150 # cortex
 # lo = 150
#  hi = Inf # wm
  myvolsi[ x ] = sum( thresholdImage( ilist[[x]], lo, hi  ) )
  myvolsj[ x ] = sum( thresholdImage( template, lo, hi ) * jList2[[ x ]] )
  }
# Validate the jacobian calculations against the "real" volume change
# We must look at both the deformable and affine component, together.
plot( myvolsi, myvolsj , main='correlation in image-based and jacobian-based volume')

hist( myvolsi - myvolsj , main='difference in volume' )

print( paste("cor:", cor( myvolsi, myvolsj  ), "err:", mean( abs( myvolsi - myvolsj ) ), "serr:", mean( myvolsi - myvolsj ) ) )
## [1] "cor: 0.992319514493284 err: 54.6502276211977 serr: -53.3713649809361"
#  "cor: 0.99426596560049 err: 21.2326674573123 serr: 11.3190616928041" grad 0.05 synsampling = 1 flow 3 total 0

Quick look at the template …

plot( template )

## NULL

And a registration

plot( tx$warpedmovout  )

## NULL

And a registration overlay

plot( template, tx$warpedmovout %>% iMath("Canny",1,5,12),
  window.overlay=c(0.5,1)  )

## NULL

Jacobian

Now look at a log-jacobian image …

plot( jlist[[1]])

## NULL

another log-jacobian image …

plot( jlist[[5]])

## NULL

Now we can do statistics on either warped segmentations or the log-jacobian

…. but we should look at the demographics first.

Statistics

Always the first thing …

Look at the demographics file …

pander( demog )
vol gender
Bigger M
Bigger F
Bigger M
Bigger F
Smaller M
Smaller F
Smaller F
Smaller M

Let’s do a quick t-test on the jacobian

templatemask=getMask( template )
jmat=imageListToMatrix( jlist, templatemask )
# simple t-test to look at overall trends
rjmeans = rowMeans( jmat )
g1 = demog$vol == "Bigger"
g2 = demog$vol == "Smaller"
print( t.test( rjmeans[g1], rjmeans[g2] , paired=FALSE ) ) 
## 
##  Welch Two Sample t-test
## 
## data:  rjmeans[g1] and rjmeans[g2]
## t = 2.066, df = 5.702, p-value = 0.0868
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.009539579  0.105164016
## sample estimates:
##   mean of x   mean of y 
## -0.02601552 -0.07382774

Let’s do a quick voxel-wise test on the jacobian

Check out bigLMStats. It is a great little function for large (on the left) statistical studies …

demog$vol = factor( demog$vol )
umdl=lm( jmat ~ vol + gender, data=demog )
mdl=bigLMStats( umdl )
volstats=mdl$beta.pval[ "volSmaller", ]
qv=p.adjust(volstats,'BH')
qvimg=makeImage( templatemask, 1 - qv )

for large on the right studies, use a quick implementation in RcppEigen or RcppArmadillo

Look at results at one voxel

Use visreg to see the regression at the voxel with the maximum significant difference.

vdemog = data.frame( demog, voxel = jmat[,which.min(qv)])
vmdl=lm( voxel ~ vol + gender, data=vdemog )
visreg::visreg( vmdl, "vol", main="Regression at a voxel" )

Look at results in the image space

plot( template, qvimg, window.overlay=c(0.95,1))

## NULL

Exercise … repeat for the GM …

…

templateCortexMask = thresholdImage( template, lo, hi )
gmat=imageListToMatrix( plist, templateCortexMask  )
mdl=lm( gmat ~ vol + gender, data=demog )
mdl=bigLMStats( mdl )
gmstats=mdl$beta.pval[ "volSmaller", ]
gqv=p.adjust(gmstats,'BH')
gqvimg=makeImage( templatemask, 1-gqv )
## Warning in makeImage(templatemask, 1 - gqv): Number of voxels not the same
## as the positive values

Exercise … use a global covariate …

… this is always important to think about …

Multivariate statistics

Try a multivariate version of the same study

Do we believe that many aspects of cognition are fundamentally multivariate?

…. multivariate statistics in R

volClasses = c( rep( 1, 4 ), rep(0, 4 ) )
mypreds<-scale( as.matrix( cbind( volClasses, as.numeric( demog$gender ) ) ) )
scaledMatrix = scale( jmat )
sccan<-sparseDecom2(
  inmatrix  = list( scaledMatrix, mypreds ),
  inmask    = c( templatemask , NA ) ,
  mycoption = 0,
  sparseness=c( 0.25, 0.5 ),
  nvecs=3, its=5,
  cthresh=c(50,0),
  smooth=0.5, perms=200, robust=1 ) # might increase perms to be > 100
sccanimgs = matrixToImages( t(sccan$eig1) , templatemask )
sccansol <-  abs(sccanimgs[[1]]) %>% iMath("Normalize")
sccansol2 <- abs(sccanimgs[[2]]) %>% iMath("Normalize")

Significance based on permutation testing …

pander( sccan$ccasummary )
corrs pvalues
0.9954 0.005
0.9146 0.51
0.1518 0.6

Exercise: What happens to the permutation based significance when you vary parameters?

Visualizing multivariate statistics in R

plot( template, sccansol, window.overlay=c(0.1,1.01) )

## NULL
plot( template, sccansol2, window.overlay=c(0.1,1.01) )

## NULL

Another version of dual modality sparse regression

No explicit use of demographics here ….

spatmat = t( imageDomainToSpatialMatrix( templatemask, templatemask ) )
smoomat = knnSmoothingMatrix( spatmat, k = 27, sigma = 25 )
## Note: method with signature 'CsparseMatrix#Matrix#missing#replValue' chosen for function '[<-',
##  target signature 'dsCMatrix#lsCMatrix#missing#numeric'.
##  "Matrix#lsparseMatrix#missing#replValue" would also be valid
spatmat = t( imageDomainToSpatialMatrix( templateCortexMask, templateCortexMask ) )
smoomat2 = knnSmoothingMatrix( spatmat, k = 27, sigma = 100 )
# map these - via matrix - to observed features
params = matrix( nrow = 2, ncol = 3 )
params[1,] = c(1,2,1)
params[2,] = c(2,1,1)
x = list( scale( jmat, scale = F ), scale( gmat, scale = F ) )
jj = jointSmoothMatrixReconstruction( x, 2, params,
  gamma = 0.01, sparsenessQuantile=0.66, iterations=25,
  subIterations=20, positivity = TRUE,
  smoothingMatrix = list(smoomat2,smoomat), verbose=F )
p1 = jmat %*% jj$v[[2]]
p2 = gmat  %*% jj$v[[1]]
diag( cor( p1, p2 ) )
##        u.1        u.2 
## -0.5822389 -0.5025775
mdl1=lm( p1 ~ vol + gender, data=demog )
mdl2=lm( p2 ~ vol + gender, data=demog )
mdllist = lappend( summary( mdl1 ), summary( mdl2 ) )
pander( mdllist[[1]] )
  Estimate Std. Error t value Pr(>|t|)
volSmaller -0.0101 0.001946 -5.19 0.003496
genderM 0.0005721 0.001946 0.2939 0.7806
(Intercept) 0.004767 0.001686 2.828 0.03675
Fitting linear model: u.1 ~ vol + gender
Observations Residual Std. Error \(R^2\) Adjusted \(R^2\)
8 0.002753 0.8439 0.7814 
pander( mdllist[[2]] )
  Estimate Std. Error t value Pr(>|t|)
volSmaller -0.002941 0.0005752 -5.113 0.003728
genderM 0.0005423 0.0005752 0.9427 0.3891
(Intercept) 0.0002254 0.0004982 0.4525 0.6699
Fitting linear model: u.2 ~ vol + gender
Observations Residual Std. Error \(R^2\) Adjusted \(R^2\)
8 0.0008135 0.8439 0.7815 
pander( mdllist[[3]] )
  Estimate Std. Error t value Pr(>|t|)
volSmaller 0.0002413 0.0001212 1.992 0.103
genderM 0.0001996 0.0001212 1.648 0.1603
(Intercept) 0.003344 0.0001049 31.87 5.715e-07
Fitting linear model: u.1 ~ vol + gender
Observations Residual Std. Error \(R^2\) Adjusted \(R^2\)
8 0.0001713 0.572 0.4008 
pander( mdllist[[4]] )
  Estimate Std. Error t value Pr(>|t|)
volSmaller 5.381e-05 2.444e-05 2.202 0.07893
genderM 3.92e-05 2.444e-05 1.604 0.1697
(Intercept) 0.0008137 2.117e-05 38.44 2.245e-07
Fitting linear model: u.2 ~ vol + gender
Observations Residual Std. Error \(R^2\) Adjusted \(R^2\)
8 3.457e-05 0.5974 0.4364 
# look at the predictor image
temp = makeImage( templatemask, abs(jj$v[[2]][,1]) ) %>% iMath( "Normalize" )
plot( template, temp, window.overlay = c( 0.1, 1 ) )

## NULL
temp = makeImage( templateCortexMask, abs(jj$v[[1]][,2]) ) %>% iMath( "Normalize" )
plot( template, temp, window.overlay = c( 0.1, 1 ) )

## NULL

And MNI space ….

Mapping to MNI space with ANTsR

if ( ! exists("tem") )
  tem<-antsImageRead( getANTsRData("ch2") )

Mapping to MNI space with ANTsR

clust <- thresholdImage( tem, 82, 90 ) %>%
   iMath("ME",1)  %>% labelClusters()

Mapping to MNI space with ANTsR

if ( ! exists("mymni") ) {
# try getANTsRData if you have www access
  mymni<-list( antsImageRead(getANTsRData("mni") ),
               antsImageRead(getANTsRData("mnib") ),
               antsImageRead(getANTsRData("mnia") ) )
}

Mapping to MNI space with ANTsR

template_cluster_pair<-list(tem,clust)
t1=Sys.time()
gcoords<-getTemplateCoordinates( template_cluster_pair ,
    mymni , convertToTal = TRUE )
t2=Sys.time()
print(t2-t1)
## Time difference of 33.02606 secs

Reporting cluster coordinates in MNI space with ANTsR

Show the table of results.

pander( gcoords$templatepoints  )
x y z t Brodmann AAL
-12 13 -1 0 0 Caudate_L
13 15 6 0 0 Caudate_R
-8 -24 7 0 0 Thalamus_L
19 12 -7 0 48 Putamen_R
7 -23 6 0 0 Thalamus_R
28 2 -16 0 34 Amygdala_R
0 -58 -25 0 0 Vermis_8
40 20 -5 0 47 Insula_R
17 -80 -34 0 0 Cerebelum_Crus2_R

Review