# Causal sign propagation.
#
# Start the perturbation at +1 and push the sign forward along DIRECTED edges:
# a node's state is its upstream state multiplied by the edge sign. A node
# reached with two conflicting signs is ambiguous (0); a node not reachable from
# the perturbation at all is unknown (NA). Both render grey.
#
# This replaces an earlier version that summed a node's incoming AND outgoing
# edge signs. That is not sign propagation - it is a degree statistic - and in a
# predominantly activating prior network it made almost every node come out +1.
propagate_sign_from_perturbation <- function(edges, root = "PERTURBATION") {
nodes <- sort(unique(c(edges$source, edges$target)))
state <- setNames(rep(NA_real_, length(nodes)), nodes)
if (!root %in% nodes) {
return(tibble(node = nodes, predicted_activity = NA_real_))
}
state[[root]] <- 1
e <- edges %>%
filter(!is.na(sign)) %>%
distinct(source, target, sign) %>%
mutate(sign = sign(sign))
# State only ever moves NA -> +/-1 -> 0, and 0 is absorbing, so this converges
# in at most one pass per node even when the graph contains cycles.
for (iter in seq_along(nodes)) {
changed <- FALSE
for (i in seq_len(nrow(e))) {
upstream <- state[[e$source[[i]]]]
if (is.na(upstream)) next
proposed <- upstream * e$sign[[i]]
current <- state[[e$target[[i]]]]
updated <- if (is.na(current)) proposed else if (current == proposed) current else 0
if (is.na(current) || current != updated) {
state[[e$target[[i]]]] <- updated
changed <- TRUE
}
}
if (!changed) break
}
tibble(node = names(state), predicted_activity = unname(state))
}
upstream_closure <- function(edges, seeds) {
selected <- edges[0, ]; frontier <- unique(seeds); visited_edges <- character()
repeat {
hit <- edges %>% filter(target %in% frontier, !(paste(source, target) %in% visited_edges))
if (!nrow(hit)) break
selected <- bind_rows(selected, hit) %>% distinct()
visited_edges <- unique(c(visited_edges, paste(hit$source, hit$target)))
frontier <- setdiff(unique(hit$source), unique(c(seeds, selected$target)))
if (!length(frontier)) break
}
selected
}
write_network_sidecars <- function(edges, nodes, out_prefix) {
dir.create(dirname(out_prefix), recursive = TRUE, showWarnings = FALSE)
write_tsv(edges %>% transmute(source, interaction = ifelse(sign > 0, "activates", ifelse(sign < 0, "inhibits", "interacts")), target),
paste0(out_prefix, ".sif"), col_names = FALSE)
write_csv(nodes, paste0(out_prefix, "_node_attributes.csv"))
write_csv(edges, paste0(out_prefix, "_edge_attributes.csv"))
dot_nodes <- nodes %>% mutate(shape_dot = case_when(node_type == "TF" ~ "box", node_type == "perturbation" ~ "plaintext", TRUE ~ "ellipse")) %>%
transmute(line = paste0(' "', node, '" [label="', node, '", shape=', shape_dot, ', style=filled, fillcolor="', fill_color, '", color="', border_color, '"];'))
dot_edges <- edges %>% transmute(line = paste0(' "', source, '" -> "', target, '" [color="', edge_color(sign), '"];'))
# Helvetica is one of Graphviz's built-in PostScript font names, so the .dot
# renders identically everywhere; "Arial" resolves only where it is installed.
writeLines(c("digraph causal_network {", " graph [rankdir=LR, fontname=Helvetica];", " node [fontname=Helvetica];", " edge [fontname=Helvetica];", dot_nodes$line, dot_edges$line, "}"), paste0(out_prefix, ".dot"))
}
plot_network <- function(edges, nodes, title) {
if (!nrow(edges) || !nrow(nodes)) return(ggplot() + annotate("text", 0, 0, label = "No network edges after filtering") + theme_void())
graph <- graph_from_data_frame(edges %>% select(source, target), directed = TRUE, vertices = nodes %>% select(node))
# layout_with_fr() is stochastic. Without a fixed seed every render produces new
# coordinates and rewrites every PNG/SVG, which is permanent binary churn in git.
set.seed(NETWORK_LAYOUT_SEED)
coords <- as_tibble(layout_with_fr(graph), .name_repair = "minimal") %>% setNames(c("x", "y")) %>% mutate(node = V(graph)$name) %>% left_join(nodes, by = "node")
edge_xy <- edges %>%
left_join(coords %>% select(node, x_source = x, y_source = y), by = c("source" = "node")) %>%
left_join(coords %>% select(node, x_target = x, y_target = y), by = c("target" = "node"))
ggplot() +
geom_segment(data = edge_xy, aes(x = x_source, y = y_source, xend = x_target, yend = y_target, color = edge_color(sign)), linewidth = 0.7,
arrow = grid::arrow(length = grid::unit(0.12, "inches"), type = "closed"), lineend = "round") +
geom_point(data = coords, aes(x, y, shape = node_type, fill = fill_color, color = border_color), size = 5.5, stroke = 0.9) +
geom_text(data = coords, aes(x, y, label = node), size = 3, family = PLOT_FONT, vjust = -1.05, check_overlap = TRUE) +
scale_shape_manual(values = c(TF = 22, gene = 21, intermediate = 21, perturbation = 24), drop = FALSE) +
scale_fill_identity() + scale_color_identity() +
guides(shape = guide_legend(title = "Node type", override.aes = list(fill = "white", color = "grey30", size = 4))) +
labs(
title = title,
# The fill channel carries two different things, so say which is which.
caption = paste(
"Node fill - target genes: measured log2FC (blue = down, red = up).",
"TFs / intermediates / perturbation: causal sign propagated from PERTURBATION",
"(orange = activated, blue = inhibited, grey = ambiguous or unreachable).",
sep = "\n"
)
) +
coord_equal() + theme_void(base_family = PLOT_FONT) +
theme(plot.title = element_text(face = "bold", hjust = 0.5, size = 12),
plot.caption = element_text(hjust = 0, size = 7, color = "grey30"),
legend.position = "bottom")
}
build_module_network <- function(comparison_name, module_name) {
sig_targets <- deg_all %>% filter(as.character(Comparison) == comparison_name, significant) %>% mutate(target = harmonize_symbol(Gene)) %>% inner_join(modules %>% mutate(target = harmonize_symbol(gene)), by = "target") %>% filter(module == module_name)
if (!nrow(sig_targets)) return(NULL)
active_tfs <- tf_activity %>% filter(as.character(Comparison) == comparison_name, abs(activity_z) >= TF_ACTIVITY_ABS_CUTOFF) %>% pull(tf)
tf_edges <- collectri %>% semi_join(sig_targets, by = "target") %>% filter(tf %in% active_tfs) %>% transmute(source = tf, target, sign = mor, recurrence = 1, edge_class = "TF-target")
if (!nrow(tf_edges)) return(NULL)
upstream_edges <- carnival_edges %>% filter(as.character(Comparison) == comparison_name) %>% select(source, target, sign, recurrence) %>% mutate(edge_class = "CARNIVAL-PKN") %>% upstream_closure(unique(tf_edges$source))
edges <- bind_rows(upstream_edges, tf_edges) %>% distinct(source, target, .keep_all = TRUE)
pred <- propagate_sign_from_perturbation(edges)
rna_here <- deg_all %>% filter(as.character(Comparison) == comparison_name) %>% mutate(node = harmonize_symbol(Gene))
nodes <- tibble(node = sort(unique(c(edges$source, edges$target)))) %>%
mutate(node_type = case_when(node == "PERTURBATION" ~ "perturbation", node %in% unique(tf_edges$source) ~ "TF", node %in% unique(tf_edges$target) ~ "gene", TRUE ~ "intermediate")) %>%
left_join(pred, by = "node") %>% left_join(rna_here %>% select(node, log2FoldChange, padj, significant), by = "node") %>%
mutate(fill_color = ifelse(node_type == "gene" & !is.na(log2FoldChange), rna_color(log2FoldChange), activity_color(predicted_activity)),
border_color = ifelse(SHOW_RNA_EVIDENCE_BORDER & node_type != "gene" & !is.na(log2FoldChange), pal$border_rna, pal$border_default))
if (REMOVE_AMBIGUOUS_NODES_STATIC) {
keep <- nodes %>% filter(node_type %in% c("perturbation", "gene") | (!is.na(predicted_activity) & predicted_activity != 0)) %>% pull(node)
edges <- edges %>% filter(source %in% keep, target %in% keep)
nodes <- nodes %>% filter(node %in% unique(c(edges$source, edges$target)))
}
list(edges = edges, nodes = nodes, targets = unique(sig_targets$Gene), tfs = unique(tf_edges$source))
}
