DGM-RSN10

Install required packages

# install.packages("multdyn")
# install.packages("testit")
# install.packages("ggplot2")
# install.packages("cowplot")
# install.packages("reshape2")
# install.packages("data.table")

Load libraries

library(DGM)
library(ggplot2)
library(cowplot)
library(reshape2)
library(data.table)
library(testit)

Main variables

PATH_HOME = "/home/simon"
PATH_NET  = file.path(PATH_HOME, 'Drive', 'DGM_NetMats10_results')
PATH_NET_BO=file.path(PATH_HOME, 'Drive', 'DGM_NetMats10_results_stepwise')
PATH_TS   = file.path(PATH_HOME, 'Drive', 'HCP900_Parcellation_Timeseries_Netmats_ICAd25')
PATH      = file.path(PATH_HOME, 'Dropbox', 'Data', 'DGM-Sim')
PATH_DATA = file.path(PATH, 'data')
PATH_DEMOG =file.path(PATH_HOME, 'Drive', 'HCP')
PATH_FIG  = file.path(PATH, 'figures')
N      = 500
Nn = 10 # RSNs
N_t    = 1200 # Volumes
N_runs = 4
INFO   = 'NetMats10'
SUBJECTS = as.matrix(read.table(file.path(PATH_TS, 'subjectIDs.txt')))
TR=0.72
labels = read.table(file.path(PATH, 'results', 'RSN10Labels.txt'), header = T)
labels = labels[order(labels$NoNetMats25),]
COMP=labels$NoNetMats25
print(labels)

Demographics

We use the HCP 900 release

Load data

subj = as.matrix(read.table(file.path(PATH_TS, 'subjectIDs.txt')))
subj=subj[1:N]
           
demog=read.table(file.path(PATH_DEMOG, 'unrestricted_schwab_9_19_2017_8_41_11.csv'), 
                 sep = ",", header = TRUE)
idx = is.element(demog$Subject, subj)
assert(sum(idx) == N)
demog = subset(demog, idx)
assert(nrow(demog) == N)

summary(demog$Gender)/N

    F     M 
0.566 0.434 

summary(demog$Age)/N

22-25 26-30 31-35   36+ 
0.206 0.422 0.362 0.010 

Loading network data

# subj_r1 = vector(mode = "list", length = N)
# subj_r2 = vector(mode = "list", length = N)
# subj_r3 = vector(mode = "list", length = N)
# subj_r4 = vector(mode = "list", length = N)
# 
# for (s in 1:N) {
#   subj_r1[[s]] = read.subject(PATH_NET, sprintf("%s_Run_%03d", SUBJECTS[s], 1), Nn)
#   subj_r2[[s]] = read.subject(PATH_NET, sprintf("%s_Run_%03d", SUBJECTS[s], 2), Nn)
#   subj_r3[[s]] = read.subject(PATH_NET, sprintf("%s_Run_%03d", SUBJECTS[s], 3), Nn)
#   subj_r4[[s]] = read.subject(PATH_NET, sprintf("%s_Run_%03d", SUBJECTS[s], 4), Nn)
#   print(sprintf("Subject %03d loaded",s))
# }
# 
# dgm.net10.r1 = dgm.group(subj_r1)
# dgm.net10.r2 = dgm.group(subj_r2)
# dgm.net10.r3 = dgm.group(subj_r3)
# dgm.net10.r4 = dgm.group(subj_r4)
# 
# # save some space
# #dgm.net10.r1$models = NULL
# dgm.net10.r2$models = NULL
# dgm.net10.r3$models = NULL
# dgm.net10.r4$models = NULL
# 
# f=file(file.path(PATH,"results", "DGM-RSN10.RData"))
# save(dgm.net10.r1, dgm.net10.r2, dgm.net10.r3, dgm.net10.r4, file = f, compress = T)
# close(f)
load(file.path(PATH, 'results', 'DGM-RSN10.RData'))

Plot: discount factor delta distribution per node

