Example neurodegeneration study with tensor based morphometry

Brian B. Avants

01 August, 2017

example human brain study of neurodegeneration

Preprocessing.

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
rdir = "/Users/bavants/Box Sync/Mice_FAS_study/Workshop_Related/"
fns = Sys.glob( paste( rdir, "data/neurodegen/ndgen*png", sep='' ) )
layout( matrix( 1:length(fns), ncol=2, byrow=T) )
layout( matrix( 1:length(fns), nrow=2, byrow=F) )
template = antsImageRead( paste( rdir,"data/neurodegen/template2.png", sep='' ) )
templateMask = thresholdImage( template, 1, Inf ) %>% iMath( "FillHoles" )
templateMask = thresholdImage( kmeansSegmentation( template, 3 )$segmentation, 2, 2 )
rlist = list( )
jlist = list( )
plist = list( )
for ( i in 1:length( fns ) ) {
  tar = antsImageRead( fns[ i ] )
  tarseg = kmeansSegmentation( tar, 3 )
  plot( tar )
  rlist[[ i ]] = antsRegistration(  template, tar, typeofTransform = "SyNCC", totalSigma = 0.0 )
  plot(  rlist[[ i ]]$warpedmovout )
  jlist[[ i ]] = createJacobianDeterminantImage( template,
    rlist[[ i ]]$fwdtransforms[1], 1, 1 )
  plist[[ i ]] = antsApplyTransforms( template, tarseg$probabilityimages[[2]], 
                                      transformlist = rlist[[ i ]]$fwdtransforms )
  }

Voxel-based mass univariate statistics

layout( matrix( 1:6, nrow=2, byrow=F) )
jmat = imageListToMatrix( plist, templateMask )
ndemog = read.csv( paste( rdir,"data/neurodegen/demog.csv", sep='' ) )
ndemog$dx = factor( ndemog$dx )
umdl=lm( jmat ~ . , data=ndemog )
mdl=bigLMStats( umdl )
volstats=mdl$beta.pval[ "dx1", ]
qv=p.adjust(volstats,'BH')
print( range( qv ) )
## [1] 0.07762248 0.99971200
ndemog$globalMean = rowMeans( jmat )
umdl=lm( jmat ~ . , data=ndemog )
mdl=bigLMStats( umdl )
volstats=mdl$beta.pval[ "dx1", ]
qv=p.adjust(volstats,'BH')
qvimg=makeImage( templateMask, 1 - qv )
print( range( qv ) )
## [1] 0.009067891 0.999881109
timg=makeImage( templateMask, mdl$beta.t[ "dx1", ] )
timgMap = timg * thresholdImage( timg, 2, Inf )
plot( template, timgMap, window.overlay = range( timgMap )  )
## NULL
timgMap = ( timg * (-1) ) * thresholdImage( timg * (-1), 2, Inf )
plot( template, timgMap, window.overlay = range( timgMap )  )
## NULL
spatmat = t( imageDomainToSpatialMatrix( templateMask, templateMask ) )
smoomat = knnSmoothingMatrix( spatmat, k = 100, sigma = 10 )
## 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
sjmat = as.matrix( jmat %*% smoomat )
plot( makeImage( templateMask, jmat[1,] ) )
## NULL
plot( makeImage( templateMask, sjmat[1,] ) )
## NULL
umdl=lm( sjmat ~ . , data=ndemog )
mdl=bigLMStats( umdl )
volstats=mdl$beta.pval[ "dx1", ]
qv=p.adjust(volstats,'BH')
timg=makeImage( templateMask, mdl$beta.t[ "dx1", ] )
timgMap =  thresholdImage( timg * (1), 2, 100 ) * timg
plot( template, timgMap, window.overlay = c( 0, max(timgMap) )  )
## NULL
timgMap = ( timg * (-1) ) * thresholdImage( timg * (-1), 2, 100 )
plot( template, timgMap, window.overlay = c( 0, max(timgMap) )  )

## NULL

Eigenanatomy statistics

spatmat = t( imageDomainToSpatialMatrix( templateMask, templateMask ) )
smoomat = knnSmoothingMatrix( spatmat, k = 100, sigma = 2 )
rjmat = as.matrix( residuals( lm( jmat ~ globalMean , data=ndemog) ) %*% smoomat )
nc = 10
eanatStruct <- sparseDecom( 
  inmatrix   = rjmat,
  inmask     = templateMask,
  sparseness = 0.01, # set this explicitly, negative
  # values allow signed pseudo-eigenvectors in results
  nvecs      = nc,
  smooth     = 0.5,  # components will be smooth
  # negative values use edge-preserving methods (WIP)
  cthresh    = 50,  # get large components
  its        = 3,   # optimize for a "long" time
  mycoption  = 0)   # this controls orthogonality constraints
