DGM-Mouse

Install required packages

# install.packages("devtools")
# install.packages("testit")
# install.packages("ggplot2")
# install.packages("cowplot")
# install.packages("reshape2")

Load libraries

library(DGM)
library(testit)
library(ggplot2)
library(cowplot)
library(reshape2)

Main variables

N=16 # Number of simulated subjects/datasets
Nn=14 # Number of nodes, does not apply to all datasets
Nn1=8
Nn2=6
N_t= 2000 # Volumes/no. of timepoints
PATH_HOME = "/home/simon"
PATH = file.path(PATH_HOME, 'Data', 'DGM-Sim')
PATH_DATA = file.path(PATH_HOME, "Drive", "Mouse")

Loading network data

# subj=list()
# for (s in 1:N) {
#   subj[[s]] = read.subject(file.path(PATH_DATA, 'Mouse_14ROIs_CP_results'), sprintf("Mouse_%03d",s), 14)
# }
# dgm.mouse14=dgm.group(subj)
# 
# subj=list()
# for (s in 1:8) {
#   subj[[s]] = read.subject(file.path(PATH_DATA, 'Mouse_N816ROIs_results'), sprintf("Mouse_%03d",s), 16)
# }
# dgm.mouse16=dgm.group(subj)
# 
# subj=list()
# for (s in 1:N) {
#   subj[[s]] = read.subject(file.path(PATH_DATA, 'Net1_WT_results'), sprintf("Mouse_%03d",s), 8)
# }
# dgm.net1wt=dgm.group(subj)
# 
# subj=list()
# for (s in 1:N) {
#   subj[[s]] = read.subject(file.path(PATH_DATA, 'Net2_WT_results'), sprintf("Mouse_%03d",s), 6)
# }
# dgm.net2wt=dgm.group(subj)
# 
# f=file(file.path(PATH,"results", "DGM-Mouse.RData"))
# save(dgm.mouse14, dgm.mouse16, dgm.net1wt, dgm.net2wt, file =f, compress = T)
# close(f)
# takes quite long to load above data
load(file.path(PATH, "results", "DGM-Mouse.RData"))

Plot: discount factor delta distribution per node

node=as.factor(sort(rep(1:14,N)))
d1=data.frame(df=c(dgm.mouse14$df_), node=node)
node=as.factor(sort(rep(1:16,8)))
d2=data.frame(df=c(dgm.mouse16$df_), node=node)
node=as.factor(sort(rep(1:8,N)))
d3=data.frame(df=c(dgm.net1wt$df_),  node=node)
node=as.factor(sort(rep(1:6,16)))
d4=data.frame(df=c(dgm.net2wt$df_),  node=node)
p1 = ggplot(d1, aes(x=node, y=df)) + geom_boxplot(width=0.4, color="blue1") +
  geom_point(shape=1, color="gray50", size=1, position = position_jitter(width = 0.2, height = 0.0003)) + ggtitle("Mouse 14 ROIs (N=16)")
p2 = ggplot(d2, aes(x=node, y=df)) + geom_boxplot(width=0.4, color="blue1") +
  geom_point(shape=1, color="gray50", size=1, position = position_jitter(width = 0.2, height = 0.0003)) + ggtitle("Mouse 16 ROIs (N=8)")
p3 = ggplot(d3, aes(x=node, y=df)) + geom_boxplot(width=0.4, color="blue1") +
  geom_point(shape=1, color="gray50", size=1, position = position_jitter(width = 0.2, height = 0.0003)) + ggtitle("Net1_WT (N=16)")
p4 = ggplot(d4, aes(x=node, y=df)) + geom_boxplot(width=0.4, color="blue1") +
  geom_point(shape=1, color="gray50", size=1, position = position_jitter(width = 0.2, height = 0.0003)) + ggtitle("Net2_WT (N=16)")
plot_grid(p1, p3, p2, p4, ncol = 2, nrow = 2, rel_widths = c(1, 0.5))

DFs with parents only

df.1 = dgm.net1wt$df_
df.2 = dgm.net2wt$df_
df.1[t(apply(dgm.net1wt$am, 3, colSums)) == 0] = NA
df.2[t(apply(dgm.net2wt$am, 3, colSums)) == 0] = NA
summary(colMeans(df.1, na.rm = T))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
 0.9900  0.9922  0.9930  0.9926  0.9933  0.9940       2 