d.r1=melt(dgm.net10.r1$df_)
d.r2=melt(dgm.net10.r2$df_)
d.r3=melt(dgm.net10.r3$df_)
d.r4=melt(dgm.net10.r4$df_)
set.seed(1980) # because jitter will always create a slightly different image
p1 = ggplot(d.r1, aes(x=as.factor(Var2), y=value)) +
  geom_point(shape=1, color="gray50", size=2, position = position_jitter(width = NULL, height = 0.01)) +
  geom_boxplot(width=0.6) +
  ggtitle(sprintf("Run 1", N)) + ylab("df") + xlab("node")
p2 = ggplot(d.r2, aes(x=as.factor(Var2), y=value)) +
  geom_point(shape=1, color="gray50", size=2, position = position_jitter(width = NULL, height = 0.01)) +
  geom_boxplot(width=0.6) +
  ggtitle(sprintf("Run 2", N)) + ylab("df") + xlab("node")
p3 = ggplot(d.r3, aes(x=as.factor(Var2), y=value)) +
  geom_point(shape=1, color="gray50", size=2, position = position_jitter(width = NULL, height = 0.01)) +
  geom_boxplot(width=0.6) +
  ggtitle(sprintf("Run 3", N)) + ylab("df") + xlab("node")
p4 = ggplot(d.r4, aes(x=as.factor(Var2), y=value)) +
  geom_point(shape=1, color="gray50", size=2, position = position_jitter(width = NULL, height = 0.01)) +
  geom_boxplot(width=0.6) +
  ggtitle(sprintf("Run 4", N)) + ylab("df") + xlab("node")
plot_grid(p1, p2, p3, p4, ncol = 2, nrow = 2)

DFs with parents only

df.10r1 = dgm.net10.r1$df_
df.10r2 = dgm.net10.r2$df_
df.10r3 = dgm.net10.r3$df_
df.10r4 = dgm.net10.r4$df_
df.10r1[t(apply(dgm.net10.r1$am, 3, colSums)) == 0] = NA
df.10r2[t(apply(dgm.net10.r2$am, 3, colSums)) == 0] = NA
df.10r3[t(apply(dgm.net10.r3$am, 3, colSums)) == 0] = NA
df.10r4[t(apply(dgm.net10.r4$am, 3, colSums)) == 0] = NA
summary(colMeans(df.10r1, na.rm = T))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.6951  0.7418  0.7986  0.8082  0.8766  0.9480 

summary(colMeans(df.10r2, na.rm = T))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.6947  0.7266  0.7995  0.8010  0.8580  0.9554 

summary(colMeans(df.10r3, na.rm = T))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.6649  0.7396  0.8057  0.8083  0.8876  0.9481 

summary(colMeans(df.10r4, na.rm = T))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.6606  0.7258  0.8295  0.8099  0.8806  0.9478 

Load time series data

# f = list.files(file.path(PATH_TS, "node_timeseries", "3T_HCP820_MSMAll_d25_ts2"), "*.txt")
# # read data
# ts = array(NA, dim=c(N_t, Nn, N_runs, N))
# for (s in 1:N) {
#   d = scaleTs(as.matrix(read.table(file.path(PATH_TS, "node_timeseries", "3T_HCP820_MSMAll_d25_ts2", f[s]))))
#   assert(nrow(d) == N_t*N_runs)
#   ts[,,1,s] = d[1:1200,COMP]
#   ts[,,2,s] = d[1201:2400,COMP]
#   ts[,,3,s] = d[2401:3600,COMP]
#   ts[,,4,s] = d[3601:4800,COMP]
# }
# 
# f=file(file.path(PATH,"RSN10-ts.RData"))
# save(ts, file = f)
# close(f)
load(file.path(PATH, 'results', 'RSN10-ts.RData'))

Variance of nodes

We calculate the SD across time for the subject’s node.

# ts 1200 x 10 x 4 x 500
r1=apply(ts[,,1,], c(3,2), sd)
r2=apply(ts[,,2,], c(3,2), sd)
r3=apply(ts[,,3,], c(3,2), sd)
r4=apply(ts[,,4,], c(3,2), sd)

