library(rmarkdown)
library(knitr)
opts_chunk$set(fig.width=9, fig.height=6, fig.path='Figs/',
               warning=FALSE, message=FALSE)
zscore<-TRUE

Introduction

Decoding from blood oxygenation-level dependent (BOLD) magnetic resonance imaging (MRI) data remains an esoteric art. Relatively few public datasets exist and these are rarely accompanied by accessible and validated, general-purpose decoding methodology. This package seeks to document some relatively standard approaches to processing BOLD data and follow with a few novel extensions. These are implemented via R, its packages and ANTsR.

RKRNS Algorithms

This package integrates several frameworks for BOLD processing:

In combination, these tools enable one to go from near-raw medical imaging data to a BOLD decoding experiment.

Data organization expected by RKRNS

RKRNS makes several assumptions about data organization.

Contributions of the package

RKRNS includes:

Decoding in simulated data

Here, we outline test data that is similar to what one may acquire in a real event related design. However, in this data, there is no need for motion correction.

Example BOLD Data

We adapted methods from the neuRosim (@neuRosim) package in order to simulate event-related BOLD data. This data involves a 4-class decoding task. One seeks to identify brain responses to novel faces, famous faces and their first and second presentations to the subject. We distribute the underlying anatomical BOLD data along with a cortical mask within the package. The data may be loaded by:

library(RKRNS)
fn<-paste(path.package("RKRNS"),"/extdata/111157_mocoref_masked.nii.gz",sep="") 
eximg<-antsImageRead(fn,3)
fn<-paste(path.package("RKRNS"),"/extdata/subaal.nii.gz",sep="") 
mask<-antsImageRead(fn,3)

which gives you both the cortical mask and the example BOLD signal image eximg. ANTsR also provides AAL label (@Tzourio-Mazoyer2002) names via:

data(aal,package='ANTsR')
labs<-1:90

with cortical labs defined by labs.

Run the simulation

We generate simulated bold data as in neuRosim.
The simulateBOLD function results in a bold image in the antsImage class. We subsequently convert this image to a $n\times p$ matrix where there are $n$ time volumes and $p$ voxels in the mask.

bb<-simulateBOLD(option="henson",eximg=eximg,mask=mask)
boldImage<-bb$simbold
mat<-timeseries2matrix( boldImage, bb$mask )

The data frame containing the bold image also contains the design matrix. One should inspect that matrix to get an idea of the structure expected by RKRNS. Row names and column names are expected to exist in the design matrix dataframe.

Simulated BOLD time series

Perform some diagnostic visualizations on the data.

randvox<-sample(1:ncol(mat),4)
plot(ts(mat[,randvox]))

The time series is noisy, as expected.

Frequency filtering

Let's get rid of high frequency noise.

# only filter high frequencies 1/(2*tr)
filtmat<-filterfMRI4KRNS( data.matrix(mat), tr=1, NA, NA, trendfrequency=2, trendfrequency2=NA )
plot(ts(filtmat[,randvox]))

We can eliminate low frequency contamination with polynomial regressors and other denoising techniques at a later step.

Autoregression Correction

armat<-arCorrection( mat )
print(armat$arcoefs)
plot(ts(armat$outmat[,randvox]))

The armat$arcoefs indicate auto-correlation levels of the input data. This function tries to shrink these levels toward zero.

Estimate the hemodynamic response function

Use glmDenoiseR to estimate a FIR HRF.

mypolydegree<-4
glmda<-glmDenoiseR( data.matrix(mat), bb$desmat[,1:4], hrfBasis=NA, 
  hrfShifts = 12, maxnoisepreds=0 , selectionthresh=0.1 , 
  collapsedesign=T, polydegree=mypolydegree, baseshift=0 )

Visualize the HRF

Take a quick look at the estimated hemodynamic response function.

plot(ts(glmda$hrf))

If this doesn't look "normal," then there may be a problem. One might also do an SVD on the run-wise HRFs and derive a more flexible basis set for your study.

Estimate $\beta$ values for each event

We estimate a $\beta^i$ vector for each of $i \in 1 \cdots q$ events. This vector contains a scalar $\beta$ value for each voxel, i.e. $\beta_j^i$ where $j \in 1 \cdots p$. This example can be found via ?bold2betas.

runs<-bb$desmat$Run; 
btsc<-bold2betas( boldmatrix=data.matrix(mat), 
      designmatrix=bb$desmat[,1:4], baseshift=0,
      blockNumb=runs, maxnoisepreds=0, hrfBasis=glmda$hrf,
      hrfShifts=0, polydegree=mypolydegree, selectionthresh=0.2 )
if ( zscore ) {  # z-score the betas
  betameans<-colMeans(btsc$eventbetas)
  betasds<-apply(btsc$eventbetas,FUN=sd,MARGIN=2)
  btsc$eventbetas<-(btsc$eventbetas-betameans)/betasds
}
mylabs<-rep("",nrow(btsc$eventbetas))
for ( i in 1:nrow(btsc$eventbetas) ) 
  mylabs[i]<-substr( rownames(btsc$eventbetas)[i],1,2)
