-
Notifications
You must be signed in to change notification settings - Fork 0
/
3_plot_dataset_example.R
161 lines (142 loc) · 6.08 KB
/
3_plot_dataset_example.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
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
library(tidyverse)
library(dyngen)
library(dyngen.manuscript)
library(dynutils)
library(tidygraph)
library(ggraph)
exp <- start_analysis("usecase_network_inference")
##################################################
### Load in example ###
##################################################
dataset <- read_rds(exp$dataset_file("bifurcating_loop_1"))
set.seed(1)
feature_info <- dataset$feature_info
feature_network <-
dataset$regulatory_network %>%
arrange(as.character(regulator) == as.character(target))
# add extra edges invisible between regulators from the same module
feature_network_tmp <-
bind_rows(
feature_network,
feature_info %>%
filter(is_tf) %>%
select(module_id, feature_id) %>%
group_by(module_id) %>%
do({
crossing(regulator = .$feature_id, target = .$feature_id) %>%
mutate(effect = -2)
}) %>%
ungroup() %>%
filter(regulator < target)
)
gr <- tbl_graph(nodes = feature_info, edges = feature_network_tmp)
layout <- igraph::layout_with_fr(gr) %>%
dynutils::scale_minmax() %>%
magrittr::set_rownames(feature_info$feature_id) %>%
magrittr::set_colnames(c("x", "y")) %>%
as.data.frame()
cell_wps <- do.call(dynwrap::select_waypoint_cells, c(dataset[c("milestone_ids", "milestone_network", "milestone_percentages", "progressions", "divergence_regions")], list(num_cells_selected = 1)))
cells <- sample(cell_wps, 2) %>% {.[order(match(., dataset$cell_ids))]}
node_df <- bind_cols(feature_info %>% select(-mol_premrna:-mol_protein), layout)
cap <- .015
capend <- .02
edge_df <-
bind_rows(
feature_network %>% mutate(group = "Global GRN"),
dataset$regulatory_network_sc %>%
filter(cell_id %in% cells) %>%
arrange(as.character(regulator) == as.character(target)) %>%
left_join(feature_network %>% select(regulator, target, effect), by = c("regulator", "target")) %>%
mutate(effect = strength, group = paste0("GRN of cell ", match(cell_id, cells)))
) %>%
mutate(
regulator = as.character(regulator),
target = as.character(target),
x_ = layout[regulator, "x"],
y_ = layout[regulator, "y"],
xend_ = layout[target, "x"],
yend_ = layout[target, "y"],
len = sqrt( (x_ - xend_)^2 + (y_ - yend_)^2 ),
al = (len - cap) / len,
alend = (len - capend) / len,
x = al * x_ + (1 - al) * xend_,
y = al * y_ + (1 - al) * yend_,
xend = alend * xend_ + (1 - alend) * x_,
yend = alend * yend_ + (1 - alend) * y_
)
##################################################
### Supp Fig subplot ###
##################################################
arrow_up <- grid::arrow(type = "closed", angle = 30, length = grid::unit(1.3, "mm"))
arrow_down <- grid::arrow(type = "closed", angle = 89, length = grid::unit(1.3, "mm"))
g <- ggplot() +
geom_point(aes(x, y), colour = "gray", node_df) +
geom_segment(aes(x = x, y = y, xend = xend, yend = yend, colour = effect), edge_df %>% filter(effect >= 0), arrow = arrow_up) +
geom_segment(aes(x = x, y = y, xend = xend, yend = yend, colour = effect), edge_df %>% filter(effect < 0), arrow = arrow_down) +
theme_graph(base_family = "Helvetica") +
coord_equal() +
labs(colour = "Regulatory\nactivity") +
scale_colour_gradientn(colours = c(rev(RColorBrewer::brewer.pal(5, "Reds")), RColorBrewer::brewer.pal(5, "Greens")), limits = c(-1, 1)) +
facet_wrap(~group) +
theme(
plot.margin = margin(0, 0, 0, 0, "cm")
)
g
write_rds(g, exp$result("casewise_grn.rds"), compress = "gz")
##################################################
### Figure 2 subplot ###
##################################################
# make plot of ground truth
g1 <- ggplot() +
geom_point(aes(x, y), colour = "gray", node_df) +
geom_segment(aes(x = x, y = y, xend = xend, yend = yend, colour = effect), edge_df %>% filter(effect >= 0, group == "GRN of cell 1"), arrow = arrow_up) +
geom_segment(aes(x = x, y = y, xend = xend, yend = yend, colour = effect), edge_df %>% filter(effect < 0, group == "GRN of cell 1"), arrow = arrow_down) +
theme_graph(base_family = "Helvetica") +
coord_equal() +
labs(colour = "Regulatory\nactivity") +
scale_colour_gradientn(colours = c(rev(RColorBrewer::brewer.pal(5, "Reds")), RColorBrewer::brewer.pal(5, "Greens")), limits = c(-1, 1)) +
theme(
plot.margin = margin(0, 0, 0, 0, "cm")
)
# make plot of prediction
dataset$prior_information$regulators <- dataset$regulators
dataset$prior_information$targets <- dataset$targets
pred4 <- dynwrap::infer_trajectory(dataset, cni_pyscenic_sgbm())
edge_pred_df <-
pred4$regulatory_network_sc %>%
filter(cell_id %in% cells[[1]]) %>%
arrange(desc(abs(strength))) %>%
head(nrow(feature_network)) %>%
arrange(as.character(regulator) == as.character(target)) %>%
mutate(
effect = pmin(strength / quantile(abs(strength), .8), 1),
group = paste0("GRN of cell ", match(cell_id, cells)),
regulator = as.character(regulator),
target = as.character(target),
x_ = layout[regulator, "x"],
y_ = layout[regulator, "y"],
xend_ = layout[target, "x"],
yend_ = layout[target, "y"],
len = sqrt( (x_ - xend_)^2 + (y_ - yend_)^2 ),
al = (len - cap) / len,
alend = (len - capend) / len,
x = al * x_ + (1 - al) * xend_,
y = al * y_ + (1 - al) * yend_,
xend = alend * xend_ + (1 - alend) * x_,
yend = alend * yend_ + (1 - alend) * y_
)
g2 <- ggplot() +
geom_point(aes(x, y), colour = "gray", node_df) +
geom_segment(aes(x = x, y = y, xend = xend, yend = yend, colour = effect), edge_pred_df %>% filter(effect >= 0), arrow = arrow_up) +
geom_segment(aes(x = x, y = y, xend = xend, yend = yend, colour = effect), edge_pred_df %>% filter(effect < 0), arrow = arrow_down) +
theme_graph(base_family = "Helvetica") +
coord_equal() +
labs(colour = "Regulatory\nactivity") +
scale_colour_gradientn(
colours = c(rev(RColorBrewer::brewer.pal(5, "Reds")), RColorBrewer::brewer.pal(5, "Greens")),
limits = c(-1, 1)
) +
theme(
plot.margin = margin(0, 0, 0, 0, "cm")
)
write_rds(list(groundtruth = g1, prediction = g2), exp$result("usecase_separateplots.rds"), compress = "gz")