Prepare data

Read in cifti files

# Load CIFTI data files
cii <- read_cifti(mscfile, drop_data = FALSE, trans_data = T) 
parcel <- as.matrix(read_cifti(parcel_file)$data)
comm <- as.matrix(read_cifti(comm_file)$data)

u_parcel <- unique(parcel)
u_parcel <- u_parcel[u_parcel!=0] # Remove parcel 0 and order parcel by their number

Make brainstructure index

The brainstrucure index lets us filter out anatomical structures based on an index (this mirrors the cifti packages in MATLAB). As of June 2019, the cifti package on CRAN would ignore some subcortical labels, so make sure to use the development version on github.

cii$brainstructureindex <- as.matrix(NA, dim(cii$data)[1])
for(i in 1:length(cii$BrainModel)){
  startindx <- attributes(cii$BrainModel[[i]])$IndexOffset + 1
  endindx <- attributes(cii$BrainModel[[i]])$IndexOffset + 
             attributes(cii$BrainModel[[i]])$IndexCount
  
  cii$brainstructureindex[startindx:endindx] <- i
}

Check dimension of cifti data (volume/frame x vertices)

  • Dimension of BOLD, Parcel, and Community:
dim(cii$data) # ~ 64k vertices, includes subcortical volumes
## [1] 65890   818     1
dim(parcel)   # surface only, excluded medial wall
## [1] 59412     1
dim(comm)     # surface only, excluded medial wall
## [1] 59412     1

Mismatch can be due to inclusion of subcortical/medial wall

CIFTI data contains the surface (cortx left and right) and subcortical structures based on volumetric data. The labels should contain the left & right coritcal surface, AND subcortical labels.

What are the labeled brain structures in the BOLD cifti file?

cifti_brain_structs(cii)
##  [1] "CIFTI_STRUCTURE_CORTEX_LEFT"      
##  [2] "CIFTI_STRUCTURE_CORTEX_RIGHT"     
##  [3] "CIFTI_STRUCTURE_ACCUMBENS_LEFT"   
##  [4] "CIFTI_STRUCTURE_ACCUMBENS_RIGHT"  
##  [5] "CIFTI_STRUCTURE_AMYGDALA_LEFT"    
##  [6] "CIFTI_STRUCTURE_AMYGDALA_RIGHT"   
##  [7] "CIFTI_STRUCTURE_CAUDATE_LEFT"     
##  [8] "CIFTI_STRUCTURE_CAUDATE_RIGHT"    
##  [9] "CIFTI_STRUCTURE_CEREBELLUM_LEFT"  
## [10] "CIFTI_STRUCTURE_CEREBELLUM_RIGHT" 
## [11] "CIFTI_STRUCTURE_HIPPOCAMPUS_LEFT" 
## [12] "CIFTI_STRUCTURE_HIPPOCAMPUS_RIGHT"
## [13] "CIFTI_STRUCTURE_PALLIDUM_LEFT"    
## [14] "CIFTI_STRUCTURE_PALLIDUM_RIGHT"   
## [15] "CIFTI_STRUCTURE_PUTAMEN_LEFT"     
## [16] "CIFTI_STRUCTURE_PUTAMEN_RIGHT"    
## [17] "CIFTI_STRUCTURE_THALAMUS_LEFT"    
## [18] "CIFTI_STRUCTURE_THALAMUS_RIGHT"

Subcortical data are included. Since subcortical data are not sorted into the community assignments provided by MSC data, only the cortical surface data will be extracted and analyzed.

cdata <- as.matrix(cii$data[cii$brainstructureindex==1 | cii$brainstructureindex==2,,])

Check the dimension of the BOLD data again

dim(cdata)
## [1] 59412   818

Remove motion contaminated data

tmask <- read.table(tmask_file)$V1
ctmask <- cdata[,as.logical(tmask)]
sprintf("Number of high-motion frames = %s (%s%% removed)", sum(tmask==0), round(sum(tmask==0)/length(tmask)*100))
## [1] "Number of high-motion frames = 295 (36% removed)"

Extract mean time series from each parcel into a matrix (node x volume/frame)

tp <- matrix(0, length(u_parcel), sum(tmask)) # initialize empty matrix

for(i in 1:length(u_parcel)){               
  tp[i,]<- colMeans(ctmask[which(parcel==u_parcel[i]),])
}

tp <- tp[order(u_parcel),] # Order matrix by parcel number

Plot processed mean time series of each node

superheat::superheat(tp,
                     heat.lim = c(-20, 20), 
                     heat.pal = c("black","white"),
                     grid.hline = FALSE,
                     grid.vline = FALSE,
                     title="Mean Time series of each parcel (high-motion frames removed)")