p1 = ggplot(melt(r1), aes(x=as.factor(Var2), y=value)) +
  geom_point(shape=1, color="gray50", size=2) +
  geom_boxplot(width=0.3) +
  ggtitle(sprintf("Run 1", N)) + ylab("SD") + xlab("RSNs")
p2 = ggplot(melt(r2), aes(x=as.factor(Var2), y=value)) +
  geom_point(shape=1, color="gray50", size=2) +
  geom_boxplot(width=0.3) +
  ggtitle(sprintf("Run 2", N)) + ylab("SD") + xlab("RSNs")
p3 = ggplot(melt(r3), aes(x=as.factor(Var2), y=value)) +
  geom_point(shape=1, color="gray50", size=2) +
  geom_boxplot(width=0.3) +
  ggtitle(sprintf("Run 3", N)) + ylab("SD") + xlab("RSNs")
p4 = ggplot(melt(r4), aes(x=as.factor(Var2), y=value)) +
  geom_point(shape=1, color="gray50", size=2) +
  geom_boxplot(width=0.3) +
  ggtitle(sprintf("Run 4", N)) + ylab("SD") + xlab("RSNs")
plot_grid(p1, p2, p3, p4, ncol = 2, nrow = 2)

ggsave(path = PATH_FIG, "Sd_10RSNs.png")

Saving 6 x 4 in image

Correlation Plots

LIM=c(-0.65,0.65)
p1 = gplotMat(rmdiag(corTs(ts[,,1,])), title=sprintf("Run 1", N), lim=LIM, nodeLabels = labels$Label,
              axisTextSize = 8, xAngle = 90, colMapLabel=expression("Pearson\'s"~italic(r)),
              gradient = c("blue", "white", "red")) +
  xlab("Node") + ylab("Node") + scale_x_continuous(breaks = 0.5:9.5, labels = labels$Label )
p2 = gplotMat(rmdiag(corTs(ts[,,2,])), title=sprintf("Run 2", N), lim=LIM, nodeLabels = labels$Label,
              axisTextSize = 8, xAngle = 90, colMapLabel=expression("Pearson\'s"~italic(r)),
              gradient = c("blue", "white", "red")) +
  xlab("Node") + ylab("Node") + scale_x_continuous(breaks = 0.5:9.5, labels = labels$Label )
p3 = gplotMat(rmdiag(corTs(ts[,,3,])), title=sprintf("Run 3", N), lim=LIM, nodeLabels = labels$Label,
              axisTextSize = 8, xAngle = 90, colMapLabel=expression("Pearson\'s"~italic(r)),
              gradient = c("blue", "white", "red")) +
  xlab("Node") + ylab("Node") + scale_x_continuous(breaks = 0.5:9.5, labels = labels$Label )
p4 = gplotMat(rmdiag(corTs(ts[,,4,])), title=sprintf("Run 4", N), lim=LIM, nodeLabels = labels$Label,
              axisTextSize = 8, xAngle = 90, colMapLabel=expression("Pearson\'s"~italic(r)),
              gradient = c("blue", "white", "red")) +
  xlab("Node") + ylab("Node") + scale_x_continuous(breaks = 0.5:9.5, labels = labels$Label )
plot_grid(p1, p2, p3, p4, ncol = 2, nrow = 2, rel_widths = c(1, 1))

ggsave(path = PATH_FIG, "R_10RSNs.png")

Plot example timeseries

idx = 1:round(120/TR) # first 2 minutes
# random sampling of subject s and nodes n
nodes = 5 # no. of nodes to plot
d = ts[idx,sample(Nn,nodes),,sample(N,1)]
p1 = ggplot(melt(d[,1,]), aes(x = Var1, y = value, group=Var2, color=as.factor(Var2))) + 
  geom_line() + theme_minimal() + ggtitle("Run 1")
p2 = ggplot(melt(d[,2,]), aes(x = Var1, y = value, group=Var2, color=as.factor(Var2))) +
  geom_line() + theme_minimal() + ggtitle("Run 2")
p3 = ggplot(melt(d[,3,]), aes(x = Var1, y = value, group=Var2, color=as.factor(Var2))) +
  geom_line() + theme_minimal() + ggtitle("Run 3")
