Co-expression Network

Build a co-expression network; find modules and hub genes.

Overview

Problem. Move from single genes to co-varying modules.

Learning goals

• From gene lists to network thinking
• High connectivity ≠ causation

Figures

Tutorial

Overview

Build weighted gene co-expression networks to identify modules of coordinately expressed genes and discover hub genes that may be key regulators. This workflow uses WGCNA (Weighted Gene Co-expression Network Analysis) to group genes into modules based on their expression patterns across samples, then correlates these modules with experimental conditions or traits.

Key Concept: Unlike single-gene analysis, WGCNA identifies groups of genes that behave similarly across samples, revealing biological pathways and potential regulatory relationships.

Default Prompt: "Build a co-expression network to identify gene modules and hub genes from my RNA-seq data"

When to Use This Skill

Use WGCNA when you want to:

• Identify gene modules associated with experimental conditions or phenotypes
• Discover hub genes that are highly connected within modules and may be key regulators
• Find co-expressed genes with similar expression patterns to known genes of interest
• Reduce dimensionality of large gene expression datasets for downstream analysis
• Generate hypotheses about gene function based on co-expression patterns

Installation

Core WGCNA packages:

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("WGCNA")

Visualization packages:

install.packages(c("ggplot2", "ggprism"))
BiocManager::install("ComplexHeatmap")

Enrichment analysis (optional):

BiocManager::install(c("clusterProfiler", "org.Hs.eg.db"))  # Human
# BiocManager::install("org.Mm.eg.db")  # Mouse
# BiocManager::install("org.Rn.eg.db")  # Rat

Inputs

2. Sample metadata (CSV/TSV): - Sample IDs matching expression matrix columns - Traits/conditions for module-trait correlation

Outputs

Clarification Questions

1. Input Files (ASK THIS FIRST):

   • Do you have specific normalized expression data and sample metadata files to analyze?
   • If uploaded: Are these the expression matrix and metadata you'd like to analyze?
   • Expected formats: CSV or TSV with genes as rows, samples as columns
   • Or use example data? Female mouse liver dataset (135 samples, liver tissue, multiple traits)

2. What is your normalized expression data format? - VST (variance stabilizing transformation) from DESeq2

   • rlog (regularized log) from DESeq2
   • TPM (transcripts per million)
   • FPKM/RPKM
   • If unsure or raw counts: normalize first using DESeq2

3. How many samples do you have? - 15-30 samples (minimum for WGCNA, results may be less robust)

   • 30-50 samples (good power for network detection)
   • 50+ samples (excellent power, most reliable results)

4. What traits/conditions do you want to correlate with modules? - Treatment vs control (binary)

   • Disease status or phenotype
   • Continuous variables (age, dose, time, weight)
   • Multiple traits (all will be tested)

5. Gene filtering strategy? - Top 5,000 most variable genes (default, recommended)

   • Top 10,000-15,000 genes (for larger datasets)
   • All genes passing expression threshold
   • Pre-filtered gene list (e.g., from DE analysis)

Standard Workflow

🚨 MANDATORY: USE SCRIPTS EXACTLY AS SHOWN - DO NOT WRITE INLINE CODE 🚨

Step 1 - Load data:

library(WGCNA)
allowWGCNAThreads()

source("scripts/load_example_data.R")
wgcna_data <- load_example_wgcna_data()
datExpr <- wgcna_data$datExpr
meta <- wgcna_data$meta

# For your own data:
# source("scripts/prepare_wgcna_data.R")
# data <- prepare_wgcna_data("expression.csv", "metadata.csv", top_n_genes = 5000)
# datExpr <- data$datExpr
# meta <- data$meta

Step 2 - Run WGCNA analysis:

source("scripts/wgcna_workflow.R")
results <- run_wgcna_analysis(
  datExpr,
  meta,
  traits = c("weight_g", "Glucose_mg_dl"),  # Adjust to your traits
  organism = "mouse"  # or "human", "rat", or NULL to skip enrichment
)

Step 3 - Generate visualizations:

source("scripts/plot_all_wgcna.R")
plot_all_wgcna(results, output_dir = "wgcna_results")

The script handles PNG + SVG export with graceful fallback for SVG dependencies.

Step 4 - Export results:

source("scripts/export_wgcna_results.R")
export_all(results, output_dir = "wgcna_results")

NEVER skip directly to writing inline code without trying the script first.

Parameter Customization

When customization is needed:

• Soft power selection: Read references/parameter-tuning-guide.md to understand how to choose appropriate power values for your data
• Module detection parameters: See references/parameter-tuning-guide.md#module-detection for guidance on min_module_size and merge_cut_height
• Complete custom workflow: Read references/wgcna-reference.md for detailed code examples with explanations (only if you need full control)

Common Issues

Interpretation Guidelines

Suggested Next Steps

After identifying modules and hub genes:

1. Functional validation - Validate hub genes experimentally (qPCR, knockdown, overexpression)
2. Enrichment analysis - Test modules for GO/KEGG enrichment to understand biological processes
3. Compare with DE results - Overlay DE genes on network to see which modules are enriched
4. Network visualization - Export to Cytoscape for detailed network visualization
5. Cross-dataset validation - Test module preservation in independent datasets

Related Skills

• bulk-rnaseq-counts-to-de-deseq2 - Normalize counts and perform differential expression analysis (run before WGCNA)
• de-results-to-gene-lists - Extract gene lists from DE results to overlay on network
• functional-enrichment-from-degs - Perform GO/KEGG enrichment on modules

References

Code preview

scripts/build_network.R

# Build co-expression network and detect modules

library(WGCNA)

#' Build co-expression network and detect modules
#'
#' @param datExpr Expression matrix (samples x genes)
#' @param power Soft-thresholding power
#' @param min_module_size Minimum genes per module
#' @param merge_cut_height Height for merging similar modules
#' @return List containing network object, module colors, and module labels
build_network <- function(datExpr, power, min_module_size = 30, merge_cut_height = 0.25) {

  cat("Building network with power =", power, "\n")

  # One-step network construction and module detection
  net <- blockwiseModules(
    datExpr,
    power = power,
    TOMType = "signed",
    networkType = "signed",
    minModuleSize = min_module_size,
    reassignThreshold = 0,
    mergeCutHeight = merge_cut_height,
    numericLabels = TRUE,
    pamRespectsDendro = FALSE,
    saveTOMs = FALSE,
    verbose = 3
  )

  # Convert numeric labels to colors
  module_colors <- labels2colors(net$colors)

  cat("\nModule detection complete:\n")
  cat("Number of modules:", length(unique(module_colors)) - 1, "(excluding grey/unassigned)\n")
  print(table(module_colors))

  return(list(
    net = net,
    module_colors = module_colors,
    module_labels = net$colors
  ))
}

scripts/correlate_modules_traits.R

# Correlate module eigengenes with sample traits

library(WGCNA)
library(ComplexHeatmap)
library(circlize)