summary(colMeans(df.2, na.rm = T))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.9550  0.9779  0.9914  0.9835  0.9941  0.9947 

Correlation plot

f14 = list.files(file.path(PATH_DATA, "Timeseries_14ROIs_CP"), "*.txt")
f16 = list.files(file.path(PATH_DATA, "WT_fmr1ko_120d_dr_stage1"), "*.txt")
f1  = list.files(file.path(PATH_DATA, "Net1_WT"), "*.txt")
f2  = list.files(file.path(PATH_DATA, "Net2_WT"), "*.txt")
assert(length(f14) == N)
assert(length(f16) == 8)
assert(length(f1)  == N)
assert(length(f2)  == N)
# read data
d14 = array(NA, dim=c(N_t, Nn, N))
d16 = array(NA, dim=c(N_t, 16, 8))
d1  = array(NA, dim=c(N_t, Nn1, N))
d2  = array(NA, dim=c(N_t, Nn2, N))
for (s in 1:N) {
  d14[,,s]= scaleTs(as.matrix(read.table(file.path(PATH_DATA, "Timeseries_14ROIs_CP", f14[s]))))
  if (s <=8) {
    d16[,,s]= scaleTs(as.matrix(read.table(file.path(PATH_DATA, "WT_fmr1ko_120d_dr_stage1", f16[s]))))
  }
  d1[,,s] = scaleTs(as.matrix(read.table(file.path(PATH_DATA, "Net1_WT", f1[s]))))
  d2[,,s] = scaleTs(as.matrix(read.table(file.path(PATH_DATA, "Net2_WT", f2[s]))))
  #d=scaleTs(d)
  assert(nrow(d14) == N_t)
  assert(nrow(d16) == N_t)
  assert(nrow(d1) == N_t)
  assert(nrow(d2) == N_t)
  assert(ncol(d14) == Nn)
  assert(ncol(d16) == 16)
  assert(ncol(d1) == Nn1)
  assert(ncol(d2) == Nn2)
}
p1 = gplotMat(rmdiag(corTs(d14)), title='14 ROIs (N=16)', lim=c(0, 1),
              colMapLabel=expression("Pearson\'s"~italic(r))) + xlab("Node") + ylab("Node")
p2 = gplotMat(rmdiag(corTs(d16)), title='16 ROIs (N=8)', lim=c(-0.4, 0.6),
              colMapLabel=expression("Pearson\'s"~italic(r)), gradient = c("blue", "white", "red")) + xlab("Node") + ylab("Node")
p3 = gplotMat(rmdiag(corTs(d1)), title='Net1_WT (N=16)', lim=c(0, 0.5),
              colMapLabel=expression("Pearson\'s"~italic(r))) + xlab("Node") + ylab("Node")
p4 = gplotMat(rmdiag(corTs(d2)), title='Net2_WT (N=16)', lim=c(0, 0.5),
              colMapLabel=expression("Pearson\'s"~italic(r))) + xlab("Node") + ylab("Node")
plot_grid(p1, p2, p3, p4, ncol = 2, nrow = 2, rel_widths = c(1, 1))

Plot random timeseries

t = 1:200 # interval to plot
# random sampling mice m and nodes n
nodes = 5 # no. of nodes to plot
d = melt(d14[t,sample(14,nodes),sample(N,1)])
p1=ggplot(d, aes(x = Var1, y = value, group=Var2, color=as.factor(Var2))) + geom_line() + theme_minimal() + 
  ggtitle("Mouse 14 ROIs")
d = melt(d16[t,sample(16,nodes),sample(8,1)])
p2=ggplot(d, aes(x = Var1, y = value, group=Var2, color=as.factor(Var2))) + geom_line() + theme_minimal() + ggtitle("Mouse 16 ROIs")
d = melt(d1[t,sample(8,nodes),sample(N,1)])
p3=ggplot(d, aes(x = Var1, y = value, group=Var2, color=as.factor(Var2))) + geom_line() + theme_minimal() + ggtitle("Net1_WT")
d = melt(d2[t,sample(6,nodes),sample(N,1)])
p4=ggplot(d, aes(x = Var1, y = value, group=Var2, color=as.factor(Var2))) + geom_line() + theme_minimal() + ggtitle("Net2_WT")
plot_grid(p1, p2, p3, p4, ncol = 1, nrow = 4, rel_widths = c(1, 1))

Network consistency across mice