p4 = ggplot(melt(d[,4,]), aes(x = Var1, y = value, group=Var2, color=as.factor(Var2))) + 
  geom_line() + theme_minimal() + ggtitle("Run 4")
plot_grid(p1, p2, p3, p4, ncol = 1, nrow = 4, rel_widths = c(1, 1))

Network consistency across subjects

stats.r1 = binom.nettest(dgm.net10.r1$am, alter = "greater")
stats.r2 = binom.nettest(dgm.net10.r2$am, alter = "greater")
stats.r3 = binom.nettest(dgm.net10.r3$am, alter = "greater")
stats.r4 = binom.nettest(dgm.net10.r4$am, alter = "greater")

Difference between e=0 and e=20

x = c(sum(dgm.net10.r1$tam != dgm.net10.r1$am),
      sum(dgm.net10.r2$tam != dgm.net10.r2$am),
      sum(dgm.net10.r3$tam != dgm.net10.r3$am),
      sum(dgm.net10.r4$tam != dgm.net10.r4$am)
)
print(x)

[1] 300 243 237 229

print(x/(N*Nn*(Nn-1)))

[1] 0.006666667 0.005400000 0.005266667 0.005088889

Figure 8

mylim = c(0.1, 0.68)
pos = 0.5:9.5
s = 0.2
p1 = gplotMat(stats.r1$adj, title = "run 1", titleTextSize = 12, nodeLabels=labels$Label,
              axisTextSize=8, xAngle=90, lim = mylim, barWidth = 0.2, textSize = 11) +
  scale_x_continuous(breaks = pos, labels = labels$Label) +
  theme(plot.margin = unit(c(s, 0, 0, 0), "cm"))
p2 = gplotMat(rmna(stats.r1$adj_fdr), title = "binomial test", titleTextSize = 12, nodeLabels=labels$Label,
              axisTextSize=8, xAngle=90, lim = mylim, barWidth = 0.2, textSize = 11) +
  scale_x_continuous(breaks = pos, labels = labels$Label) +
  theme(plot.margin = unit(c(s, 0, 0, 0), "cm"))
p3 = gplotMat(stats.r2$adj, title = "run 2", titleTextSize = 12, nodeLabels=labels$Label,
              axisTextSize=8, xAngle=90, lim = mylim, barWidth = 0.2, textSize = 11) +
  scale_x_continuous(breaks = pos, labels = labels$Label) +
  theme(plot.margin = unit(c(s, 0, 0, 0), "cm"))
p4 = gplotMat(rmna(stats.r2$adj_fdr), title = "binomial test", titleTextSize = 12, nodeLabels=labels$Label,
              axisTextSize=8, xAngle=90, lim = mylim, barWidth = 0.2, textSize = 11) +
  scale_x_continuous(breaks = pos, labels = labels$Label) +
  theme(plot.margin = unit(c(s, 0, 0, 0), "cm"))
p5 = gplotMat(stats.r3$adj, title = "run 3", titleTextSize = 12, nodeLabels=labels$Label,
              axisTextSize=8, xAngle=90, lim = mylim, barWidth = 0.2, textSize = 11) +
  scale_x_continuous(breaks = pos, labels = labels$Label) +
  theme(plot.margin = unit(c(s, 0, 0, 0), "cm"))
p6 = gplotMat(rmna(stats.r3$adj_fdr), title = "binomial test", titleTextSize = 12, nodeLabels=labels$Label,
              axisTextSize=8, xAngle=90, lim = mylim, barWidth = 0.2, textSize = 11) +
  scale_x_continuous(breaks = pos, labels = labels$Label) +
  theme(plot.margin = unit(c(s, 0, 0, 0), "cm"))
p7 = gplotMat(stats.r4$adj, title = "run 4", titleTextSize = 12, nodeLabels=labels$Label,
              axisTextSize=8, xAngle=90, lim = mylim, barWidth = 0.2, textSize = 11) +
  scale_x_continuous(breaks = pos, labels = labels$Label) +
  theme(plot.margin = unit(c(s, 0, 0, 0), "cm"))