#' Correlate module eigengenes with sample traits
#'
#' @param datExpr Expression matrix
#' @param module_colors Module color assignments
#' @param traits Data frame of sample traits (numeric)
#' @param output_file Path to save heatmap
#' @return List with module eigengenes, correlations, and p-values
correlate_modules_traits <- function(datExpr, module_colors, traits,
                                      output_file = "module_trait_correlation.svg") {

  # Calculate module eigengenes
  MEs <- moduleEigengenes(datExpr, colors = module_colors)$eigengenes
  MEs <- orderMEs(MEs)

  # Convert traits to numeric if needed
  traits_numeric <- data.frame(lapply(traits, function(x) {
    if (is.factor(x) || is.character(x)) {
      as.numeric(as.factor(x))
    } else {
      as.numeric(x)
    }
  }))
  rownames(traits_numeric) <- rownames(traits)

  # Calculate correlations
  module_trait_cor <- cor(MEs, traits_numeric, use = "pairwise.complete.obs")
  module_trait_pval <- corPvalueStudent(module_trait_cor, nrow(datExpr))

  # Create text matrix for heatmap (correlation and p-value)
  text_matrix <- matrix(
    paste0(signif(module_trait_cor, 2), "\n(", signif(module_trait_pval, 1), ")"),
    nrow = nrow(module_trait_cor),
    ncol = ncol(module_trait_cor)
  )

  # Create color function
  col_fun <- colorRamp2(c(-1, 0, 1), c("blue", "white", "red"))

  # Create heatmap with ComplexHeatmap
  ht <- Heatmap(
    module_trait_cor,
    name = "Correlation",
    col = col_fun,
    cluster_rows = FALSE,
    cluster_columns = FALSE,
    show_row_names = TRUE,
    show_column_names = TRUE,
    row_names_side = "left",
    column_names_side = "bottom",
    column_title = "Module-Trait Relationships",
    cell_fun = function(j, i, x, y, width, height, fill) {
      grid.text(text_matrix[i, j], x, y, gp = gpar(fontsize = 8))
    },
    heatmap_legend_param = list(
      title = "Correlation",
      at = c(-1, -0.5, 0, 0.5, 1)
    ),
    width = unit(max(8, ncol(traits) * 1.5), "cm"),
    height = unit(max(8, nrow(module_trait_cor) * 0.5), "cm")
  )

  # Save to file
  svg(output_file, width = max(8, ncol(traits) * 1.5), height = max(8, nrow(module_trait_cor) * 0.4))
  draw(ht)
  dev.off()

  cat("Saved:", output_file, "\n")

  # Return results
  results <- data.frame(
    module = rownames(module_trait_cor),
    stringsAsFactors = FALSE
  )

scripts/export_wgcna_results.R

# Export all WGCNA results

library(WGCNA)

#' Export all WGCNA results (including RDS objects for downstream skills)
#'
#' @param results Results object from run_wgcna_analysis()
#' @param output_dir Directory to save results (default: "wgcna_results")
#' @param output_prefix Prefix for output files (default: "wgcna")
export_all <- function(results, output_dir = "wgcna_results", output_prefix = "wgcna") {

  cat("\n=== Exporting WGCNA Results ===\n\n")

  # Create output directory if needed
  if (!dir.exists(output_dir)) {
    dir.create(output_dir, recursive = TRUE)
  }

  # Extract components from results
  gene_info <- results$hub_results$gene_info
  hub_genes <- results$hub_results$hub_genes
  trait_results <- results$trait_results

  # 1. Export gene-module assignments
  csv_path <- file.path(output_dir, paste0(output_prefix, "_gene_modules.csv"))
  write.csv(gene_info, csv_path, row.names = FALSE)
  cat("  Saved:", csv_path, "\n")

  # 2. Export hub genes
  hub_df <- do.call(rbind, lapply(names(hub_genes), function(mod) {
    df <- hub_genes[[mod]]
    df$module <- mod
    df
  }))
  csv_path <- file.path(output_dir, paste0(output_prefix, "_hub_genes.csv"))
  write.csv(hub_df, csv_path, row.names = FALSE)
  cat("  Saved:", csv_path, "\n")

  # 3. Export module-trait correlations
  csv_path <- file.path(output_dir, paste0(output_prefix, "_module_trait_cor.csv"))
  write.csv(trait_results$results, csv_path, row.names = FALSE)
  cat("  Saved:", csv_path, "\n")

  # 4. Export module eigengenes
  csv_path <- file.path(output_dir, paste0(output_prefix, "_eigengenes.csv"))
  write.csv(trait_results$MEs, csv_path)
  cat("  Saved:", csv_path, "\n")

  # 5. Save analysis objects as RDS for downstream skills (CRITICAL)
  cat("\n  Saving analysis objects for downstream use:\n")

  rds_path <- file.path(output_dir, paste0(output_prefix, "_network.rds"))
  saveRDS(results$network, rds_path)
  cat("    • wgcna_network.rds\n")
  cat("      (Load with: net <- readRDS('", basename(rds_path), "'))\n", sep = "")

  rds_path <- file.path(output_dir, paste0(output_prefix, "_module_colors.rds"))
  saveRDS(results$module_colors, rds_path)
  cat("    • wgcna_module_colors.rds\n")
  cat("      (Load with: colors <- readRDS('", basename(rds_path), "'))\n", sep = "")

  rds_path <- file.path(output_dir, paste0(output_prefix, "_expression_matrix.rds"))
  saveRDS(results$datExpr, rds_path)
  cat("    • wgcna_expression_matrix.rds\n")
  cat("      (Load with: expr <- readRDS('", basename(rds_path), "'))\n", sep = "")

  rds_path <- file.path(output_dir, paste0(output_prefix, "_full_results.rds"))
  saveRDS(results, rds_path)
  cat("    • wgcna_full_results.rds (complete results object)\n")
  cat("      (Load with: results <- readRDS('", basename(rds_path), "'))\n", sep = "")

  # 6. Create summary report
  cat("\n  Creating summary report...\n")
  summary_lines <- c(
    "# WGCNA Co-expression Network Analysis Summary\n",
    paste("**Total genes analyzed:**", nrow(gene_info)),
    paste("**Total samples:**", nrow(trait_results$MEs)),
    paste("**Number of modules:**", length(unique(gene_info$module)) - 1, "(excluding grey)\n"),
    "## Module Sizes"
  )

Companion files