Inter_Intra_Variation/cond_ind_graph_from_fits.R at master · MattNolanLab/Inter_Intra_Variation · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
library(ggraph)
library(igraph)
library(corrplot)


data_summary.mm_vsris <- prep_int_slopes_PCA(data.sc_r, "property", "mm_vsris")

data_intercepts <-  spread(data_summary.mm_vsris[,1:3], property, ind_intercept)
data_slopes <-  spread(data_summary.mm_vsris[,c(1:2, 4)], property, ind_slope)


data.sc_fits <- data_intercepts %>% dplyr::select(vm:fi) %>%
  na.omit


N               <- nrow(data_summary.mm_vsris)
Rho_neurons     <- cor(data.sc_fits)
tol <- 0
Q_neurons       <- corpcor::cor2pcor(Rho_neurons)
Q_neurons       <- Q_neurons %>% as.data.frame

## bootstrap, resampling rows of data matrix
M <- 1000
index <- 1:nrow(data.sc_fits)
## index <- 1:length(unique(id))
## uniqueid <- unique(id)
Q_neurons_star <- array(dim=c(ncol(data.sc_fits), ncol(data.sc_fits), M))
tmp <- NULL
for(i in 1:M){
  index_star <- base::sample(x=index, length(index), replace=TRUE)
  tmp <- rbind(tmp, data_summary.mm_vsris[index_star,])
  ## for(j in 1:length(uniqueid[index_star])){
  ##     tmp <- rbind(tmp, data.sc[which(id == uniqueid[index_star][j]),])
  ## }
  Rho_neurons_star <- cor(data.sc_fits[index_star,])
  Q_neurons_star[,,i] <- as.matrix(corpcor::cor2pcor(Rho_neurons_star))
  tmp <- NULL
}

## bootstrap, resampling animals
if(FALSE){
  M <- 1000
  index <- 1:nrow(data.sc_fits)
  index <- 1:length(unique(id))
  uniqueid <- unique(id)
  Q_neurons_star <- array(dim=c(ncol(data.sc_fits), ncol(data.sc_fits), M))
  tmp <- NULL
  for(i in 1:M){
    index_star <- base::sample(x=index, length(index), replace=TRUE)
    for(j in 1:length(uniqueid[index_star])){
      tmp <- rbind(tmp, data.sc[which(id == uniqueid[index_star][j]),])
    }
    Rho_neurons_star <- cor(data.sc_fits[index_star,])
    Q_neurons_star[,,i] <- as.matrix(corpcor::cor2pcor(Rho_neurons_star))
    tmp <- NULL
  }
}


Q_neurons_low <- apply(Q_neurons_star, c(1,2), quantile, 0.025)
Q_neurons_upp <- apply(Q_neurons_star, c(1,2), quantile, 0.975)
CI <- Q_neurons_low * Q_neurons_upp
CI[CI<0]  <- 0
CI[CI!=0] <- 1
CI <- as.data.frame(CI)
Q_neurons <- CI*Q_neurons
g              <- igraph::graph.adjacency(abs(Q_neurons)>0,
                                          mode="undirected", diag=FALSE)
## g              <- igraph::graph.adjacency(abs(CI)>0,
##                                   mode="undirected", diag=FALSE)
#igraph::V(g)$class     <- c(data.sc_r$property, "dvlocmm")
igraph::V(g)$degree    <- igraph::degree(g)
edge_width    <- NULL
edge_colour      <- NULL
k <- 0
for(i in 1:(ncol(Q_neurons)-1))
{
  for(j in (i+1):ncol(Q_neurons))
    if(abs(Q_neurons[i,j])>0)
    {
      k<-k+1
      edge_width[k]    <- abs(Q_neurons[i,j])/2
      edge_colour[k]   <- ifelse(Q_neurons[i,j] > 0 , "black", "red")
    }
}
igraph::E(g)$width     <- edge_width
igraph::E(g)$colour    <- edge_colour

## graph          <- tidygraph::as_tbl_graph(g)
## ##
## p <- ggraph(graph, 'igraph', algorithm = 'circle') +
##   ## ggraph(graph, layout= "kk") +
##   geom_node_circle(aes( r=0.1), fill="orange", colour="black")+
##   geom_edge_fan(aes(width = edge_width, color= colour),
##                 show.legend = FALSE,
##                 start_cap = circle(1.4, 'cm'),
##                 end_cap = circle(1.4,  'cm'),
##                 lineend = "butt") +
##   ## scale_edge_color_hue(colour=c("black","gray"))+
##   geom_node_text(aes(label=class),
##                  size=6, nudge_x=.0, nudge_y=0.0, colour="blue") +
##   coord_fixed()+
##   theme_graph(foreground = 'steelblue', fg_text_colour = 'white')
## ##

Q <- as.matrix(Q_neurons)
colnames(Q) <- colnames(data.sc_fits)
rownames(Q) <- colnames(data.sc_fits)


## pdf("CIG_neurons.pdf", height=13, width=13)
par(mfrow=c(1,2))
set.seed(111020)
corrplot(Q, type = "upper", order = "FPC")
g$layout <- layout_with_fr
plot(g, vertex.label=colnames(data.sc_fits), vertex.size=20,
     edge.width=20*edge_width, edge.color=edge_colour)
## dev.off()


## pdf("CIG_neurons_circle.pdf", height=13, width=13)
par(mfrow=c(1,2))
set.seed(111020)
corrplot(Q)
g$layout <- layout_in_circle
plot(g, vertex.label=colnames(data.sc_fits), vertex.size=20,
     edge.width=20*edge_width, edge.color=edge_colour)
## dev.off()