eimgs = matrixToImages( eanatStruct$eigenanatomyimages, templateMask )
eseg = eigSeg( templateMask, eimgs, FALSE )
plot( template, eseg, window.overlay = c( 0, max(eseg) )  )

## NULL
proj = jmat %*% t( eanatStruct$eigenanatomyimages )
umdl=lm( proj ~ . , data=ndemog )
mdl=bigLMStats( umdl )
volstats=mdl$beta.pval[ "dx1", ]
print( min( volstats ) )
## [1] 0.003745953
wm = which.min( volstats )
plot( template, eimgs[[wm]], window.overlay = c( 0, max( eimgs[[wm]]) )  )

## NULL
smoomat = knnSmoothingMatrix( spatmat, k = 27, sigma = 20 )
mysvd <- fastICA::fastICA( t( jmat ), n.comp = nc, verbose = F )
tdf = data.frame( mysvd$K )
tdf = data.frame(  jmat %*% t( icawhiten( jmat, nc )) )
# tdf = data.frame(  irlba::partial_eigen( jmat )$vec )
names( tdf ) = colnames( tdf ) = paste("u",1:ncol(tdf),sep='')
myform =  as.formula( paste( "x~",paste( colnames( tdf ), collapse='+' ) ))
x = scale( jmat )
smoo = smoothMatrixPrediction( x=x,
    basisDf = tdf, myform,  positivity = T,
    verbose = T, gamma = 0.001, iterations = 50, 
    sparsenessQuantile = 0.9,
    smoothingMatrix = smoomat  )
## [1] "iteration 0 err 0.717904956268178"
## [1] "2 0.721607551854146"
## reducing gradient step: 5e-04
## [1] "3 0.722508914650693"
## [1] "4 0.721920404353719"
## reducing gradient step: 0.00025
## [1] "5 0.722040891134609"
## [1] "6 0.721665634271011"
## [1] "7 0.721661188137733"
## reducing gradient step: 0.000125
## [1] "8 0.721723672107896"
## [1] "9 0.721545206518841"
## [1] "10 0.721521955599531"
## reducing gradient step: 6.25e-05
## [1] "11 0.72152637177949"
## [1] "12 0.721434967077574"
## [1] "13 0.721402541181637"
## [1] "14 0.721383954169074"
## [1] "15 0.721373118557886"
## [1] "16 0.721368909569382"
## [1] "17 0.721367551242716"
## reducing gradient step: 3.125e-05
## [1] "18 0.721369824445415"
## [1] "19 0.721324395097761"
## [1] "20 0.721303687076019"
## [1] "21 0.72128742461419"
## [1] "22 0.721276520748143"
## [1] "23 0.721267782838082"
## [1] "24 0.721260584347475"
## [1] "25 0.721253781702028"
## [1] "26 0.721248844884618"
## [1] "27 0.721245424550287"
## [1] "28 0.721242334662672"
## [1] "29 0.721240180644346"
## [1] "30 0.721238254296973"
## [1] "31 0.721236521004685"
## [1] "32 0.721236121614615"
## [1] "33 0.721236016856406"
## reducing gradient step: 1.5625e-05
## [1] "34 0.721236210939405"
## [1] "35 0.721213772518164"
## [1] "36 0.72120122405646"
## [1] "37 0.721190027715526"
## [1] "38 0.721180665360363"
## [1] "39 0.721172218823165"
## [1] "40 0.721164452217383"
## [1] "41 0.721157589752523"
## [1] "42 0.721150849919199"
## [1] "43 0.721144879132705"
## [1] "44 0.721139394912955"
## [1] "45 0.721134182975171"
## [1] "46 0.721129087120465"
## [1] "47 0.721124247543977"
## [1] "48 0.721119728782361"
## [1] "49 0.721115404504871"
## [1] "50 0.721111144632612"
## [1] "51 0.721106857585981"
## [1] "end 0.721106857585981"
nms =  colnames( smoo$v )
proj = jmat %*% ( smoo$v  )
umdl=lm( proj ~ . , data=ndemog )
mdl=bigLMStats( umdl )
volstats=mdl$beta.pval[ "dx1", ]
print( min( volstats ) )
## [1] 0.001788297
wm = which.min( volstats )
temp = makeImage( templateMask, abs( smoo$v[,wm] ) ) %>% iMath( "Normalize" )
plot( template, temp, window.overlay = c( 0.1, max( temp ) )  )

## NULL
# temp = makeImage( templateMask, abs( mysvd$S[,wm] ) ) %>% iMath( "Normalize" )
# plot( template, temp, window.overlay = c( 0.1, max( temp ) )  )