RestingBOLD.RmdThis document provides some examples illustrating how ANTsR may be used to analyze resting state BOLD fMRI data using established methodology. The focus is on processing the BOLD signal, here we start with the the following data
A BOLD fMRI time-series image
A brain mask image
A tissue-segmentation that identifies (at least): CSF, gray matter, white matter
A set of labels identifying anatomical regions of interest (the network to analyze)
library(ggplot2)
library(igraph)
library(pracma)
library(dplyr)
library(mFilter)
fdthresh = 0.1 # threshold for throwing out bad data, typically between 0.1 and 0.5
# based on ACT output or sample data
id = "BAvants"
pre = paste("~/rsfTest/",id,"/",sep='')
fn = path.expand( paste( pre, "rsfMRI0/",id,"_rsfMRI0.nii.gz", sep='' ) )
testData = 1
if ( file.exists( fn ) )
{
fdthresh = 0.2 # threshold for throwing out bad data, typically between 0.1 and 0.5
testData = 3
imgIn = antsImageRead( fn )
meanbold = getAverageOfTimeSeries( imgIn )
mask = getMask( meanbold )
sfn = path.expand( paste( pre, "/act/", id,"_STRUCTURAL_BrainSegmentation.nii.gz", sep='' ) )
seg = antsImageRead( sfn )
t1fn = path.expand( paste( pre, "/act/", id,"_STRUCTURAL_BrainSegmentation0N4.nii.gz", sep='' ) )
t1 = antsImageRead( t1fn )
} else { # use sample data
imgIn = antsImageRead(getANTsRData("rsbold"))
meanbold = getAverageOfTimeSeries( imgIn )
mask = antsImageRead(getANTsRData("rsboldmask"))
seg = antsImageRead(getANTsRData("rsboldseg"))
t1 = antsImageClone( seg )
t1[ t1 > 3 ] = 2
}
t1brain = t1 * thresholdImage( seg, 1, 6 )
if ( ! exists("boldmap") )
boldmap = antsRegistration( meanbold, t1brain, typeofTransform='SyNBoldAff' )
if ( ! exists("mnimap") )
{
mni = antsImageRead( getANTsRData( "mni" ) )
mnimap = antsRegistration( t1brain, mni, typeofTransform='SyN' )
}
mni2boldmaps = c( boldmap$fwdtransforms, mnimap$fwdtransforms )
mni2boldmapsInv = c( mnimap$invtransforms , boldmap$invtransforms )
mni2bold = antsApplyTransforms( meanbold, mni, mni2boldmaps )
seg2bold = antsApplyTransforms( meanbold, seg, boldmap$fwdtransforms, interpolator = "NearestNeighbor" )
plot( meanbold , boldmap$warpedmovout %>% iMath("Canny", 10, 1, 1) )
plot( meanbold , mni2bold %>% iMath("Canny", 10, 1, 1) )
plot( meanbold , maskImage( seg2bold, seg2bold, 2 ) )Obtaining these for your own data set is a non-trivial matter, but will be the topic of a future document as the process is the same for both resting and task-based BOLD.
The processing here is largely based upon a recent review of methods for dealing with motion in resting fMRI (Power et al. 2014).
The Preprocessing section includes methods used for processing fMRI for a variety of purposes. This is followed by a section on Connectivity processing that includes steps specific to the processing of resting state data. Finally, a section on building graphs and calculating graph metrics is presented.
One step that is omitted here is slice timing correction as this can vary greatly between acquisition sequences. The included steps are:
Removal of initial time points that occur before magnetization steady state is reached
Rigid registration of time points to correct for subject head motion
Plotting of motion parameters for inspection
Identification of “bad” time points for exclusion
The first step is the removal of pre steady state time points. It is typical to exclude any data obtained during the first 10 seconds as shown below. The choice of 10s is based on an informal review of current literature of 3T human data (Power et al. 2012). For other applications, be sure to check the relevant literature.
# Find first steady state timepoint
tr = antsGetSpacing(imgIn)[4]
steady = floor(10.0 / tr) + 1
# Global signal before cropping (save for visualization)
origmean = getAverageOfTimeSeries( imgIn )
fullmean = rowMeans(timeseries2matrix(imgIn, mask))
allTimes = dim(imgIn)[4]
# Eliminate non steady-state timepoints
img = cropIndices(imgIn, c(1,1,1,steady), dim(imgIn) )In the plot below, mean global signal in the brain is plotted with a red box indicating the points at which the system has non yet reached the steady state. Data in this range are discarded.
The temporal mean of the original time series data:
To correct for rigid body motion that occurs during acquisition:
Find the mean volume of all time points. This is done with
apply.antsImage which is an extension of the R
method apply with additional functionality to maintain
image header info integrity.
Align all time-points to the mean. This is accomplished with
antsMotionCalculation, the fixed parameter is
used to set the reference image to which all time-points are aligned and
the txtype parameter indicates the type of transform to be
estimated. The default for txtype is “Affine”, but for this
type of analyses it is typical to use “Rigid”.
To help ensure an accurate fit, 3 iterations of the above steps are used
Examine motion correction parameters for quality control.
Identify “bad” time points for removal
if ( ! exists("moco") )
{
meanbold <- getAverageOfTimeSeries( img )
for ( i in 1:testData )
{
moco <- antsrMotionCalculation( img, fixed=meanbold, typeofTransform="Rigid" )
meanbold = getAverageOfTimeSeries( moco$moco_img )
}
}It can also be informative to plot the data as a matrix, where each row is the time-series for a voxels. Due to the large number of voxels however, using just a sample of the voxels is much faster
nVox = length(which(as.array(mask)==1))
vox = sample(1:nVox, 1000)
invisible(plot(as.antsImage( t(timeseries2matrix(img,mask)[,vox]))))
invisible(plot(as.antsImage( t(timeseries2matrix(moco$moco_img,mask)[,vox]))))Plotting the registration parameters from the motion correction provides a qualitative feel for how much motion is in the data. In addition to the registration parameters, we plot the mean framewise displacement, which measures the average displacement of voxels, between consecutive time points.
Another metric that is typically plotted is the DVARS (D for temporal derivative, VARS for RMS variance over voxels) which illustrates the BOLD signal change across the brain between consecutive time points. First the BOLD signal is adjusted to have mean=1000 so that 10 units of BOLD value change = 1% signal change. While not necessary for a single-subject examination, this helps in inter-subject comparisons.
The motion parameters are often used to identify “bad” timepoints. A mean framewise displacement greater than 0.2mm is a common threshold used in human data. For illustrative purposes, we will use a threshold of 0.1 mm. Because the displacement is a measure of motion between two timepoints, both timepoints associated with the displacement are marked as bad.
The steps here are:
Demean and detrend the data
Regress out nuisance parameters
Frequency filtering
Spatial smoothing
The time-series data is detrended while excluding the bad timepoints identified earlier. The global signal (mean signal over the whole brain) is used to illustrate the effect of the demeaning & detrending.
global_pre <- rowMeans(timeseries2matrix(img, mask))
global_moco <- rowMeans(timeseries2matrix(moco$moco_img, mask))
boldMat = timeseries2matrix(moco$moco_img, mask)
boldMat[goodtimes,] = detrend(boldMat[goodtimes,])
boldMat[badtimes,] = NA
global_moco_detrend = rowMeans(boldMat)
global_pre[badtimes] = NA
global_moco[badtimes] = NASome typical nuisance parameters are
detrended motion parameters, their squares, and the derivatives of both
mean signal in white matter & it’s derivative
mean signal in CSF & it’s derivative
physiologocial noise estimated via compcor
global mean signal in brain & it’s derivative
There are two camps when it comes to global signal, some leave it in, others regress it out. It is unclear which is best. Here we will include it as a nuisance parameter as done in the paper on which the methods are based. This should not be interpreted as an implied endorsement of one camp over the other. The values at “bad” timepoints are interpolated with splines so that derivatives may be calculated without “spreading” the influence of the bad time points.
# white matter is labeled as 3
wmMask = seg2bold*1*mask
wmMask[ wmMask != 3] = 0
wmMask[ wmMask == 3 ] = 1
wmMask = iMath( wmMask, "ME", 1)
wmVox = which(subset(wmMask, mask > 0 )==1)
wmMean = rowMeans(boldMat[,wmVox])
# CSF is labeled as 1
csfMask = seg2bold*1
csfMask[ csfMask != 1] = 0
#csfMask = iMath( csfMask, "ME", 1)
csfVox = which(subset(csfMask, mask > 0)==1)
csfMean= rowMeans(boldMat[,csfVox])
#csfMean = rowMeans(timeseries2matrix(detrendImg, csfMask))
globalMean = rowMeans(boldMat)
ncompcor = 4
compcorTemp = compcor(boldMat[goodtimes,], ncompcor = ncompcor)
compcorNuis = matrix(0, nTimes, ncompcor )
compcorNuis[goodtimes, ] = compcorTemp
compcorNuis[badtimes, ] = NA
colnames( compcorNuis ) = paste("compcor",1:ncol(compcorNuis), sep='' )
tissueNuis = cbind(globalMean, wmMean, csfMean)
if ( length(badtimes) > 0 ) {
for ( v in c(1:dim(tissueNuis)[2]) ) {
tissueInterp = spline( c(1:nTimes)[goodtimes], tissueNuis[goodtimes,v],
method='natural', xout=badtimes )$y
tissueNuis[badtimes,v]=tissueInterp
}
}
tissueDeriv = rbind( rep(0,dim(tissueNuis)[2]), diff(tissueNuis,1) )The nuisance parameters are now regressed out the signal. This is illustrated by looking at the mean signal in the cortex before and after the regression.
mocoNuis = cbind(reg_params, reg_params*reg_params)
mocoNuis = detrend(mocoNuis)
mocoDeriv = rbind( rep(0,dim(mocoNuis)[2]), diff(mocoNuis,1) )
nuissance = cbind( mocoNuis, mocoDeriv, tissueNuis, tissueDeriv, compcorNuis )
boldMat[goodtimes,] <- residuals( lm( boldMat[goodtimes,] ~ nuissance[goodtimes,] ) )The next step is frequency filtering. However, first we want to fill in the “bad” timepoints with interpolated data. This is to avoid artifacts that would result from having non-evenly sampled time-series data. This bad timepoints will be again removed after the frequency filtering. Frequencies with the range of 0.009 Hz - 0.08 Hz are retained. In some cases it may be interesting to examine smaller subranges of frequencies to look for phenomena that occur at specific frequencies.
if ( length(badtimes) > 0 ) {
for ( v in c(1:nVox) ) {
boldMat[badtimes,v]=spline( c(1:nTimes)[goodtimes], boldMat[goodtimes,v],
method='natural', xout=badtimes )$y
}
}
# save interpolated values for plotting
ctxMeanSpline = rowMeans(boldMat[,ctxVox])
fboldMat <- frequencyFilterfMRI( boldMat, tr=tr, freqLo=0.009, freqHi=0.08, opt="trig" )
# save filtered values for plotting
ctxMeanFiltered = rowMeans(fboldMat[,ctxVox])
ctxMeanFiltered[badtimes] = NASmoothing should be applied to all spatial dimensions, but the time dimension should be left alone. It is common to smooth with a Gaussian kernel with FWHM=6.0mm.
fimg = matrix2timeseries( img, mask, fboldMat )
sptl = sqrt( sum( antsGetSpacing(img)[1:3]^2 ))
simg = smoothImage( fimg, c(rep(sptl,3),0), FWHM=TRUE )
sboldMat = timeseries2matrix(simg, mask)First, a set of ROIs is needed. While each individual voxel could be treated as an ROI, it is more common to define ROIs that contain many voxels and represent regions for which there is some a priori knowledge regarding the functional network to which each region belongs. The provided network definition used here lists the center point of each ROI so we must first create an image where each ROI is a sphere of radius=5mm centered on the provided point.
data(powers_areal_mni_itk)
pts = antsApplyTransformsToPoints( 3, powers_areal_mni_itk, transformlist = mni2boldmapsInv )
pts[ , 4:ncol(pts) ] = powers_areal_mni_itk[ , 4:ncol(pts) ]
labelImg = makePointsImage( pts, mask )
nPts = dim(pts)[1]
plot( meanbold, labelImg, axis=3, nslices=30, ncolumns=10,
window.overlay = c( 1, max(labelImg) ) )Now that we have roi labels, we want to find the mean time signal for each ROI and then find the correlation matrix that give the correlation between each of these ROI signals. We need to be careful in the case where ROIs did not map into the image space and thus have 0 voxels.
labelMask = labelImg*1
labelMask[labelMask > 0] = 1
labelMask[mask == 0] = 0
labelVox = which(subset(labelMask, mask > 0)==1)
labeledBoldMat = boldMat[goodtimes,labelVox]
labels = labelImg[labelMask > 0]
nLabels = max(labels)
roiMat = matrix(0, nrow=dim(labeledBoldMat)[1], ncol=nLabels)
for ( i in c(1:nLabels) ) {
if (length(which(labels==i)) > 1 ) {
roiMat[,i] = rowMeans(labeledBoldMat[,(labels==i)])
}
}
nActualTimes = dim(roiMat)[1]Plot of the ROI averaged signals
If system labels are provided for the ROIs, it is possible to look at the mean and standard deviation for the BOLD signal within a system. Comparing the mean signal between systems provides some hints about how much those systems are working together, while the standard deviation ribbons indicate how cohesively the components of a given system are working.
systemNames = unique(pts$SystemName)
nSystems = length(systemNames)
sysMatMean = matrix(0, nrow=dim(labeledBoldMat)[1], ncol=nSystems)
sysMatSD = matrix(0, nrow=dim(labeledBoldMat)[1], ncol=nSystems)
systems = pts$SystemName[labels]
for ( i in 1:nSystems ) {
sys = systemNames[i]
sysIdx = which(systems==sys)
if ( length(sysIdx) > 0)
{
sysMatMean[,i] = rowMeans(labeledBoldMat[,sysIdx])
sysMatSD[,i] = apply(labeledBoldMat[,sysIdx], 1, sd)
}
}Create correlation matrix
missingROIs = which(colMeans(roiMat)==0)
goodROIs = (1:nLabels)
if ( length(missingROIs) > 0 ) {
goodROIs = goodROIs[-missingROIs]
}
connMat = suppressWarnings(cor(roiMat))
diag(connMat) = rep(0, length(diag(connMat)) )
if ( length(missingROIs) > 0 ) {
connMat[missingROIs,] = 0
connMat[,missingROIs] = 0
}Create mutual information matrix
if ( usePkg( "entropy" ) ) {
miMat = connMat * 0
nb = 10 # can have a strong influence on results
rr1 = rr2 = range( roiMat )
for ( i in goodROIs )
for ( j in goodROIs[ goodROIs > i ] ) {
if ( i != j) {
rr1 = range(roiMat[,i])
rr2 = range(roiMat[,j])
y2d = discretize2d( roiMat[,i], roiMat[,j] , nb, nb, r1=rr1, r2=rr2 )
# miMat[ i, j ] = mi.Dirichlet( y2d, a=1/8)
miMat[ i, j ] = mi.empirical( y2d )
}
}
image( miMat + t(miMat) )
}Networks are often created by using a constant density. For example, to create a network with density=0.1, we binarize the correlation matrix to retain 10% of the edges, favoring the egdes with the highest correlation values. The resulting adjacency matrix is then used to create a graph.
density = 0.1
nEdges = length(upper.tri(connMat))*density
thresh = sort( connMat[upper.tri(connMat)], decreasing=T)[nEdges]
adj = 1*(connMat >= thresh)
bingraph = graph.adjacency(adj, mode="undirected", weighted=NULL)
components = clusters(bingraph)
maxID = which(components$csize == max(components$csize))[1]
adj[components$membership!=maxID,] = 0
adj[,components$membership!=maxID] = 0
bingraph = graph.adjacency(adj, mode="undirected", weighted=NULL)
invisible(plot(as.antsImage(adj)))A more interesting way to visualize the adjacency matrix is to color the components by the system to which they belong. Of course, this requires that the ROIs used to define the network, include system identifies. In the plot, connections within a system are color coded while connections between systems are gray. For this type of plot, the order of the rows are sorted, so that nodes in the same system are clustered together.
It is also possible to export the data to a graphml file for visualization with an application such as GEPHI. In the example, the node colors are retained and edges for intra-system connections are colored the same as their nodes, while inter-system edges are gray.
# Retain only the largest connected component
bingraph = graph.adjacency(adj, mode="undirected", weighted=NULL)
components = clusters(bingraph)
maxID = which(components$csize == max(components$csize))[1]
adj[components$membership!=maxID,] = 0
adj[,components$membership!=maxID] = 0
graph = graph.adjacency( adj, mode="undirected", weighted=NULL )
# Set node colors
graph = set.vertex.attribute(graph, "r", index=V(graph), value=as.double(pts$r))
graph = set.vertex.attribute(graph, "g", index=V(graph), value=as.double(pts$g))
graph = set.vertex.attribute(graph, "b", index=V(graph), value=as.double(pts$b))
# Set edge colors
edges = get.edges( graph, E(graph) )
nEdges = dim(edges)[1]
er = rep(200, nEdges)
eg = rep(200, nEdges)
eb = rep(200, nEdges)
# colors for intra-system connections
# gray for inter-system connections
for ( e in c(1:nEdges) )
{
if ( pts$SystemName[edges[e,1]] == pts$SystemName[edges[e,2]] )
{
er[e] = pts$r[edges[e,1]]
eg[e] = pts$g[edges[e,1]]
eb[e] = pts$b[edges[e,1]]
}
}
graph = set.edge.attribute(graph, "r", index=E(graph), value=as.double(er))
graph = set.edge.attribute(graph, "g", index=E(graph), value=as.double(eg))
graph = set.edge.attribute(graph, "b", index=E(graph), value=as.double(eb))
# uncomment line below to write out graph
# write.graph(graph, "network.graphml", format="graphml", prefixAttr=FALSE)Having derived a graph representation of the resting brain network, the next logical step is to examine graph-metrics that encapsulate various properties of the network. These are describe in further detail in a review of graph metrics for studying connectivity in the brain (Rubinov and Sporns 2010). Two types of metrics will be examined.
Node measures
Global network measures
The measures are made independently for each node/vertex in the graph. The following is a non-comprehensive list of possible node metrics
Degree - the number of connections that include the node
Clustering coefficient - local neighborhood connectivity
Path length - mean shortest distance between this node and all others
Local efficiency - measures “closeness” off nodes in a neighborhood
Page-rank - Google page rank measure
graph = graph.adjacency( adj, mode="undirected", weighted=NULL )
deg = degree(graph)
deg[deg==0] = NA
pathsmat = shortest.paths(graph, weights=NA)
pathsmat[!is.finite(pathsmat)] = NA
paths = rowMeans(pathsmat, na.rm=TRUE)
paths[paths==0] = NA
clust = transitivity(graph, type="local")
clust[deg < 2] = NA
pager = page.rank(graph)$vector
pager[deg < 2] = NA
# from http://pastebin.com/XqkEYtJS
leff <- numeric(length(deg))
goodnodes <- which(deg > 1)
leff[goodnodes] <- sapply(goodnodes, function(x) {
neighbs <- neighbors(graph, v=x)
g.sub <- induced.subgraph(graph, neighbs)
Nv <- vcount(g.sub)
lpaths <- shortest.paths(g.sub, weights=NA)
lpaths <- paths[upper.tri(lpaths)]
pathsup <- lpaths[upper.tri(lpaths)]
2 / Nv / (Nv - 1) * sum(1 / lpaths[which(is.na(lpaths)==FALSE)])
}
)
leff[deg < 2] = NA
leff[which(is.na(deg)==TRUE)] = NAThe measures summarize the entire network in a single measure. All node metrics may be meaned over all nodes to obtain a global metric. Additional metrics include
Global efficiency - closeness of nodes to all other nodes
Clustering coefficient - same as node based, but using entire graph as neighborhood
geff<-1/(shortest.paths(graph))
geff[!is.finite(geff)]<-NA
geff<-mean(geff,na.rm=TRUE)
cc = transitivity(graph)Finally, make a dense map of the default mode network, using regression.
refSignal = sysMatMean[ , systemNames == "Default Mode" ]
# get priors for different networks
if ( ! exists( "networkPriors" ) )
{
networkPriors = getANTsRData("fmrinetworks")
ilist = networkPriors$images
for ( i in 1:length(ilist) )
ilist[[i]] = antsApplyTransforms( meanbold, ilist[[i]], mni2boldmaps )
}
pr = imageListToMatrix( ilist, mask )
refSignal = ( sboldMat %*% t(pr) )[,1]
networkDf = data.frame( ROI=refSignal[goodtimes],
nuissance[goodtimes,grep("compcor",colnames(nuissance))] )
mdl = lm( sboldMat[goodtimes,] ~ . , data=networkDf )
bmdl = bigLMStats( mdl, 1.e-4 )
betas = bmdl$beta.t["ROI",]
betasI = makeImage( mask, betas )
loth = quantile( betas, probs=0.5 )
plot( meanbold, betasI, axis=3, nslices=30, ncolumns=10,
window.overlay = c( 2, max(betas) ) )
# same thing, different denoising strategy
dnz = aslDenoiseR( fboldMat[goodtimes,], refSignal[goodtimes],
polydegree = 2, crossvalidationgroups=8, verbose=T )
networkDf = data.frame(
ROI=refSignal[goodtimes],
dnz$noiseu, dnz$covariates )
mdl = lm( sboldMat[goodtimes,] ~ . , data=networkDf )
bmdl = bigLMStats( mdl, 1.e-4 )
betas = bmdl$beta.t["ROI",]
betasI = makeImage( mask, betas )
loth = quantile( betas, probs=0.5 )
plot( meanbold, betasI, axis=3, nslices=30, ncolumns=10,
window.overlay = c( 2, max(betas) ) )