Bio Gene Regulatory Networks Coexpression Networks

Version Compatibility

Reference examples tested with: WGCNA 1.72+, hdWGCNA 0.3+, CEMiTool 1.26+, GENIE3-adjacent GGM via GeneNet 1.2.16+.

Before using code patterns, verify installed versions match. If versions differ:

• R: packageVersion('<pkg>') then ?function_name to verify parameters
• Python: pip show <package> then help(module.function) to check signatures

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

WGCNA argument defaults differ by entry point: pickSoftThreshold() and blockwiseModules() default to networkType='unsigned'. The signed-network choice (below) must be set identically at every step or the soft power and modules silently mismatch.

Co-expression Networks

"Find co-expression modules and hub genes from my expression data" -> Build a weighted gene co-expression network, detect modules of co-regulated genes by clustering a topological-overlap dissimilarity, summarize each module by its eigengene, and relate modules to sample traits.

• R: WGCNA::blockwiseModules() for network construction + module detection (bulk)
• R: hdWGCNA metacell workflow for single-cell expression
• R: GeneNet/graphical lasso when direct (not indirect) edges are required

The Single Most Important Modern Insight -- Co-expression Is Marginal Correlation; Regulation Is Conditional Dependence

A WGCNA edge between genes B and C exists whenever a common driver A correlates with both -- even if B and C never interact. If A regulates B and C, then B and C are conditionally independent given A, yet a co-expression network still draws a strong B-C edge. A co-expression module is a descriptive object ("these genes vary together"), not a regulatory network, even when its hub is a transcription factor. The dividing line is marginal vs partial correlation: WGCNA, MEGENA, and CEMiTool measure marginal correlation (direct + indirect edges mixed together), while Gaussian graphical models (GeneNet, graphical lasso) estimate partial correlation (conditional dependence, i.e. direct edges only). Decide which object the biological question needs before choosing a tool. Causal/regulatory language ("module X is driven by hub TF Y") from a marginal co-expression edge is the single most common over-claim in this domain.

A secondary, load-bearing caveat: the scale-free topology assumption underpinning soft-threshold selection is empirically weak. Broido & Clauset 2019 (Nat Commun 10:1017) found scale-free structure in only ~4% of ~1000 real networks; Khanin & Wit 2006 (J Comput Biol 13:810) found none of 10 biological networks fit a pure power law. The scale-free R^2 >= 0.8 rule is a heuristic for picking the power where the connectivity curve flattens, not a hypothesis test confirming the biology is scale-free -- so it does not license "the network is scale-free, therefore hubs are master regulators."

Co-expression Method Taxonomy

Rule of thumb: module discovery and trait association -> WGCNA (signed); single-cell -> hdWGCNA; "is this edge direct or indirect?" -> a Gaussian graphical model.

Decision Tree by Scenario

WGCNA: Build a Signed Network

Goal: Detect robust co-expression modules from a bulk expression matrix and choose network parameters defensibly.

Approach: Use a signed network with biweight midcorrelation (bicor) so activators and repressors are not merged into one module and outliers are down-weighted; pick the soft power on the SAME network type used for construction; set maxBlockSize above the gene count to avoid the block artifact.

library(WGCNA)
options(stringsAsFactors = FALSE)
allowWGCNAThreads()

# WGCNA convention: genes as columns, samples as rows
expr <- t(as.matrix(read.csv('normalized_counts.csv', row.names = 1)))
gsg <- goodSamplesGenes(expr, verbose = 0)
expr <- expr[gsg$goodSamples, gsg$goodGenes]

# Soft power on the SIGNED fit (must match construction below)
powers <- 1:20
sft <- pickSoftThreshold(expr, powerVector = powers, networkType = 'signed', verbose = 0)
soft_power <- sft$powerEstimate            # signed networks usually land near 12

# Signed network, bicor with maxPOutliers to avoid spurious outlier flags at modest n.
# maxBlockSize >= n_genes keeps everything in one block (no cross-block blindness).
net <- blockwiseModules(
    expr, power = soft_power,
    networkType = 'signed', TOMType = 'signed',
    corType = 'bicor', maxPOutliers = 0.05,
    minModuleSize = 30, mergeCutHeight = 0.25, deepSplit = 2,
    maxBlockSize = ncol(expr) + 1,
    numericLabels = TRUE, pamRespectsDendro = FALSE, verbose = 0
)
module_colors <- labels2colors(net$colors)
table(module_colors)                       # large 'grey' fraction = poor fit, not a module