mylabs<-as.factor(mylabs)

Estimate $\beta$ values for each event: 2

We recommend that one looks at the estimated hemodynamic response function in order to help set parameters. This is available as output from the bold2betas function (and glmDenoiseR). Key parameters are the basis length for HRF estimation, the maximum number of noise predictors (e.g. maxnoisepreds=4 or maxnoisepreds=2:10 to search a range), the polynomial degree and the selectionthresh which should be kept below 0.5. Real data processing should proceed similarly. The latter part of the sample code extracts labels from the rownames of the output event beta dataframe. This will set us up for decoding.

Decode from event-wise $\beta$

Use effect size to select a subset of the full voxel matrix.

# sample for training
inds<-sample(1:nrow(btsc$eventbetas),size=round(nrow(btsc$eventbetas)*3./4.))
# basic voxel selection using effect size & classification
zz<-apply(btsc$eventbetas,FUN=mean,MARGIN=2)
zze<-zz/apply(btsc$eventbetas,FUN=sd,MARGIN=2)
th<-0.4
ff<-which( abs(zze) > th )
mydf<-data.frame( lab=mylabs,  vox=data.matrix(btsc$eventbetas)[,ff])
mdl<-svm( lab ~., data=mydf[inds,])
err<-sum(mydf[-inds,]$lab==predict( mdl, newdata=mydf[-inds,]))/nrow(mydf[-inds,])
print(paste("NPredVox",length(ff),"Correct",err*100))

The "high effect size" voxels are sent to a naive support vector machine to do training (on a fraction of the data) and testing on the left out data.

Cross-validated effect size

Effect size increases when a measurement has low variance and a high value. The effect size of cross-validated beta estimates may aid voxel selection.

# ccnuis<-compcor( mat, 2 ) # a denoising option not used here ....
stb<-stableEventResponse( mat,  bb$desmat[,1:4], bb$desmat$Run,  
            polydegree=mypolydegree, hrf=btsc$hrf )

Voxels of clustered effect size

Use a threshold to find clusters of voxels with large effect.

stb2<-abs(stb)
stb2[ stb2 < 5 ]<-0
ee<-eigSeg(mask, matrixToImages( stb2, mask) )
ImageMath(3,ee,'ClusterThresholdVariate',ee,mask,5)
ff<-which( ee[mask==1] > 0 )
ofn<-"Figs/eigseg.jpg"
plotANTsImage( eximg, list(ee), slices='12x56x1',
  thresh='1x4',color='red',outname=ofn)

The eigSeg function labels each region of the mask with the class that has the largest effect size.

Should return the magnitude of effects ... FIXME colors above.

eigSeg Viz

Decode with effect size clusters

Use only those voxels that exist within clusters and have maximum response for some label class.

mydf<-data.frame( lab=mylabs,  vox=data.matrix(btsc$eventbetas)[,ff])
mdl<-svm( lab ~., data=mydf[inds,])
err<-sum(mydf[-inds,]$lab==predict( mdl, newdata=mydf[-inds,]))/nrow(mydf[-inds,])
print(paste("EigSegPredVox",length(ff),"Correct",err*100))

Time series per stimulus

Investigate the time series for each stimulus.

ff<-which( ee[mask==1] == 1 )
stimbold1<-ts( scale(rowMeans(mat[,ff]) ) )
ff<-which( ee[mask==1] == 2 )
stimbold2<-ts( scale(rowMeans(mat[,ff]) ) )
ff<-which( ee[mask==1] == 3 )
stimbold3<-ts( scale(rowMeans(mat[,ff]) ) )
ff<-which( ee[mask==1] == 4 )
stimbold4<-ts( scale(rowMeans(mat[,ff]) ) )
plot(cbind(stimbold1,stimbold2,stimbold3,stimbold4))

Sparse canonical correlation between BOLD and stimuli

CCA maximizes $PearsonCorrelation( XW^T, ZY^T )$ where $X, W$ are as above and $Z$ and $Y$ are similarly defined. CCA optimizes the matrices $W, Y$ operating on $X, Z$ to find a low-dimensional representation of the data pair $( X , Z )$ in which correlation is maximal. Following ideas outlined in @Dhillon2014 and @Avants2014, this method can be extended with sparsity constraints that yield rows of $W, Y$ with a controllable number of non-zero entries.

SCCAN voxel selection & classification

Set up the CCA by pairing the beta matrix with the design matrix.

ccamats<-list( data.matrix(btsc$eventbetas)[inds,] , 
               data.matrix(bb$desmat[btsc$eventrows,1:4])[inds,] )

Initialize SCCAN

Use the SVD of the beta matrix to initialize sparse CCA.

initcca<-t( svd( btsc$eventbetas, nu=0, nv=10 )$v )
initcca<-initializeEigenanatomy( initcca, mask=mask, nreps=1 )$initlist
nv<-length(initcca)

These r nv vectors initialize the sparse optimizer in a good place.

Supervised clustering with SCCAN