p8 = gplotMat(rmna(stats.r4$adj_fdr), title = "binomial test", titleTextSize = 12, nodeLabels=labels$Label,
              axisTextSize=8, xAngle=90, lim = mylim, barWidth = 0.2, textSize = 11) +
  scale_x_continuous(breaks = pos, labels = labels$Label) +
  theme(plot.margin = unit(c(s, 0, 0, 0), "cm"))
plot_grid(p1, p2, p3, p4, p5, p6, p7, p8, ncol = 2, nrow = 4, rel_heights = c(1,1,1,1))

ggsave(path = PATH_FIG, "Fig8.png")

Networks consistency across runs

# edge reproduced in 3 runs or more
m3 = dgm.net10.r1$am + dgm.net10.r2$am + dgm.net10.r3$am + dgm.net10.r4$am > 2
# edge reproduced in all runs
m4 = dgm.net10.r1$am * dgm.net10.r2$am * dgm.net10.r3$am * dgm.net10.r4$am
# consistently no edge in all runs
mn = (1-dgm.net10.r1$am) * (1-dgm.net10.r2$am) * (1-dgm.net10.r3$am) * (1-dgm.net10.r4$am)
stats.m3 = binom.nettest(m3, alter = "greater")
stats.m4 = binom.nettest(m4, alter = "greater")
stats.mn = binom.nettest(mn, alter = "greater")

Figure 9

p1 = gplotMat(stats.m3$adj, title = "edge in 3/4 runs", nodeLabels=labels$Label,
              axisTextSize=8, xAngle=90, lim = c(0.1, 0.62), barWidth = 0.2, textSize = 11) +
  scale_x_continuous(breaks = pos, labels = labels$Label) +
  theme(plot.margin = unit(c(s, 0, 0, 0), "cm"))
p2 = gplotMat(rmna(stats.m3$adj_fdr), title = "binomial test", nodeLabels=labels$Label,
              axisTextSize=8, xAngle=90, lim = c(0.1, 0.62), barWidth = 0.2, textSize = 11) +
  scale_x_continuous(breaks = pos, labels = labels$Label) +
  theme(plot.margin = unit(c(s, 0, 0, 0), "cm"))
p3 = gplotMat(stats.m4$adj, title = "edge in 4/4 runs", nodeLabels=labels$Label,
              axisTextSize=8, xAngle=90, lim = c(0.1, 0.62), barWidth = 0.2, textSize = 11) +
  scale_x_continuous(breaks = pos, labels = labels$Label) +
  theme(plot.margin = unit(c(s, 0, 0, 0), "cm"))
p4 = gplotMat(rmna(stats.m4$adj_fdr), title = "binomial test", nodeLabels=labels$Label, 
              axisTextSize=8, xAngle=90, lim = c(0.1,0.62), barWidth = 0.2, textSize = 11) +
  scale_x_continuous(breaks = pos, labels = labels$Label) +
  theme(plot.margin = unit(c(s, 0, 0, 0), "cm"))
p5 = gplotMat(rmdiag(stats.mn$adj), title = "no edge in 4/4 runs", nodeLabels=labels$Label,
              gradient = c("white", "violet", "blue"), axisTextSize=8, xAngle=90, 
              barWidth = 0.2, textSize = 11) +
  scale_x_continuous(breaks = pos, labels = labels$Label) +
  theme(plot.margin = unit(c(s, 0, 0, 0), "cm"))
p6 = gplotMat(rmdiag(rmna(stats.mn$adj_fdr)), title = "binomial test", nodeLabels=labels$Label,
              gradient = c("white", "violet", "blue"), axisTextSize=8, xAngle=90,
              barWidth = 0.2, textSize = 11) +
  scale_x_continuous(breaks = pos, labels = labels$Label) +
  theme(plot.margin = unit(c(s, 0, 0, 0), "cm"))
plot_grid(p1, p2, p3, p4, p5, p6, labels=c("A", "", "", "", "B", ""),
          ncol = 2, nrow = 3, rel_heights = c(1,1,1))

ggsave(path = PATH_FIG, "Fig9.png")