stats1 = binom.nettest(dgm.mouse14$am)
stats2 = binom.nettest(dgm.mouse16$am)
p1 = gplotMat(stats1$adj, title = "Mouse 14 ROIs")
p2 = gplotMat(rmna(stats1$adj_fdr), title = "binom test & FDR")
p3 = gplotMat(stats2$adj, title = "Mouse 16 ROIs")
p4 = gplotMat(rmna(stats2$adj_fdr), title = "binom test & FDR")
plot_grid(p1, p2, p3, p4, ncol = 2, nrow = 2, rel_heights = c(1,1), labels = c("A", "", "B") )

Figure 11

stats3 = binom.nettest(dgm.net1wt$am)
stats4 = binom.nettest(dgm.net2wt$am)
pos = 0.5:8.4
mylabels = c("OrbM", "OrbL", "NAcc", "Cing",  "Amyg", "CA1", "DG", "Entorh")
p5 = gplotMat(stats3$adj, title = "amyg-hipp-enthorin", nodeLabels=mylabels, axisTextSize=9,
              xAngle=90, barWidth = 0.2) + 
  scale_x_continuous(breaks = pos, labels = mylabels)

Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.

p6 = gplotMat(rmna(stats3$adj_fdr), title = "binomial test", nodeLabels=mylabels, axisTextSize=9,
              xAngle=90, barWidth = 0.2) +
  scale_x_continuous(breaks = pos, labels = mylabels)

Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.

pos = 0.5:6.4
mylabels = c("SomSensR", "MotorR", "PutamR", "SomSensL", "MotorL", "PutamL")
p7 = gplotMat(stats4$adj, title = "cort-striat-pallid", nodeLabels=mylabels, axisTextSize=9, 
              xAngle=90, barWidth = 0.2) + 
  scale_x_continuous(breaks = pos, labels = mylabels)

Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.

p8 = gplotMat(rmna(stats4$adj_fdr), title = "binomial test", nodeLabels=mylabels, axisTextSize=9, 
              xAngle=90, barWidth = 0.2) +
  scale_x_continuous(breaks = pos, labels = mylabels)

Scale for 'x' is already present. Adding another scale for 'x', which will replace the existing scale.

plot_grid(p5, p6, p7, p8, ncol = 2, nrow = 2,
          rel_heights = c(1,1), labels = c("A", "", "B") )

ggsave(path = file.path(PATH, 'figures'), "Fig11.png")

Saving 5.2 x 4.5 in image

Difference between e=0 and e=20

x= c(sum(dgm.net1wt$am != dgm.net1wt$tam),
     sum(dgm.net2wt$am != dgm.net2wt$tam))
print(x)

[1] NA NA

print(x/c(N*Nn1*(Nn1-1), N*Nn2*(Nn2-1)))

[1] NA NA

Proportions

stats3$adj

       [,1]   [,2]   [,3]   [,4]   [,5] [,6]   [,7] [,8]
[1,] 0.0000 0.3125 0.3125 0.9375 0.3125    0 0.4375    0
[2,] 0.8125 0.0000 0.3125 0.3125 0.2500    0 0.2500    0
[3,] 0.8125 0.3125 0.0000 0.3750 0.3125    0 0.3750    0
[4,] 0.7500 0.0625 0.1250 0.0000 0.0000    0 0.5000    0
[5,] 0.5000 0.1875 0.2500 0.4375 0.0000    0 0.3125    0
[6,] 0.5000 0.1875 0.3125 0.6250 0.0625    0 0.5000    0
[7,] 0.6250 0.1250 0.1250 0.7500 0.1875    0 0.0000    0
[8,] 0.4375 0.1875 0.0625 0.3125 0.3125    0 0.1250    0

max(stats3$adj)

[1] 0.9375

stats4$adj

       [,1]   [,2]   [,3]   [,4]   [,5]   [,6]
[1,] 0.0000 0.3750 0.1875 0.6875 0.5000 0.5000
[2,] 0.1875 0.0000 0.1875 0.4375 0.9375 0.5625
[3,] 0.2500 0.4375 0.0000 0.3125 0.5000 0.8125
[4,] 0.3125 0.3750 0.1875 0.0000 0.7500 0.3125
[5,] 0.2500 0.6875 0.1875 0.4375 0.0000 0.2500
[6,] 0.2500 0.4375 0.5000 0.2500 0.3125 0.0000

max(stats4$adj)

[1] 0.9375