The design matrix guides the dimensionality reduction performed on the $\beta$ matrix.

mycca<-sparseDecom2( inmatrix=ccamats, initializationList=initcca,
  sparseness=c( -0.001, -0.95 ), nvecs=nv, its=10, cthresh=c(250,0),
  uselong=0, smooth=0.0, mycoption=1, inmask=c(mask,NA) )
ccaout<-(data.matrix(imageListToMatrix( mycca$eig1, mask )))
ff<-which(colSums(abs(ccaout))>1.e-4)

We also count the non-zero voxels which cover r length(ff)/ncol(mat)*100% of the brain.

Predictive model

Given CCA solution matrix $W$, one may employ the low-dimensional representation, $XW^T$, in multi-label classification. Currently, we employ SVM or random forests as multi-label learners for the problem:

$$L_i = f( XW^T ),$$

that is, learning a (sentence) label function from the BOLD data.

Predictive model: 2

mydf<-data.frame( lab=mylabs,  
  vox=data.matrix(btsc$eventbetas) %*% t(ccaout) )
mdl<-svm( lab ~., data=mydf[inds,])
err<-sum(mydf[-inds,]$lab==predict( mdl, newdata=mydf[-inds,]))/nrow(mydf[-inds,])
print(paste("CCA: Correct",err*100))
# here is another approach ... use cca to transform BOLD signal to stimuli 
mydf3<-data.frame( lab=mylabs, 
  vox=data.matrix(btsc$eventbetas) %*% t(ccaout) %*% t(mycca$eig2) )
mdl2<-svm( lab ~., data=mydf3[inds,])
err<-sum(mydf3[-inds,]$lab==predict( mdl2, newdata=mydf3[-inds,]))/nrow(mydf3[-inds,])
print(paste("CCA2: Correct",err*100))

The input predictors are both clustered and sparse.

Quick look at CCA results

Rescale the cca results and make a picture.

for ( img in mycca$eig1 ) {
  img[mask==1]<-abs(img[mask==1])
  img[mask==1]<-img[mask==1]/max(img[mask==1])
}
ofn<-"Figs/temp.jpg"
plotANTsImage( eximg, mycca$eig1, slices='12x56x1' , 
 thresh='0.01x1.1', color=rainbow(nv), outname=ofn)

CCA results in image space

Interpreting SCCAN results

Confusion matrix

heatmap of co-occurrence of predictions

library(caret) 
truth<-mydf[-inds,]$lab
pred<-predict( mdl, newdata=mydf[-inds,])
xtab <- table(pred, truth)
cm<-confusionMatrix(xtab)
print(cm$table)

Confusion matrix 2

pheatmap(cm$table, cluster_rows = F, cluster_cols = F)

Anatomical coordinates of SCCAN predictors

ANTs propagates AAL labels to the cortex (@Avants2014a,@Tustison2014)

fn<-paste(path.package("RKRNS"),"/extdata/111157_aal.nii.gz",sep="") 
aalimg<-antsImageRead(fn,3)
ccaanat<-list()
for ( img in mycca$eig1 ) {
  nzind<-img[ mask == 1 ] > 0 
  aalvals<-aalimg[ mask == 1 ][ nzind ]
  aalmax<-which.max( hist( aalvals, breaks=1:90, plot=FALSE )$counts )+1 
  ccaanat<-lappend( ccaanat, aal$label_name[aalmax] )
}

The SCCAN predictors include: r unlist(ccaanat).

How much of the known network do we actually find?

Associating classes to SCCAN predictors

Recalling: CCA maximizes $PearsonCorrelation( XW^T, ZY^T )$, we can study matrix $Y$ which contrasts or combines columns of the design matrix.

rownames(mycca$eig2)<-levels(mylabs)
temp<-(mycca$eig2)
temp[ abs(mycca$eig2) < 0.03 ]<-0

Associating classes to SCCAN predictors: 2

pheatmap(temp)

Conclusions

With the current RKRNS, one may:

The package needs evaluation at this level of detail on real data.

See RKRNS for all source code and documentation and RKRNS-talk for html slides.

Real data processing with RKRNS

We recommend that you motion correct the data before you use RKRNS tools. This should be done once and done well. Good approaches exist in ANTsR which yield both motion matrices and relevant summary measurements such as FD and DVARS. See ?antsPreprocessfMRI for a simplified utility function. This function could be used on each run of an experiment and the results stored in organized fashion for later use.

Motion correction

To motion correct your data, one might run:

# get an average image
averageImage <- getAverageOfTimeSeries( boldImage )
motionCorrectionResults <- motion_correction(boldImage, 
   fixed = averageImage, moreaccurate = 0 )

Set the moreaccurate flag to 1 or 2 for usable (not test) results. You might also estimate FD and DVARS from these results. One might use antsPreprocessfMRI to get these directly. Note, however, to turn this function's frequency filtering off if you want to do decoding.

FD and DVARS

hey hey ...

References



stnava/RKRNS documentation built on May 30, 2019, 7:21 p.m.