Figure connectivity strength

# find representative subject showing group effect [DMN VisMed and Aud]->Cereb
# s=sample(N,1) # subject
s=35
idx=array(FALSE, dim=c(Nn, Nn))
idx[2,8] = idx[3,8] = idx[9,8] = TRUE
#stats.r1$adj_fdr[idx]
a=dgm.net10.r1$am[,,s]
a[idx]

[1] 1 1 1

r=1 # run
node=8
# get model space
sub = read.subject(PATH_NET, sprintf("%s_Run_%03d", SUBJECTS[s], r), Nn)
df = sub$winner[12,node]
pars = sub$winner[2:Nn,node]
pars = pars[pars != 0]
mypars=c(2,9)
idxpar=match(mypars,pars)
Nt=nrow(ts)
Ft=array(1,dim=c(Nt,length(pars)+1))
Ft[,2:ncol(Ft)]=ts[,pars,r,s] # selects parents
Yt=ts[,node,r,s]
# get df corresponding to parent model pars
#df = getModel(sub$models[[1]], pars)[Nn+2]
fit=dlm.lpl(Yt, t(Ft), delta = df)
y = dlm.retro(fit$mt, fit$CSt, fit$RSt, fit$nt, fit$dt)
theta=t(y$smt[2:(length(pars)+1),])
scaleFUN <- function(x) sprintf("%.1f", x)
start=1
seconds=500 # max is 864
idx=round(start/TR):round((start+seconds)/TR)
Nt_=length(idx)
d=melt(Yt[idx])
d$time=1:Nt_
d$labels="Cer  "
p1 = ggplot(d, aes(x = time*TR, y = value, color=labels)) + geom_line() +
  theme_minimal() + ggtitle("Node time series (child)") +
  scale_color_discrete(name = "labels") + xlab("seconds") + ylab("BOLD") + #xlim(c(0,seconds)) +
  scale_y_continuous(labels=scaleFUN) +
  theme(axis.text.x = element_text(size=10), axis.text.y = element_text(size=10),
        axis.title.y = element_text(size=10), legend.text=element_text(size=10),
        plot.title = element_text(size=10,face="bold"))
#p1 = p1 + scale_colour_brewer(palette="Dark2")
#d=melt(Ft[idx,2:(length(pars)+1)])
d=melt(Ft[idx,idxpar+1])
d$time=1:Nt_
d$labels=c(rep("DMN", Nt_),rep("Au", Nt_))
p2 = ggplot(d, aes(x = time*TR, y = value, group=Var2, color=labels)) + geom_line() +
  theme_minimal() + ggtitle("Node time series (parents)") +
  scale_color_discrete(name = "labels") + xlab("seconds") + ylab("BOLD") + #xlim(c(0,seconds)) +
  scale_y_continuous(labels=scaleFUN) +
  theme(axis.text.x = element_text(size=10), axis.text.y = element_text(size=10),
        axis.title.y = element_text(size=10), legend.text=element_text(size=10),
        plot.title = element_text(size=10,face="bold"))
d=melt(theta[idx,idxpar])
d$time=1:Nt_
d$labels=c(rep("DMN", Nt_),rep("Au", Nt_))
p3 = ggplot(d, aes(x = time*TR, y = value, group=Var2, color=labels)) + geom_line() +
  theme_minimal() + ggtitle("Influence of parents on the child node") +
  #scale_color_discrete(name = expression(paste(italic(theta^"r=1"), ""[italic("t")]))) +
  scale_color_discrete(name = "labels") +
  xlab("seconds") + ylab(expression(theta)) + #xlim(c(0,seconds)) +
  theme(axis.text.x = element_text(size=10), axis.text.y = element_text(size=10),
        axis.title.y = element_text(size=10), legend.text=element_text(size=10),
        plot.title = element_text(size=10,face="bold"))
plot_grid(p1, p2, p3, ncol = 1, nrow = 3, rel_widths = c(1,1,1), labels = c("A", "B", "C"))

ggsave(path = PATH_FIG, "Theta_RSN.png")

Saving 5 x 5 in image