forked from MattNolanLab/Inter_Intra_Variation
-
Notifications
You must be signed in to change notification settings - Fork 0
/
cond_ind_graph_from_fits.R
146 lines (111 loc) · 4.06 KB
/
cond_ind_graph_from_fits.R
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()