WGCNA: Eigengenes, Hubs, and Trait Relationships

Goal: Relate modules to sample traits and identify defensible hub genes.

Approach: Summarize each module by its eigengene (first PC), correlate eigengenes with traits, and define hubs by module membership kME (signed, bounded, comparable across modules) rather than raw connectivity.

# Recompute eigengenes from the COLOR labels so ME/kME columns are MEturquoise/kMEturquoise.
# (net$MEs uses numeric names ME0/ME1... under numericLabels=TRUE and won't match colors.)
MEs <- orderMEs(moduleEigengenes(expr, module_colors)$eigengenes)
traits <- read.csv('sample_traits.csv', row.names = 1)

module_trait_cor  <- cor(MEs, traits, use = 'p')
module_trait_pval <- corPvalueStudent(module_trait_cor, nrow(expr))

# Hub = high module membership (kME), the WGCNA-preferred definition.
# kME is signed and bounded [-1,1] with a p-value; prefer it over kWithin connectivity.
kME <- signedKME(expr, MEs)
moi <- 'turquoise'
moi_genes <- colnames(expr)[module_colors == moi]
hubs <- sort(kME[moi_genes, paste0('kME', moi)], decreasing = TRUE)
head(hubs, 20)

WGCNA: Module Preservation (the step that makes a module real)

Goal: Test whether discovered modules reproduce in an independent dataset -- the scientific claim, not module detection itself.

Approach: Run modulePreservation() with the discovery network as reference and an independent cohort as test; read Zsummary and medianRank.

# multiData/multiColor hold reference + test expression and the reference module labels
mp <- modulePreservation(
    multiData, multiColor,
    referenceNetworks = 1, nPermutations = 200,
    randomSeed = 1, verbose = 0
)
stats <- mp$preservation$Z$ref.reference$inColumnsAlsoPresentIn.test
stats[, c('moduleSize', 'Zsummary.pres')]
# Zsummary > 10 strong preservation; 2-10 weak/moderate; < 2 no evidence (Langfelder 2011).
# medianRank (mp$preservation$observed) compares modules to each other, size-independent.

hdWGCNA: Single-Cell Co-expression

Goal: Build co-expression networks from scRNA-seq without dropout-driven spurious correlation.

Approach: Aggregate transcriptionally similar cells into metacells (pseudobulk) to flatten zero-inflation; cap metacell sharing to avoid pseudo-replication; then run standard WGCNA on the metacells.

library(hdWGCNA); library(Seurat)
seurat_obj <- readRDS('clustered.rds')
seurat_obj <- SetupForWGCNA(seurat_obj, gene_select = 'fraction', fraction = 0.05,
                            wgcna_name = 'hdwgcna')

# Build metacells WITHIN cell types; max_shared caps how many cells two metacells share.
# Heavy overlap inflates downstream correlations (pseudo-replication) -- keep it low.
seurat_obj <- MetacellsByGroups(seurat_obj, group.by = 'cell_type',
                                k = 25, max_shared = 10, ident.group = 'cell_type')
seurat_obj <- NormalizeMetacells(seurat_obj)

seurat_obj <- SetDatExpr(seurat_obj, group.by = 'cell_type', group_name = 'all')
seurat_obj <- TestSoftPowers(seurat_obj, networkType = 'signed')
seurat_obj <- ConstructNetwork(seurat_obj, setDatExpr = FALSE, tom_name = 'hdwgcna')
seurat_obj <- ModuleEigengenes(seurat_obj)
seurat_obj <- ModuleConnectivity(seurat_obj)   # kME per module

Direct vs Indirect Edges: Gaussian Graphical Model

Goal: Recover edges that survive conditioning on all other genes (direct dependence), not marginal co-expression.

Approach: Estimate a shrinkage partial-correlation matrix (n << p safe) and test edges by local FDR.

library(GeneNet)
# expr: samples x genes
pcor   <- ggm.estimate.pcor(expr)              # shrinkage partial correlations
edges  <- network.test.edges(pcor, direct = FALSE, plot = FALSE)
net_gg <- extract.network(edges, method.ggm = 'prob', cutoff.ggm = 0.9)
# Far sparser than WGCNA -- that is the point: indirect edges have been removed.

CEMiTool (cemitool(expr, annot)) returns an S4 object -- use accessors (module_genes(), get_hubs(), generate_report()), never $ -- and silently variance-filters genes with filter=TRUE. MEGENA returns nested modules; do not feed them to one-module-per-gene WGCNA code.