CIFTI in R, MSC single-subject

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

Setup System Color for Plot

parlabel <- data.frame(parcel_num=sort(u_parcel), 
                       community=NA, 
                       comm_label=NA, 
                       comm_shortlabel=NA)

plotlabel <- read.csv("../../data/systemlabel_MSC.txt", header=F,
                          col.names = c("community","comm_label","color","comm_shortlabel"))

for(i in 1:length(u_parcel)){
  comm_i <- unique(comm[which(parcel==sort(u_parcel)[i])])
  parlabel$community[i] <- comm_i
  parlabel$comm_label[i] <- as.character(plotlabel$comm_label[is.element(plotlabel$community, comm_i)])
  parlabel$comm_shortlabel[i] <- as.character(plotlabel$comm_shortlabel[is.element(plotlabel$community, comm_i)])
}

superheat::superheat(X = z, 
                     y.axis.reverse = TRUE,
                     membership.rows = parlabel$comm_shortlabel,
                     membership.cols = parlabel$community,
                     left.label.col=plotlabel$color,
                     bottom.label.col=plotlabel$color,
                     extreme.values.na = FALSE,
                     heat.lim = c(-.2, .4), 
                     heat.pal = rev(brewer.rdylbu(100)),
                     heat.pal.values = c(0, 0.15, 0.25, 0.75,1),
                     title="Node x Node Correlation Matrix Organized by Systems")

Participation Coefficient (4% edge density)

• Participation coefficient measures a node’s connections within its community proportion to its conncetion to the entire network.

if (!require("NetworkToolbox", character.only=T, quietly=T)) {
  devtools::install_github("AlexChristensen/NetworkToolbox")
}

## 
## Attaching package: 'NetworkToolbox'

## The following objects are masked from 'package:igraph':
## 
##     betweenness, closeness, degree, diversity, strength,
##     transitivity

## The following object is masked from 'package:dplyr':
## 
##     desc

library(NetworkToolbox)

p <- participation(A = z4, comm = parlabel$community)

# Each node's PC calculated using positive & negative edge
# Negative edges were taken out in previous steps, so PC caculated with all-edges and positive-edges are the same. os
head(p$overall)

## [1] 0.8375559 0.3965046 0.6669420 0.5441093 0.6031861 0.3480727

head(p$positive)

## [1] 0.8375559 0.3965046 0.6669420 0.5441093 0.6031861 0.3480727

# PC based on negative edges should all be zero (not usable in this case).
head(p$negative)

## [1] 0 0 0 0 0 0

# Coloring nodes based on PC.
gcol <- grey.colors(n=nrow(z))
plot(net, layout=layout_with_fr, vertex.label=NA, vertex.size=5, 
     vertex.color=gcol[order(p$positive)], alpha=.6)

Distribution of participation coefficient across entire network and subnetwork

parlabel$pc_4td <- p$positive

ggplot(parlabel, aes(x=pc_4td)) +
  geom_histogram(bins=20) +
  xlab("Participation Coefficient (4%)") +
  ggtitle("Participation Coefficient (4%) across entire network") +
  theme_bw()

ggplot(parlabel, aes(x=pc_4td)) +
  facet_wrap(~comm_shortlabel) +
  geom_histogram(bins = 20) +
  xlab("Participation Coefficient (4%)") +
  ggtitle("Participation Coefficient (4%) across each sub-network") +
  theme_bw()

Use a density plot to visualize distributions where there are small number of nodes.

ggplot(parlabel, aes(x=pc_4td)) +
  facet_wrap(~comm_shortlabel) +
  geom_density() +
  xlab("Participation Coefficient (4%)") +
  ggtitle("Participation Coefficient (4%) across each sub-network") +
  theme_bw()

• 00Bd (UnAssigned) are parcels that don’t belong to any community, and likely have very little connections (low degree)
• Histogram of degree distribution (4% edge density) shows that nodes that are UnAssinged have very low degree

parlabel$degree_4td <- degree(A = z4) # calculate degree for each node at 4% edge density

ggplot(parlabel, aes(x=degree_4td)) +
  facet_wrap(~comm_shortlabel) +
  geom_histogram(bins = 20) +
  xlab("Degree (4%)") +
  ggtitle("Degree (4%) across each sub-network") +
  theme_bw()