WGCNA Module Perturbation Test
F. V. Rodriguez
2026-04-10
Overview
This notebook tests whether the intramodular connectivity loss observed in genotype-response modules (Inv4m DEGs) exceeds what random genotype label shuffling produces. The topology control uses leaf-gradient genes (photosynthesis + senescence hub genes from Ojeda 2026) that have real coexpression structure but are not differentially expressed between genotypes.
Both module sets come from separate consensus WGCNA pipeline runs (1000 iterations each), with bootstrap-validated stability (LFP from step 05).
Gene set terminology:
- Genotype-response genes – differentially expressed between CTRL and Inv4m plants; coexpression driven by the Inv4m introgression.
- Leaf-gradient genes – photosynthesis and senescence hub genes (Ojeda 2026); coexpression driven by the leaf developmental gradient, independent of Inv4m genotype.
Prediction: Genotype-response modules show significant kWithin loss; leaf-gradient modules do not.
Setup
library(here)
library(tidyverse)
library(WGCNA)
library(ggridges)
library(patchwork)
allowWGCNAThreads()## Allowing multi-threading with up to 10 threads.set.seed(42)N_PERM <- 10000
# Genotype-response run
GR_RUN <- here("results/inversion_paper/field_perturbation/run_20251231_201332")
# Leaf-gradient run
LG_RUN <- here("results/inversion_paper/field_perturbation/run_leaf_gradient")
source(here("scripts/utils/setup_paths.R"))
paths <- setup_project_paths("inversion_paper")Load Pre-computed Data
# Genotype-response consensus modules (from CTRL reference network)
gr_modules <- read_csv(
file.path(GR_RUN, "04_consensus_networks", "ctrl", "consensus_modules.csv"),
show_col_types = FALSE
)
gr_lfp <- read_csv(
file.path(GR_RUN, "05_bootstrap_support", "ctrl", "module_support.csv"),
show_col_types = FALSE
)
gr_preservation <- read_csv(
file.path(GR_RUN, "06_preservation", "zsummary_table.csv"),
show_col_types = FALSE
)
# Expression and metadata (shared across runs)
expr_gr <- readRDS(file.path(GR_RUN, "02_gene_filter", "filtered_expression.rds"))
sample_meta <- read_csv(
file.path(GR_RUN, "01_data_prep", "sample_metadata.csv"),
show_col_types = FALSE
)
SOFT_POWER <- readRDS(
file.path(GR_RUN, "03_reference_network", "ctrl", "selected_power.rds")
)
# GO annotations
go_bp_file <- file.path(GR_RUN, "07_module_annotation", "top_BP_per_module_rrvgo.csv")
gr_go_terms <- if (file.exists(go_bp_file)) {
go_bp <- read_csv(go_bp_file, show_col_types = FALSE)
setNames(str_trunc(go_bp$Description, 35), go_bp$module_color)
} else {
NULL
}
cat("Genotype-response:", nrow(gr_modules), "genes,",
n_distinct(gr_modules$color[gr_modules$color != "grey"]), "modules\n")## Genotype-response: 465 genes, 13 modulescat("Soft power:", SOFT_POWER, "\n")## Soft power: 18# Leaf-gradient consensus modules (from CTRL reference network)
lg_modules <- read_csv(
file.path(LG_RUN, "04_consensus_networks", "ctrl", "consensus_modules.csv"),
show_col_types = FALSE
)
lg_lfp <- read_csv(
file.path(LG_RUN, "05_bootstrap_support", "ctrl", "module_support.csv"),
show_col_types = FALSE
)
lg_preservation <- read_csv(
file.path(LG_RUN, "06_preservation", "zsummary_table.csv"),
show_col_types = FALSE
)
# Leaf-gradient expression matrix
expr_lg <- readRDS(file.path(LG_RUN, "02_gene_filter", "filtered_expression.rds"))
# Remap leaf-gradient colors to avoid collisions with genotype-response
lg_color_remap <- c(
turquoise = "sky",
blue = "peach",
brown = "jade",
yellow = "coral",
green = "lavender",
grey = "grey"
)
lg_color_hex <- c(
sky = "#5DA5DA",
peach = "#FAA43A",
jade = "#60BD68",
coral = "#F15854",
lavender = "#B276B2",
grey = "grey"
)
cat("Leaf-gradient:", nrow(lg_modules), "genes,",
n_distinct(lg_modules$color[lg_modules$color != "grey"]), "modules\n")## Leaf-gradient: 218 genes, 5 modulesGO Enrichment for Leaf-gradient Modules
Step 07 was skipped for the leaf-gradient run (DEG-specific inputs). Run GO BP enrichment here using clusterProfiler + rrvgo, matching the approach used for genotype-response modules.
library(clusterProfiler)
library(rrvgo)
library(AnnotationHub)
TERM2GENE <- readRDS(here("data", "TERM2GENE.rds"))
TERM2NAME <- readRDS(here("data", "TERM2NAME.rds"))
# Universe = all genes in the full expression matrix
full_expr <- readRDS(file.path(GR_RUN, "01_data_prep", "combined_expression.rds"))
if (inherits(full_expr, "EList")) {
go_universe <- rownames(full_expr$E)
} else {
go_universe <- rownames(as.matrix(full_expr))
}
rm(full_expr)
# Gene lists per leaf-gradient module (use remapped colors, exclude grey)
lg_remapped <- lg_color_remap[lg_modules$color]
names(lg_remapped) <- lg_modules$gene
lg_gene_lists <- split(
names(lg_remapped[lg_remapped != "grey"]),
lg_remapped[lg_remapped != "grey"]
)
cat("Leaf-gradient modules for GO:", length(lg_gene_lists), "\n")## Leaf-gradient modules for GO: 5cat("Module sizes:", paste(names(lg_gene_lists),
sapply(lg_gene_lists, length), sep = "=", collapse = ", "), "\n")## Module sizes: coral=26, jade=30, lavender=13, peach=64, sky=78lg_ego <- function(x) {
enricher(
gene = x,
pvalueCutoff = 0.05,
pAdjustMethod = "none",
universe = go_universe,
minGSSize = 10,
maxGSSize = 500,
qvalueCutoff = 0.2,
TERM2GENE = TERM2GENE[TERM2GENE$GO %in% TERM2NAME[["BP"]]$go_id, ],
TERM2NAME = TERM2NAME[["BP"]]
)
}
lg_compGO <- compareCluster(
geneCluster = lg_gene_lists,
fun = lg_ego
)
lg_go_results <- as.data.frame(lg_compGO@compareClusterResult) |>
dplyr::rename(module_color = Cluster) |>
arrange(module_color, pvalue)
cat("Total enriched BP terms:", nrow(lg_go_results), "\n")## Total enriched BP terms: 365cat("Terms per module:\n")## Terms per module:print(table(lg_go_results$module_color))##
## coral jade lavender peach sky
## 48 55 78 122 62# rrvgo semantic reduction
ah <- AnnotationHub()
org.Zm.eg.db <- ah[["AH119718"]]
lg_go_filtered <- lg_go_results |>
filter(Count > 1) |>
group_by(module_color) |>
arrange(pvalue, .by_group = TRUE) |>
ungroup()
lg_go_ids <- unique(lg_go_results$ID)
if (length(lg_go_ids) > 1) {
simMatrix <- calculateSimMatrix(
lg_go_ids,
orgdb = org.Zm.eg.db,
ont = "BP",
method = "Rel"
)
scores <- setNames(-log10(lg_go_results$pvalue), lg_go_results$ID)
reducedTerms <- reduceSimMatrix(
simMatrix, scores, threshold = 0.7, orgdb = org.Zm.eg.db
)
lg_go_rrvgo <- lg_go_filtered |>
filter(ID %in% reducedTerms$parent)
} else {
lg_go_rrvgo <- lg_go_filtered
}
# Top 1 BP term per leaf-gradient module
lg_top_bp <- lg_go_rrvgo |>
group_by(module_color) |>
slice_min(pvalue, n = 1, with_ties = FALSE) |>
ungroup() |>
dplyr::select(module_color, ID, Description, pvalue, Count)
lg_go_terms <- setNames(
str_trunc(lg_top_bp$Description, 35),
lg_top_bp$module_color
)
cat("\nTop BP term per leaf-gradient module:\n")##
## Top BP term per leaf-gradient module:print(lg_top_bp)## # A tibble: 5 × 5
## module_color ID Description pvalue Count
## <fct> <chr> <chr> <dbl> <int>
## 1 coral GO:0016042 lipid catabolic process 8.63e- 4 3
## 2 jade GO:0046395 carboxylic acid catabolic process 7.17e- 5 4
## 3 lavender GO:0015979 photosynthesis 2.77e-20 12
## 4 peach GO:0015979 photosynthesis 5.06e-40 31
## 5 sky GO:0015979 photosynthesis 6.09e-60 44Prepare Expression Matrices
# Module color vectors
gr_colors <- setNames(gr_modules$color, gr_modules$gene)
gr_colors <- gr_colors[rownames(expr_gr)]
lg_colors <- lg_color_remap[lg_modules$color]
names(lg_colors) <- lg_modules$gene
lg_colors <- lg_colors[rownames(expr_lg)]
# Genotype labels
true_labels <- setNames(sample_meta$Genotype, sample_meta$sample_id)
all_samples <- sample_meta$sample_id
n_ctrl <- sum(true_labels == "CTRL")
n_inv4m <- sum(true_labels == "Inv4m")
n_total <- length(true_labels)
# WGCNA format: samples x genes
datExpr_gr <- t(expr_gr[, all_samples])
datExpr_lg <- t(expr_lg[, all_samples])
cat("Genotype-response matrix:", nrow(datExpr_gr), "samples x",
ncol(datExpr_gr), "genes\n")## Genotype-response matrix: 43 samples x 465 genescat("Leaf-gradient matrix:", nrow(datExpr_lg), "samples x",
ncol(datExpr_lg), "genes\n")## Leaf-gradient matrix: 43 samples x 218 genescat("CTRL:", n_ctrl, "| Inv4m:", n_inv4m, "\n")## CTRL: 20 | Inv4m: 23Supplementary Figure S5: Bootstrap Support
Two-panel figure showing CTRL consensus module bootstrap support (LFP) for both gene sets.
STABILITY_STABLE <- 0.70
STABILITY_MODERATE <- 0.50
# Combine both gene sets into one data frame
gr_support_df <- gr_lfp |>
filter(color != "grey") |>
mutate(
display_color = color,
fill_hex = color,
gene_set = "Genotype Response DEGs"
)
lg_support_df <- lg_lfp |>
filter(color != "grey") |>
mutate(
display_color = lg_color_remap[color],
fill_hex = lg_color_hex[lg_color_remap[color]],
gene_set = "Leaf Gradient Hubs"
)
support_df <- bind_rows(gr_support_df, lg_support_df) |>
mutate(gene_set = recode(gene_set,
"Genotype Response DEGs" = "Genotype Response",
"Leaf Gradient Hubs" = "Leaf Gradient"))
p_support <- ggplot(
support_df,
aes(x = reorder(display_color, mean_lfp), y = mean_lfp, fill = fill_hex)
) +
geom_col(color = "black", alpha = 0.7) +
geom_errorbar(
aes(ymin = pmax(0, mean_lfp - se_lfp),
ymax = pmin(1, mean_lfp + se_lfp)),
width = 0.3
) +
geom_hline(yintercept = STABILITY_STABLE, linetype = "dashed",
color = "darkgreen") +
geom_hline(yintercept = STABILITY_MODERATE, linetype = "dashed",
color = "orange") +
scale_fill_identity() +
coord_flip() +
facet_grid(gene_set ~ ., scales = "free_y", space = "free_y") +
theme_minimal(base_size = 14) +
theme(panel.spacing.y = unit(1, "lines"),
strip.text.y = element_text(size = 14)) +
labs(x = "Module", y = "Mean Support (LFP)")
print(p_support)
ggsave(file.path(paths$figures, "module_support.pdf"),
p_support, width = 10, height = 5.5)
ggsave(file.path(paths$figures, "module_support.png"),
p_support, width = 10, height = 8, dpi = 300)Helper Function
compute_module_delta_k <- function(datExpr_a, datExpr_b, power, mod_colors) {
adj_a <- adjacency(datExpr_a, power = power, type = "signed")
adj_b <- adjacency(datExpr_b, power = power, type = "signed")
k_a <- intramodularConnectivity(adj_a, mod_colors)
k_b <- intramodularConnectivity(adj_b, mod_colors)
delta_kWithin <- k_b$kWithin - k_a$kWithin
delta_kOut <- k_b$kOut - k_a$kOut
module_ids <- unique(mod_colors[mod_colors != "grey" & !is.na(mod_colors)])
list(
kWithin = sapply(module_ids, function(m) {
idx <- which(mod_colors == m)
mean(delta_kWithin[idx], na.rm = TRUE)
}),
kOut = sapply(module_ids, function(m) {
idx <- which(mod_colors == m)
mean(delta_kOut[idx], na.rm = TRUE)
})
)
}Observed Delta kWithin
ctrl_idx <- which(true_labels == "CTRL")
inv4m_idx <- which(true_labels == "Inv4m")
obs_gr <- compute_module_delta_k(
datExpr_gr[ctrl_idx, ], datExpr_gr[inv4m_idx, ],
SOFT_POWER, gr_colors
)
observed_gr_kw <- obs_gr$kWithin
observed_gr_out <- obs_gr$kOut
obs_lg <- compute_module_delta_k(
datExpr_lg[ctrl_idx, ], datExpr_lg[inv4m_idx, ],
SOFT_POWER, lg_colors
)
observed_lg_kw <- obs_lg$kWithin
observed_lg_out <- obs_lg$kOut
cat("Observed mean delta kWithin per genotype-response module:\n")## Observed mean delta kWithin per genotype-response module:print(round(sort(observed_gr_kw), 4))## tan brown blue yellow pink purple
## -3.0065 -2.3642 -1.6366 -1.4675 -1.3624 -1.1055
## black greenyellow green red salmon turquoise
## -0.8012 -0.6313 -0.6115 -0.5079 -0.3888 -0.2768
## magenta
## -0.1432cat("\nObserved mean delta kWithin per leaf-gradient module:\n")##
## Observed mean delta kWithin per leaf-gradient module:print(round(sort(observed_lg_kw), 4))## sky jade coral peach lavender
## -13.8198 -4.3914 -2.3498 -2.0876 -1.3633cat("\nObserved mean delta kOut per genotype-response module:\n")##
## Observed mean delta kOut per genotype-response module:print(round(sort(observed_gr_out), 4))## pink purple yellow greenyellow brown black
## -1.5068 -1.3451 -0.7252 -0.4355 -0.3494 0.2222
## blue magenta red tan turquoise green
## 0.2779 0.7208 0.8851 1.0936 1.1145 1.5777
## salmon
## 1.6231cat("\nObserved mean delta kOut per leaf-gradient module:\n")##
## Observed mean delta kOut per leaf-gradient module:print(round(sort(observed_lg_out), 4))## lavender sky peach coral jade
## -8.9763 -4.0923 -3.1630 -1.6293 -1.4351Genotype Label Permutation Test
Shuffle genotype labels 10,000 times, keeping group sizes fixed (20 CTRL, 23 Inv4m). For each permutation, recompute \(\Delta k_{Within}\) and \(\Delta k_{Out}\) for both gene sets. Module definitions remain fixed throughout.
null_gr_kw <- matrix(NA, nrow = N_PERM, ncol = length(observed_gr_kw))
null_gr_out <- matrix(NA, nrow = N_PERM, ncol = length(observed_gr_out))
colnames(null_gr_kw) <- names(observed_gr_kw)
colnames(null_gr_out) <- names(observed_gr_out)
null_lg_kw <- matrix(NA, nrow = N_PERM, ncol = length(observed_lg_kw))
null_lg_out <- matrix(NA, nrow = N_PERM, ncol = length(observed_lg_out))
colnames(null_lg_kw) <- names(observed_lg_kw)
colnames(null_lg_out) <- names(observed_lg_out)
for (i in seq_len(N_PERM)) {
perm_idx <- sample(n_total)
pseudo_ctrl <- perm_idx[1:n_ctrl]
pseudo_inv4m <- perm_idx[(n_ctrl + 1):n_total]
perm_gr <- compute_module_delta_k(
datExpr_gr[pseudo_ctrl, ], datExpr_gr[pseudo_inv4m, ],
SOFT_POWER, gr_colors
)
null_gr_kw[i, ] <- perm_gr$kWithin
null_gr_out[i, ] <- perm_gr$kOut
perm_lg <- compute_module_delta_k(
datExpr_lg[pseudo_ctrl, ], datExpr_lg[pseudo_inv4m, ],
SOFT_POWER, lg_colors
)
null_lg_kw[i, ] <- perm_lg$kWithin
null_lg_out[i, ] <- perm_lg$kOut
}Results
ECDF P-values
One-tailed test: fraction of permuted values \(\leq\) observed (testing for connectivity loss more extreme than the null).
# kWithin
pvals_gr_kw <- sapply(names(observed_gr_kw), function(m) {
mean(null_gr_kw[, m] <= observed_gr_kw[m])
})
pvals_lg_kw <- sapply(names(observed_lg_kw), function(m) {
mean(null_lg_kw[, m] <= observed_lg_kw[m])
})
# kOut (one-tailed: test for gain in inter-modular connectivity)
pvals_gr_out <- sapply(names(observed_gr_out), function(m) {
mean(null_gr_out[, m] >= observed_gr_out[m])
})
pvals_lg_out <- sapply(names(observed_lg_out), function(m) {
mean(null_lg_out[, m] >= observed_lg_out[m])
})Genotype-response Modules
results_gr <- tibble(
module = names(observed_gr_kw),
gene_set = "Genotype Response DEGs",
n_genes = sapply(names(observed_gr_kw), function(m)
sum(gr_colors == m, na.rm = TRUE)),
observed_mean_delta_k = round(observed_gr_kw, 4),
null_mean = round(colMeans(null_gr_kw), 4),
null_sd = round(apply(null_gr_kw, 2, sd), 4),
ecdf_pvalue = round(pvals_gr_kw, 4)
) |>
arrange(ecdf_pvalue)
results_gr## # A tibble: 13 × 7
## module gene_set n_genes observed_mean_delta_k null_mean null_sd ecdf_pvalue
## <chr> <chr> <int> <dbl> <dbl> <dbl> <dbl>
## 1 brown Genotyp… 46 -2.36 -0.0912 0.608 0
## 2 tan Genotyp… 12 -3.01 -0.0518 0.974 0.0011
## 3 purple Genotyp… 22 -1.11 -0.0426 0.433 0.0056
## 4 blue Genotyp… 51 -1.64 -0.104 0.647 0.0078
## 5 pink Genotyp… 30 -1.36 -0.0598 0.572 0.0094
## 6 yellow Genotyp… 36 -1.47 -0.0782 0.737 0.0282
## 7 salmon Genotyp… 11 -0.389 -0.0167 0.275 0.0936
## 8 greenye… Genotyp… 16 -0.631 -0.0394 0.469 0.106
## 9 red Genotyp… 34 -0.508 -0.0643 0.391 0.128
## 10 black Genotyp… 32 -0.801 -0.0751 0.648 0.134
## 11 green Genotyp… 36 -0.612 -0.0741 0.652 0.205
## 12 magenta Genotyp… 24 -0.143 -0.0466 0.269 0.361
## 13 turquoi… Genotyp… 55 -0.277 -0.135 1.04 0.441Leaf-gradient Modules
results_lg <- tibble(
module = names(observed_lg_kw),
gene_set = "Leaf Gradient Hubs",
n_genes = sapply(names(observed_lg_kw), function(m)
sum(lg_colors == m)),
observed_mean_delta_k = round(observed_lg_kw, 4),
null_mean = round(colMeans(null_lg_kw), 4),
null_sd = round(apply(null_lg_kw, 2, sd), 4),
ecdf_pvalue = round(pvals_lg_kw, 4)
) |>
arrange(ecdf_pvalue)
results_lg## # A tibble: 5 × 7
## module gene_set n_genes observed_mean_delta_k null_mean null_sd ecdf_pvalue
## <chr> <chr> <int> <dbl> <dbl> <dbl> <dbl>
## 1 coral Leaf Gra… 26 -2.35 -0.116 1.83 0.117
## 2 jade Leaf Gra… 30 -4.39 -0.129 3.73 0.164
## 3 sky Leaf Gra… 78 -13.8 -0.299 12.4 0.190
## 4 lavender Leaf Gra… 13 -1.36 -0.0403 1.83 0.283
## 5 peach Leaf Gra… 64 -2.09 -0.259 4.53 0.344BH FDR Correction Across All Tests
all_results <- bind_rows(results_gr, results_lg)
# Joint FDR across both gene sets
n_gr <- length(pvals_gr_kw)
n_lg <- length(pvals_lg_kw)
all_kw_pvals <- c(pvals_gr_kw, pvals_lg_kw)
all_kw_fdr <- p.adjust(all_kw_pvals, method = "BH")
kw_fdr_gr <- setNames(all_kw_fdr[1:n_gr], names(pvals_gr_kw))
kw_fdr_lg <- setNames(all_kw_fdr[(n_gr + 1):(n_gr + n_lg)], names(pvals_lg_kw))
# kOut FDR
n_gr_out <- length(pvals_gr_out)
n_lg_out <- length(pvals_lg_out)
all_kout_pvals <- c(pvals_gr_out, pvals_lg_out)
all_kout_fdr <- p.adjust(all_kout_pvals, method = "BH")
kout_fdr_gr <- setNames(all_kout_fdr[1:n_gr_out], names(pvals_gr_out))
kout_fdr_lg <- setNames(all_kout_fdr[(n_gr_out + 1):(n_gr_out + n_lg_out)],
names(pvals_lg_out))
all_results <- all_results |>
mutate(
fdr = c(kw_fdr_gr[results_gr$module], kw_fdr_lg[results_lg$module]),
sig = ifelse(fdr < 0.05, "*", "")
) |>
arrange(ecdf_pvalue)
cat("Significant at FDR < 0.05:\n")## Significant at FDR < 0.05:all_results |>
filter(fdr < 0.05) |>
dplyr::select(module, gene_set, n_genes, observed_mean_delta_k,
ecdf_pvalue, fdr) |>
print()## # A tibble: 5 × 6
## module gene_set n_genes observed_mean_delta_k ecdf_pvalue fdr
## <chr> <chr> <int> <dbl> <dbl> <dbl>
## 1 brown Genotype Response DEGs 46 -2.36 0 0
## 2 tan Genotype Response DEGs 12 -3.01 0.0011 0.0099
## 3 purple Genotype Response DEGs 22 -1.11 0.0056 0.0336
## 4 blue Genotype Response DEGs 51 -1.64 0.0078 0.0338
## 5 pink Genotype Response DEGs 30 -1.36 0.0094 0.0338cat("\nAll leaf-gradient modules:\n")##
## All leaf-gradient modules:all_results |>
filter(gene_set == "Leaf Gradient Hubs") |>
dplyr::select(module, n_genes, observed_mean_delta_k, ecdf_pvalue, fdr) |>
print()## # A tibble: 5 × 5
## module n_genes observed_mean_delta_k ecdf_pvalue fdr
## <chr> <int> <dbl> <dbl> <dbl>
## 1 coral 26 -2.35 0.117 0.220
## 2 jade 30 -4.39 0.164 0.246
## 3 sky 78 -13.8 0.190 0.262
## 4 lavender 13 -1.36 0.283 0.340
## 5 peach 64 -2.09 0.344 0.383Save Permutation Results
OUTPUT_DIR <- file.path(GR_RUN, "06_preservation")
write_csv(all_results,
file.path(OUTPUT_DIR,
"permutation_delta_kWithin_results.csv"))
saveRDS(
list(
observed_gr_kw = observed_gr_kw,
observed_lg_kw = observed_lg_kw,
observed_gr_out = observed_gr_out,
observed_lg_out = observed_lg_out,
null_gr_kw = null_gr_kw,
null_lg_kw = null_lg_kw,
null_gr_out = null_gr_out,
null_lg_out = null_lg_out,
pvals_gr_kw = pvals_gr_kw,
pvals_lg_kw = pvals_lg_kw,
pvals_gr_out = pvals_gr_out,
pvals_lg_out = pvals_lg_out,
gr_colors = gr_colors,
lg_colors = lg_colors
),
file.path(OUTPUT_DIR, "permutation_delta_kWithin.rds")
)Visualization
Z-score Ridge Plots (kWithin)
Two-panel figure comparing connectivity disruption in genotype-response modules versus leaf-gradient modules. Z-scores are computed per module as \((observed - mean_{null}) / sd_{null}\). Red lines mark observed values.
# Z-score genotype-response null distributions
gr_null_mean <- colMeans(null_gr_kw)
gr_null_sd <- apply(null_gr_kw, 2, sd)
gr_null_z <- null_gr_kw
for (m in colnames(gr_null_z)) {
gr_null_z[, m] <- (gr_null_z[, m] - gr_null_mean[m]) / gr_null_sd[m]
}
gr_obs_z <- (observed_gr_kw - gr_null_mean) / gr_null_sd
# Z-score leaf-gradient null distributions
lg_null_mean <- colMeans(null_lg_kw)
lg_null_sd <- apply(null_lg_kw, 2, sd)
lg_null_z <- null_lg_kw
for (m in colnames(lg_null_z)) {
lg_null_z[, m] <- (lg_null_z[, m] - lg_null_mean[m]) / lg_null_sd[m]
}
lg_obs_z <- (observed_lg_kw - lg_null_mean) / lg_null_sd
# Module sizes for y-axis labels
gr_mod_sizes <- sapply(colnames(gr_null_z), function(m)
sum(gr_colors == m, na.rm = TRUE))
lg_mod_sizes <- sapply(colnames(lg_null_z), function(m)
sum(lg_colors == m))
# Build combined density data
gr_density <- map_dfr(colnames(gr_null_z), function(m) {
tibble(module = sprintf("%s\nn=%d", m, gr_mod_sizes[m]),
fill_hex = m, z_value = gr_null_z[, m],
gene_set = "Genotype Response DEGs")
})
lg_density <- map_dfr(colnames(lg_null_z), function(m) {
tibble(module = sprintf("%s\nn=%d", m, lg_mod_sizes[m]),
fill_hex = lg_color_hex[m], z_value = lg_null_z[, m],
gene_set = "Leaf Gradient Hubs")
})
density_kw <- bind_rows(gr_density, lg_density)
# Build combined points data
gr_points <- tibble(
module = sprintf("%s\nn=%d", names(gr_obs_z), gr_mod_sizes[names(gr_obs_z)]),
fill_hex = names(gr_obs_z),
z_observed = gr_obs_z, fdr = kw_fdr_gr[names(gr_obs_z)],
gene_set = "Genotype Response DEGs"
)
lg_points <- tibble(
module = sprintf("%s\nn=%d", names(lg_obs_z), lg_mod_sizes[names(lg_obs_z)]),
fill_hex = lg_color_hex[names(lg_obs_z)],
z_observed = lg_obs_z, fdr = kw_fdr_lg[names(lg_obs_z)],
gene_set = "Leaf Gradient Hubs"
)
points_kw <- bind_rows(gr_points, lg_points)
# Order modules by z_observed WITHIN each gene_set, compute per-facet rank
points_kw <- points_kw |>
group_by(gene_set) |>
arrange(z_observed, .by_group = TRUE) |>
mutate(facet_rank = row_number()) |>
ungroup()
# Set factor levels per facet ordering
mod_order <- points_kw |> arrange(gene_set, facet_rank) |> pull(module)
density_kw$module <- factor(density_kw$module, levels = mod_order)
points_kw$module <- factor(points_kw$module, levels = mod_order)
# Propagate facet_rank to density data
density_kw <- density_kw |>
left_join(points_kw |> dplyr::select(module, gene_set, facet_rank) |> distinct(),
by = c("module", "gene_set"))
# Build fill color map (module label -> hex)
fill_map <- setNames(points_kw$fill_hex, as.character(points_kw$module))
# GO term annotations
gr_go_df <- if (!is.null(gr_go_terms)) {
tibble(
module = sprintf("%s\nn=%d", names(gr_go_terms),
gr_mod_sizes[names(gr_go_terms)]),
go_desc = unname(gr_go_terms), gene_set = "Genotype Response DEGs"
) |> filter(module %in% levels(density_kw$module))
}
lg_go_df <- if (!is.null(lg_go_terms)) {
tibble(
module = sprintf("%s\nn=%d", names(lg_go_terms),
lg_mod_sizes[names(lg_go_terms)]),
go_desc = unname(lg_go_terms), gene_set = "Leaf Gradient Hubs"
) |> filter(module %in% levels(density_kw$module))
}
go_df <- bind_rows(gr_go_df, lg_go_df) |>
mutate(module = factor(module, levels = levels(density_kw$module))) |>
left_join(points_kw |> dplyr::select(module, gene_set, facet_rank) |> distinct(),
by = c("module", "gene_set")) |>
mutate(y_pos = facet_rank + 0.35)
p_kw <- ggplot(density_kw, aes(x = z_value, y = module, fill = module)) +
geom_density_ridges(color = "grey40", alpha = 0.7,
scale = 0.9, rel_min_height = 0.01) +
scale_fill_manual(values = fill_map) +
guides(fill = "none") +
geom_vline(xintercept = 0, linetype = "dashed", color = "grey50") +
geom_segment(
data = points_kw,
aes(x = z_observed, xend = z_observed,
y = facet_rank, yend = facet_rank + 0.8),
color = "red", linewidth = 0.8, inherit.aes = FALSE
) +
geom_text(
data = points_kw |> filter(fdr < 0.05),
aes(x = z_observed - 0.6, y = facet_rank + 0.35,
label = sprintf("%.3f", fdr)),
color = "red", size = 3, hjust = 0.5, inherit.aes = FALSE
) +
geom_text(
data = go_df,
aes(x = 2, y = y_pos, label = go_desc),
hjust = 0, size = 14 * 0.8 / .pt,
inherit.aes = FALSE
) +
scale_x_continuous(limits = c(NA, 5), breaks = seq(-10, 10, 2)) +
facet_grid(gene_set ~ ., scales = "free_y", space = "free_y",
switch = "y") +
theme_minimal(base_size = 20) +
theme(panel.spacing.y = unit(1, "lines"),
strip.placement = "outside",
strip.background = element_blank(),
strip.text.y.left = element_text(size = 20, face = "bold",
angle = 90),
axis.text.y = element_text(lineheight = 0.8, vjust = 0.5),
plot.margin = margin(t = 30, r = 5, b = 5, l = 5)) +
labs(y = NULL,
x = expression(Delta * "kWithin" ~ "(" * italic("Inv4m") *
" - CTRL, Z-score)"))
print(p_kw)
ggsave(file.path(OUTPUT_DIR, "zscore_density_kWithin.pdf"),
p_kw, width = 10, height = 7)
ggsave(file.path(OUTPUT_DIR, "zscore_density_kWithin.png"),
p_kw, width = 10, height = 7, dpi = 300)
saveRDS(p_kw, file.path(OUTPUT_DIR, "zscore_density_kWithin.rds"))Z-score Ridge Plots (kOut)
Same approach for inter-modular connectivity (\(\Delta k_{Out}\)).
# Z-score kOut
gr_null_mean_out <- colMeans(null_gr_out)
gr_null_sd_out <- apply(null_gr_out, 2, sd)
gr_null_z_out <- null_gr_out
for (m in colnames(gr_null_z_out)) {
gr_null_z_out[, m] <- (gr_null_z_out[, m] - gr_null_mean_out[m]) / gr_null_sd_out[m]
}
gr_obs_z_out <- (observed_gr_out - gr_null_mean_out) / gr_null_sd_out
lg_null_mean_out <- colMeans(null_lg_out)
lg_null_sd_out <- apply(null_lg_out, 2, sd)
lg_null_z_out <- null_lg_out
for (m in colnames(lg_null_z_out)) {
lg_null_z_out[, m] <- (lg_null_z_out[, m] - lg_null_mean_out[m]) / lg_null_sd_out[m]
}
lg_obs_z_out <- (observed_lg_out - lg_null_mean_out) / lg_null_sd_out
# Build combined density data
gr_density_out <- map_dfr(colnames(gr_null_z_out), function(m) {
tibble(module = sprintf("%s (n=%d)", m, gr_mod_sizes[m]),
fill_hex = m, z_value = gr_null_z_out[, m],
gene_set = "Genotype Response DEGs")
})
lg_density_out <- map_dfr(colnames(lg_null_z_out), function(m) {
tibble(module = sprintf("%s (n=%d)", m, lg_mod_sizes[m]),
fill_hex = lg_color_hex[m], z_value = lg_null_z_out[, m],
gene_set = "Leaf Gradient Hubs")
})
density_kout <- bind_rows(gr_density_out, lg_density_out)
# Build combined points data
gr_points_out <- tibble(
module = sprintf("%s (n=%d)", names(gr_obs_z_out), gr_mod_sizes[names(gr_obs_z_out)]),
fill_hex = names(gr_obs_z_out),
z_observed = gr_obs_z_out, fdr = kout_fdr_gr[names(gr_obs_z_out)],
gene_set = "Genotype Response DEGs"
)
lg_points_out <- tibble(
module = sprintf("%s (n=%d)", names(lg_obs_z_out), lg_mod_sizes[names(lg_obs_z_out)]),
fill_hex = lg_color_hex[names(lg_obs_z_out)],
z_observed = lg_obs_z_out, fdr = kout_fdr_lg[names(lg_obs_z_out)],
gene_set = "Leaf Gradient Hubs"
)
points_kout <- bind_rows(gr_points_out, lg_points_out)
# Order modules by z_observed WITHIN each gene_set, compute per-facet rank
points_kout <- points_kout |>
group_by(gene_set) |>
arrange(z_observed, .by_group = TRUE) |>
mutate(facet_rank = row_number()) |>
ungroup()
mod_order_out <- points_kout |> arrange(gene_set, facet_rank) |> pull(module)
density_kout$module <- factor(density_kout$module, levels = mod_order_out)
points_kout$module <- factor(points_kout$module, levels = mod_order_out)
density_kout <- density_kout |>
left_join(points_kout |> dplyr::select(module, gene_set, facet_rank) |> distinct(),
by = c("module", "gene_set"))
fill_map_out <- setNames(points_kout$fill_hex, as.character(points_kout$module))
# GO annotations with per-facet rank
go_df_out <- bind_rows(gr_go_df, lg_go_df) |>
mutate(module = factor(module, levels = levels(density_kout$module))) |>
left_join(points_kout |> dplyr::select(module, gene_set, facet_rank) |> distinct(),
by = c("module", "gene_set")) |>
mutate(y_pos = facet_rank + 0.35)
p_kout <- ggplot(density_kout, aes(x = z_value, y = module, fill = module)) +
geom_density_ridges(color = "grey40", alpha = 0.7,
scale = 0.9, rel_min_height = 0.01) +
scale_fill_manual(values = fill_map_out) +
guides(fill = "none") +
geom_vline(xintercept = 0, linetype = "dashed", color = "grey50") +
geom_segment(
data = points_kout,
aes(x = z_observed, xend = z_observed,
y = facet_rank, yend = facet_rank + 0.8),
color = "red", linewidth = 0.8, inherit.aes = FALSE
) +
geom_text(
data = points_kout |> filter(fdr < 0.05),
aes(x = z_observed - 0.6, y = facet_rank + 0.35,
label = sprintf("FDR=%.3f", fdr)),
color = "red", size = 3, hjust = 0.5, inherit.aes = FALSE
) +
geom_text(
data = go_df_out,
aes(x = 2.5, y = y_pos, label = go_desc),
hjust = 0, size = 14 * 0.8 / .pt, fontface = "italic",
inherit.aes = FALSE
) +
scale_x_continuous(limits = c(NA, 5), breaks = seq(-10, 10, 2)) +
facet_grid(gene_set ~ ., scales = "free_y", space = "free_y") +
theme_minimal(base_size = 14) +
theme(panel.spacing.y = unit(1, "lines"),
strip.text.y = element_text(size = 14)) +
labs(y = NULL, x = expression(Delta * k[Out] ~ "(Inv4m - CTRL, Z-score)"))
print(p_kout)
ggsave(file.path(OUTPUT_DIR, "zscore_density_kOut.pdf"),
p_kout, width = 10, height = 7)
ggsave(file.path(OUTPUT_DIR, "zscore_density_kOut.png"),
p_kout, width = 10, height = 7, dpi = 300)kOut: CTRL vs Inv4m Gene-level Scatter
Each point is a gene; axes are kOut in CTRL vs kOut in Inv4m. Points above the diagonal gained inter-modular connectivity. Colored by module.
# Gene-level kOut from genotype-response adjacencies
adj_ctrl_gr <- adjacency(datExpr_gr[ctrl_idx, ], power = SOFT_POWER, type = "signed")
adj_inv4m_gr <- adjacency(datExpr_gr[inv4m_idx, ], power = SOFT_POWER, type = "signed")
k_ctrl_gr <- intramodularConnectivity(adj_ctrl_gr, gr_colors)
k_inv4m_gr <- intramodularConnectivity(adj_inv4m_gr, gr_colors)
kout_scatter <- tibble(
gene = rownames(k_ctrl_gr),
module = gr_colors[rownames(k_ctrl_gr)],
kOut_CTRL = k_ctrl_gr$kOut,
kOut_Inv4m = k_inv4m_gr$kOut
) |>
filter(module != "grey", !is.na(module))
p_kout_scatter <- ggplot(
kout_scatter,
aes(x = kOut_CTRL, y = kOut_Inv4m, color = module)
) +
geom_point(alpha = 0.6, size = 2) +
geom_abline(slope = 1, intercept = 0,
linetype = "dashed", color = "grey50") +
scale_color_identity() +
theme_minimal(base_size = 14) +
labs(
x = expression(k[Out] ~ "(CTRL)"),
y = expression(k[Out] ~ "(Inv4m)"),
title = "Inter-modular connectivity: CTRL vs Inv4m"
) +
coord_equal()
print(p_kout_scatter)
cat("Genes above diagonal (kOut gain):",
sum(kout_scatter$kOut_Inv4m > kout_scatter$kOut_CTRL),
"of", nrow(kout_scatter), "\n")## Genes above diagonal (kOut gain): 218 of 405wt_kout <- wilcox.test(kout_scatter$kOut_Inv4m,
kout_scatter$kOut_CTRL,
paired = TRUE)
cat("Paired Wilcoxon p:", format.pval(wt_kout$p.value), "\n")## Paired Wilcoxon p: 0.012449ggsave(file.path(OUTPUT_DIR, "kOut_CTRL_vs_Inv4m_scatter.pdf"),
p_kout_scatter, width = 8, height = 8)
ggsave(file.path(OUTPUT_DIR, "kOut_CTRL_vs_Inv4m_scatter.png"),
p_kout_scatter, width = 8, height = 8, dpi = 300)|logFC| vs Delta kWithin
Scatterplot of Inv4m expression effect size against connectivity change, for genotype-response and leaf-gradient genes separately.
# Load DEG effects
effects <- read_csv(
file.path(paths$intermediate, "predictor_effects_leaf_interaction_model.csv"),
show_col_types = FALSE
)
inv4m_effects <- effects |>
filter(grepl("Inv4m", predictor), !grepl(":", predictor)) |>
dplyr::select(gene, logFC, adj.P.Val)
# Gene-level delta kWithin for genotype-response genes
adj_ctrl_gr <- adjacency(datExpr_gr[ctrl_idx, ], power = SOFT_POWER, type = "signed")
adj_inv4m_gr <- adjacency(datExpr_gr[inv4m_idx, ], power = SOFT_POWER, type = "signed")
k_ctrl_gr <- intramodularConnectivity(adj_ctrl_gr, gr_colors)
k_inv4m_gr <- intramodularConnectivity(adj_inv4m_gr, gr_colors)
gr_dk <- tibble(
gene = rownames(k_ctrl_gr),
delta_kWithin = k_inv4m_gr$kWithin - k_ctrl_gr$kWithin,
group = "Genotype Response DEGs"
)
# Gene-level delta kWithin for leaf-gradient genes
adj_ctrl_lg <- adjacency(datExpr_lg[ctrl_idx, ], power = SOFT_POWER, type = "signed")
adj_inv4m_lg <- adjacency(datExpr_lg[inv4m_idx, ], power = SOFT_POWER, type = "signed")
k_ctrl_lg <- intramodularConnectivity(adj_ctrl_lg, lg_colors)
k_inv4m_lg <- intramodularConnectivity(adj_inv4m_lg, lg_colors)
lg_dk <- tibble(
gene = names(lg_colors),
delta_kWithin = k_inv4m_lg$kWithin - k_ctrl_lg$kWithin,
group = "Leaf Gradient Hubs"
)
# Merge with logFC
scatter_df <- bind_rows(gr_dk, lg_dk) |>
left_join(inv4m_effects |> dplyr::select(gene, logFC, adj.P.Val), by = "gene") |>
mutate(
abs_logFC = abs(logFC),
neg_log10p = -log10(adj.P.Val)
) |>
filter(!is.na(logFC))
# Correlations per group
cor_logfc <- scatter_df |>
group_by(group) |>
summarise(rho = cor(abs_logFC, delta_kWithin, method = "spearman"),
.groups = "drop")
cat("Spearman correlations (|logFC| vs delta_k):\n")## Spearman correlations (|logFC| vs delta_k):print(cor_logfc)## # A tibble: 2 × 2
## group rho
## <chr> <dbl>
## 1 Genotype Response DEGs 0.162
## 2 Leaf Gradient Hubs -0.0241# Panel A: |logFC| vs delta kWithin
p_logfc <- ggplot(scatter_df, aes(x = abs_logFC, y = delta_kWithin, color = group)) +
geom_point(alpha = 0.4, size = 1.5) +
geom_smooth(method = "lm", se = TRUE, linewidth = 1) +
geom_hline(yintercept = 0, linetype = "dashed", color = "grey50") +
scale_color_manual(values = c("Genotype Response DEGs" = "grey40",
"Leaf Gradient Hubs" = "steelblue")) +
theme_minimal(base_size = 14) +
labs(
x = expression("|" * log[2] * "FC" * "| (Inv4m vs CTRL)"),
y = expression(Delta ~ k[Within] ~ "(Inv4m - CTRL)"),
color = NULL
) +
theme(legend.position = c(0.2, 0.15))
# Panel B: -log10P vs delta kWithin
p_pval <- ggplot(scatter_df, aes(x = neg_log10p, y = delta_kWithin, color = group)) +
geom_point(alpha = 0.4, size = 1.5) +
geom_smooth(method = "lm", se = TRUE, linewidth = 1) +
geom_hline(yintercept = 0, linetype = "dashed", color = "grey50") +
scale_color_manual(values = c("Genotype Response DEGs" = "grey40",
"Leaf Gradient Hubs" = "steelblue")) +
theme_minimal(base_size = 14) +
labs(
x = expression("-" * log[10] * "(FDR)"),
y = expression(Delta ~ k[Within] ~ "(Inv4m - CTRL)"),
color = NULL
) +
theme(legend.position = "none")
p_combined <- p_logfc + p_pval + plot_annotation(tag_levels = "A")
print(p_combined)
ggsave(file.path(OUTPUT_DIR, "logFC_pval_vs_delta_kWithin.pdf"),
p_combined, width = 14, height = 4)
ggsave(file.path(OUTPUT_DIR, "logFC_pval_vs_delta_kWithin.png"),
p_combined, width = 14, height = 6, dpi = 300)Linear model: does expression effect size predict connectivity loss?
# Genotype-response only
gr_only <- scatter_df |> filter(group == "Genotype Response DEGs")
fit_gr <- lm(delta_kWithin ~ abs_logFC, data = gr_only)
cat("=== Genotype-response genes: lm(delta_kWithin ~ |logFC|) ===\n")## === Genotype-response genes: lm(delta_kWithin ~ |logFC|) ===summary(fit_gr)##
## Call:
## lm(formula = delta_kWithin ~ abs_logFC, data = gr_only)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.1833 -0.4561 0.4889 0.8940 2.6496
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.95453 0.08608 -11.088 <2e-16 ***
## abs_logFC 0.01514 0.03322 0.456 0.649
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.389 on 463 degrees of freedom
## Multiple R-squared: 0.0004487, Adjusted R-squared: -0.00171
## F-statistic: 0.2078 on 1 and 463 DF, p-value: 0.6487cat(sprintf("\nR2 = %.4f (|logFC| explains %.1f%% of variance in delta_k)\n",
summary(fit_gr)$r.squared,
summary(fit_gr)$r.squared * 100))##
## R2 = 0.0004 (|logFC| explains 0.0% of variance in delta_k)# Leaf-gradient only
lg_only <- scatter_df |> filter(group == "Leaf Gradient Hubs")
fit_lg <- lm(delta_kWithin ~ abs_logFC, data = lg_only)
cat("\n=== Leaf-gradient genes: lm(delta_kWithin ~ |logFC|) ===\n")##
## === Leaf-gradient genes: lm(delta_kWithin ~ |logFC|) ===summary(fit_lg)##
## Call:
## lm(formula = delta_kWithin ~ abs_logFC, data = lg_only)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.558 -4.842 2.183 5.271 14.811
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.7408 0.7182 -9.385 <2e-16 ***
## abs_logFC 1.3936 3.4412 0.405 0.686
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.2 on 216 degrees of freedom
## Multiple R-squared: 0.0007587, Adjusted R-squared: -0.003867
## F-statistic: 0.164 on 1 and 216 DF, p-value: 0.6859cat(sprintf("\nR2 = %.4f (|logFC| explains %.1f%% of variance in delta_k)\n",
summary(fit_lg)$r.squared,
summary(fit_lg)$r.squared * 100))##
## R2 = 0.0008 (|logFC| explains 0.1% of variance in delta_k)# Combined model with interaction
fit_interaction <- lm(delta_kWithin ~ abs_logFC * group, data = scatter_df)
cat("\n=== Combined: lm(delta_kWithin ~ |logFC| * group) ===\n")##
## === Combined: lm(delta_kWithin ~ |logFC| * group) ===summary(fit_interaction)##
## Call:
## lm(formula = delta_kWithin ~ abs_logFC * group, data = scatter_df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -17.5580 -0.8541 0.5735 1.1300 14.8108
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.95453 0.26148 -3.650 0.000282 ***
## abs_logFC 0.01514 0.10090 0.150 0.880744
## groupLeaf Gradient Hubs -5.78622 0.49555 -11.676 < 2e-16 ***
## abs_logFC:groupLeaf Gradient Hubs 1.37841 2.01935 0.683 0.495093
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.22 on 679 degrees of freedom
## Multiple R-squared: 0.2783, Adjusted R-squared: 0.2751
## F-statistic: 87.26 on 3 and 679 DF, p-value: < 2.2e-16cat(sprintf("\nR2 = %.4f\n", summary(fit_interaction)$r.squared))##
## R2 = 0.2783Combined Statistics Table
# Genotype-response table
gr_pres <- gr_preservation |>
filter(module_label != "gold", module_label != "grey") |>
dplyr::select(module = module_label, Zsummary, Zdensity, Zconnectivity)
gr_table <- gr_pres |>
left_join(
tibble(
module = names(observed_gr_kw),
delta_kWithin = round(observed_gr_kw, 3),
kWithin_FDR = round(kw_fdr_gr[names(observed_gr_kw)], 4),
delta_kOut = round(observed_gr_out, 3),
kOut_FDR = round(kout_fdr_gr[names(observed_gr_out)], 4)
),
by = "module"
) |>
mutate(
gene_set = "Genotype Response DEGs",
n_genes = sapply(module, function(m)
sum(gr_colors == m, na.rm = TRUE)),
LFP = round(gr_lfp$mean_lfp[match(module, gr_lfp$color)], 3),
Zsummary = round(Zsummary, 2),
Zdensity = round(Zdensity, 2),
Zconnectivity = round(Zconnectivity, 2)
)
# Leaf-gradient table
lg_pres <- lg_preservation |>
filter(module_label != "gold", module_label != "grey") |>
dplyr::select(module = module_label, Zsummary, Zdensity, Zconnectivity)
# Map leaf-gradient module colors through the remap
lg_pres$module <- lg_color_remap[lg_pres$module]
lg_table <- lg_pres |>
left_join(
tibble(
module = names(observed_lg_kw),
delta_kWithin = round(observed_lg_kw, 3),
kWithin_FDR = round(kw_fdr_lg[names(observed_lg_kw)], 4),
delta_kOut = round(observed_lg_out, 3),
kOut_FDR = round(kout_fdr_lg[names(observed_lg_out)], 4)
),
by = "module"
) |>
mutate(
gene_set = "Leaf Gradient Hubs",
n_genes = sapply(module, function(m)
sum(lg_colors == m)),
LFP = round(lg_lfp$mean_lfp[match(
names(lg_color_remap)[match(module, lg_color_remap)],
lg_lfp$color
)], 3),
Zsummary = round(Zsummary, 2),
Zdensity = round(Zdensity, 2),
Zconnectivity = round(Zconnectivity, 2)
)
# Combine
combined_table <- bind_rows(gr_table, lg_table) |>
dplyr::select(
gene_set, module, n_genes, LFP,
Zsummary, Zdensity, Zconnectivity,
delta_kWithin, kWithin_FDR,
delta_kOut, kOut_FDR
) |>
arrange(gene_set, kWithin_FDR)
cat("=== COMBINED MODULE STATISTICS ===\n\n")## === COMBINED MODULE STATISTICS ===print(combined_table, n = 30, width = Inf)## # A tibble: 18 × 11
## gene_set module n_genes LFP Zsummary Zdensity
## <chr> <chr> <int> <dbl> <dbl> <dbl>
## 1 Genotype Response DEGs brown 46 0.683 3.85 3.93
## 2 Genotype Response DEGs tan 12 0.684 2.24 2.47
## 3 Genotype Response DEGs purple 22 0.74 3.84 4.42
## 4 Genotype Response DEGs blue 51 0.795 3.11 1.86
## 5 Genotype Response DEGs pink 30 0.715 5.77 5.12
## 6 Genotype Response DEGs yellow 36 0.682 3.49 2.58
## 7 Genotype Response DEGs salmon 11 0.52 1.2 0.66
## 8 Genotype Response DEGs greenyellow 16 0.798 2.96 3.51
## 9 Genotype Response DEGs red 34 0.63 3.31 4.41
## 10 Genotype Response DEGs black 32 0.694 4.51 5.07
## 11 Genotype Response DEGs green 36 0.688 5.01 7.82
## 12 Genotype Response DEGs magenta 24 0.713 5.33 7.74
## 13 Genotype Response DEGs turquoise 55 0.537 6.05 6.52
## 14 Leaf Gradient Hubs coral 26 0.951 4.28 3.76
## 15 Leaf Gradient Hubs jade 30 0.948 4.37 2.49
## 16 Leaf Gradient Hubs sky 78 0.83 7.91 5.11
## 17 Leaf Gradient Hubs lavender 13 0.805 4.09 3.64
## 18 Leaf Gradient Hubs peach 64 0.864 8.33 8.36
## Zconnectivity delta_kWithin kWithin_FDR delta_kOut kOut_FDR
## <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 3.77 -2.36 0 -0.349 0.721
## 2 2.01 -3.01 0.0099 1.09 0.721
## 3 3.26 -1.11 0.0336 -1.34 0.721
## 4 4.37 -1.64 0.0338 0.278 0.721
## 5 6.41 -1.36 0.0338 -1.51 0.721
## 6 4.4 -1.47 0.0846 -0.725 0.721
## 7 1.74 -0.389 0.220 1.62 0.721
## 8 2.41 -0.631 0.220 -0.436 0.721
## 9 2.2 -0.508 0.220 0.885 0.721
## 10 3.95 -0.801 0.220 0.222 0.721
## 11 2.21 -0.612 0.264 1.58 0.721
## 12 2.92 -0.143 0.383 0.721 0.721
## 13 5.59 -0.277 0.441 1.11 0.721
## 14 4.79 -2.35 0.220 -1.63 0.721
## 15 6.25 -4.39 0.246 -1.44 0.721
## 16 10.7 -13.8 0.262 -4.09 0.721
## 17 4.55 -1.36 0.340 -8.98 0.721
## 18 8.3 -2.09 0.383 -3.16 0.721write_csv(combined_table,
file.path(OUTPUT_DIR, "combined_module_stats_table.csv"))
cat("\nSaved: combined_module_stats_table.csv\n")##
## Saved: combined_module_stats_table.csvSummary
n_gr_sig <- sum(all_results$gene_set == "Genotype Response DEGs" & all_results$fdr < 0.05)
n_lg_sig <- sum(all_results$gene_set == "Leaf Gradient Hubs" & all_results$fdr < 0.05)
cat("=== PERMUTATION TEST SUMMARY ===
Genotype-response modules significant (FDR < 0.05):", n_gr_sig, "of",
sum(all_results$gene_set == "Genotype Response DEGs"), "
Leaf-gradient modules significant (FDR < 0.05):", n_lg_sig, "of",
sum(all_results$gene_set == "Leaf Gradient Hubs"), "
Significant genotype-response modules:\n")## === PERMUTATION TEST SUMMARY ===
##
## Genotype-response modules significant (FDR < 0.05): 5 of 13
## Leaf-gradient modules significant (FDR < 0.05): 0 of 5
##
## Significant genotype-response modules:all_results |>
filter(gene_set == "Genotype Response DEGs", fdr < 0.05) |>
dplyr::select(module, n_genes, observed_mean_delta_k, ecdf_pvalue, fdr) |>
print()## # A tibble: 5 × 5
## module n_genes observed_mean_delta_k ecdf_pvalue fdr
## <chr> <int> <dbl> <dbl> <dbl>
## 1 brown 46 -2.36 0 0
## 2 tan 12 -3.01 0.0011 0.0099
## 3 purple 22 -1.11 0.0056 0.0336
## 4 blue 51 -1.64 0.0078 0.0338
## 5 pink 30 -1.36 0.0094 0.0338n_gr_sig_out <- sum(kout_fdr_gr < 0.05)
n_lg_sig_out <- sum(kout_fdr_lg < 0.05)
cat("
=== kOut (inter-modular) ===
Genotype-response modules with significant kOut gain (FDR < 0.05):",
n_gr_sig_out, "of", length(kout_fdr_gr), "
Leaf-gradient modules with significant kOut gain (FDR < 0.05):",
n_lg_sig_out, "of", length(kout_fdr_lg), "
=== CONCLUSIONS ===
1. kWithin: Genotype-response modules show significant intramodular
connectivity loss beyond the permutation null.
Leaf-gradient modules do not.
2. kOut: No significant gain in inter-modular connectivity
for genotype-response or leaf-gradient modules. The perturbation
degrades within-module coexpression structure without causing
genes to gain spurious cross-module connections.
3. This asymmetry (kWithin loss without kOut gain) argues
against a global artifact of mean expression shifts:
if the connectivity change were driven by correlated
mean shifts inflating all pairwise differences, both
kWithin and kOut would be affected. Instead, only
within-module relationships are disrupted, consistent
with degradation of specific coexpression programs.
4. The linear model confirms this: |logFC| does not
predict delta kWithin, ruling out differential expression
magnitude as a driver of connectivity loss.
")##
## === kOut (inter-modular) ===
## Genotype-response modules with significant kOut gain (FDR < 0.05): 0 of 13
## Leaf-gradient modules with significant kOut gain (FDR < 0.05): 0 of 5
##
## === CONCLUSIONS ===
##
## 1. kWithin: Genotype-response modules show significant intramodular
## connectivity loss beyond the permutation null.
## Leaf-gradient modules do not.
##
## 2. kOut: No significant gain in inter-modular connectivity
## for genotype-response or leaf-gradient modules. The perturbation
## degrades within-module coexpression structure without causing
## genes to gain spurious cross-module connections.
##
## 3. This asymmetry (kWithin loss without kOut gain) argues
## against a global artifact of mean expression shifts:
## if the connectivity change were driven by correlated
## mean shifts inflating all pairwise differences, both
## kWithin and kOut would be affected. Instead, only
## within-module relationships are disrupted, consistent
## with degradation of specific coexpression programs.
##
## 4. The linear model confirms this: |logFC| does not
## predict delta kWithin, ruling out differential expression
## magnitude as a driver of connectivity loss.Session Info
sessionInfo()## R version 4.5.2 (2025-10-31)
## Platform: aarch64-apple-darwin20
## Running under: macOS Tahoe 26.3.1
##
## Matrix products: default
## BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
##
## locale:
## [1] C.UTF-8/C.UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8
##
## time zone: America/New_York
## tzcode source: internal
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] AnnotationDbi_1.72.0 IRanges_2.44.0 S4Vectors_0.48.0
## [4] Biobase_2.70.0 AnnotationHub_4.0.0 BiocFileCache_3.0.0
## [7] dbplyr_2.5.1 BiocGenerics_0.56.0 generics_0.1.4
## [10] rrvgo_1.22.0 clusterProfiler_4.18.4 patchwork_1.3.2
## [13] ggridges_0.5.7 WGCNA_1.73 fastcluster_1.3.0
## [16] dynamicTreeCut_1.63-1 lubridate_1.9.4 forcats_1.0.1
## [19] stringr_1.6.0 dplyr_1.1.4 purrr_1.2.0
## [22] readr_2.1.6 tidyr_1.3.2 tibble_3.3.0
## [25] ggplot2_4.0.1 tidyverse_2.0.0 here_1.0.2
##
## loaded via a namespace (and not attached):
## [1] splines_4.5.2 later_1.4.4 filelock_1.0.3
## [4] ggplotify_0.1.3 R.oo_1.27.1 polyclip_1.10-7
## [7] preprocessCore_1.72.0 rpart_4.1.24 httr2_1.2.2
## [10] lifecycle_1.0.4 doParallel_1.0.17 rprojroot_2.1.1
## [13] vroom_1.6.7 NLP_0.3-2 lattice_0.22-7
## [16] MASS_7.3-65 backports_1.5.0 magrittr_2.0.4
## [19] Hmisc_5.2-4 sass_0.4.10 rmarkdown_2.30
## [22] jquerylib_0.1.4 yaml_2.3.12 httpuv_1.6.16
## [25] otel_0.2.0 ggtangle_0.0.9 askpass_1.2.1
## [28] reticulate_1.44.1 cowplot_1.2.0 DBI_1.2.3
## [31] RColorBrewer_1.1-3 R.utils_2.13.0 yulab.utils_0.2.3
## [34] nnet_7.3-20 tweenr_2.0.3 rappdirs_0.3.3
## [37] gdtools_0.4.4 enrichplot_1.30.4 tm_0.7-16
## [40] ggrepel_0.9.6 tidytree_0.4.6 pheatmap_1.0.13
## [43] umap_0.2.10.0 RSpectra_0.16-2 codetools_0.2-20
## [46] DOSE_4.4.0 xml2_1.5.1 ggforce_0.5.0
## [49] tidyselect_1.2.1 aplot_0.2.9 farver_2.1.2
## [52] matrixStats_1.5.0 base64enc_0.1-3 Seqinfo_1.0.0
## [55] jsonlite_2.0.0 Formula_1.2-5 survival_3.8-3
## [58] iterators_1.0.14 systemfonts_1.3.1 foreach_1.5.2
## [61] tools_4.5.2 ggnewscale_0.5.2 ragg_1.5.0
## [64] treeio_1.34.0 Rcpp_1.1.0 glue_1.8.0
## [67] gridExtra_2.3 mgcv_1.9-4 xfun_0.55
## [70] qvalue_2.42.0 withr_3.0.2 BiocManager_1.30.27
## [73] fastmap_1.2.0 openssl_2.3.4 digest_0.6.39
## [76] timechange_0.3.0 R6_2.6.1 mime_0.13
## [79] gridGraphics_0.5-1 textshaping_1.0.4 colorspace_2.1-2
## [82] GO.db_3.22.0 dichromat_2.0-0.1 RSQLite_2.4.5
## [85] R.methodsS3_1.8.2 utf8_1.2.6 fontLiberation_0.1.0
## [88] data.table_1.18.0 httr_1.4.7 htmlwidgets_1.6.4
## [91] scatterpie_0.2.6 pkgconfig_2.0.3 gtable_0.3.6
## [94] blob_1.2.4 S7_0.2.1 impute_1.84.0
## [97] XVector_0.50.0 htmltools_0.5.9 fontBitstreamVera_0.1.1
## [100] fgsea_1.36.0 scales_1.4.0 png_0.1-8
## [103] wordcloud_2.6 ggfun_0.2.0 knitr_1.51
## [106] rstudioapi_0.17.1 tzdb_0.5.0 reshape2_1.4.5
## [109] curl_7.0.0 checkmate_2.3.3 nlme_3.1-168
## [112] cachem_1.1.0 BiocVersion_3.22.0 parallel_4.5.2
## [115] foreign_0.8-90 treemap_2.4-4 pillar_1.11.1
## [118] grid_4.5.2 vctrs_0.6.5 slam_0.1-55
## [121] promises_1.5.0 tidydr_0.0.6 xtable_1.8-4
## [124] cluster_2.1.8.1 htmlTable_2.4.3 evaluate_1.0.5
## [127] cli_3.6.5 compiler_4.5.2 rlang_1.1.6
## [130] crayon_1.5.3 labeling_0.4.3 plyr_1.8.9
## [133] fs_1.6.6 ggiraph_0.9.2 stringi_1.8.7
## [136] gridBase_0.4-7 BiocParallel_1.44.0 Biostrings_2.78.0
## [139] lazyeval_0.2.2 GOSemSim_2.36.0 fontquiver_0.2.1
## [142] Matrix_1.7-4 hms_1.1.4 bit64_4.6.0-1
## [145] KEGGREST_1.50.0 shiny_1.12.1 igraph_2.2.1
## [148] memoise_2.0.1 bslib_0.9.0 ggtree_4.0.3
## [151] fastmatch_1.1-6 bit_4.6.0 ape_5.8-1
## [154] gson_0.1.0