GPTomics/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

| Method | Edge type | Strength | Use / fails when | |--------|-----------|----------|------------------| | WGCNA (Langfelder & Horvath 2008) | marginal, soft-threshold + TOM | de facto standard; eigengenes; preservation | bulk, n>=15-20; fails on raw scRNA-seq | | hdWGCNA (Morabito 2023) | marginal, on metacells | flattens dropout/zero-inflation | single-cell; fails if metacells overlap (pseudo-replication) | | MEGENA (Song & Zhang 2015) | marginal, planar filtered network | principled sparsification; multiscale/nested modules | hierarchical biology; nested output breaks WGCNA-style one-module-per-gene code | | CEMiTool (Russo 2018) | marginal, auto-beta | fast first pass + GSEA/ORA + HTML report | exploratory; S4 object (accessors, not $); silently variance-filters genes | | GeneNet / graphical lasso | partial (direct edges) | distinguishes direct from indirect | when "who regulates whom directly" matters; collapses at small n / large p |

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

| Scenario | Recommended | Why | |----------|-------------|-----| | Bulk RNA-seq, >=15-20 samples, want modules + trait links | WGCNA signed, bicor | stable marginal modules; eigengene-trait correlation | | Single-cell / sparse counts | hdWGCNA on metacells | naive WGCNA fails on zero-inflation; metacells flatten dropout | | Need direct vs indirect edges | GeneNet / graphical lasso | partial correlation removes confounded indirect edges | | Hierarchical / multiscale structure expected | MEGENA | nested modules at multiple resolutions | | Want directed TF -> target regulons | -> scenic-regulons (single-cell) / grn-inference (bulk) | co-expression is undirected; motif/regression priors add direction | | Compare networks across conditions | -> differential-networks | rewiring is a different question from module discovery | | <15 samples | none reliable | correlation estimates too noisy; report this, do not force WGCNA |

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.

Per-Method Failure Modes

Unsigned-by-default network

Trigger: leaving networkType='unsigned' (the default) then interpreting modules as co-regulated programs. Mechanism: unsigned uses |cor|, so a gene and its strong anti-correlate land in the same module. Symptom: modules mixing clearly opposing programs (e.g. proliferation and quiescence). Fix: use networkType='signed' (and TOMType='signed') at pickSoftThreshold AND blockwiseModules.

Soft power chosen on the wrong network type

Trigger: pickSoftThreshold() run unsigned, then blockwiseModules(networkType='signed'). Mechanism: the scale-free fit and chosen power differ by network type. Symptom: near-disconnected network or one giant module. Fix: set networkType identically in both calls (signed power is roughly double unsigned).

blockwiseModules block artifact

Trigger: n_genes > maxBlockSize (historical default 5000) with no mention of blocking. Mechanism: TOM is computed within each block only; cross-block co-expression is invisible and assignments shift with block size/RAM. Symptom: module membership changes when rerun on a different machine. Fix: set maxBlockSize >= n_genes, or report the block size.

Naive WGCNA on raw scRNA-seq

Trigger: running WGCNA on sparse single-cell counts. Mechanism: zero-inflation creates a spike of zero correlations and dropout-driven spurious correlation. Symptom: a pathological correlation distribution and unstable modules. Fix: use hdWGCNA metacells (and cap max_shared).

Detection without preservation

Trigger: reporting modules from one cohort with no replication test. Mechanism: any dendrogram can be cut into modules. Symptom: modules that do not reproduce in independent data. Fix: run modulePreservation(); report Zsummary/medianRank.

Quantitative Thresholds

| Threshold | Source | Rationale | |-----------|--------|-----------| | Samples >= 15 (ideally >= 20) | Horvath/Langfelder WGCNA FAQ | correlation estimates too noisy below 15; modules become random | | Scale-free R^2 >= 0.8 (0.8-0.9) | WGCNA convention | heuristic for power where connectivity flattens -- NOT proof of scale-free biology | | Signed soft power ~12 (vs ~6 unsigned) | WGCNA FAQ table | signed adjacency = ((1+cor)/2)^power needs higher power for the same connectivity | | maxPOutliers = 0.05-0.10 | Langfelder & Horvath 2012 | prevents bicor flagging legitimate observations as outliers at modest n | | minModuleSize = 30 (default); mergeCutHeight = 0.25 (chosen; default is 0.15) | WGCNA | smaller modules are usually noise; 0.25 merges eigengenes correlating > 0.75 | | Zsummary > 10 / 2-10 / < 2 | Langfelder 2011 | strong / weak-moderate / no preservation evidence |

Common Errors

| Error / symptom | Cause | Solution | |-----------------|-------|----------| | One giant module + large grey | power too low or wrong network type | re-pick power on the correct networkType | | cemitool results via $ return NULL | CEMiTool returns an S4 object | use accessors module_genes(), get_hubs() | | hub list dominated by housekeeping genes | hubness tracking mean expression | define hubs by kME, control for expression level | | top module is the sequencing batch | batch correlates with many genes | correct batch before network construction (-> differential-expression/batch-correction) | | modules irreproducible across reruns | block artifact | set maxBlockSize >= n_genes |

References

• Langfelder P, Horvath S. 2008. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 9:559.
• Langfelder P, Horvath S. 2012. Fast R functions for robust correlations and hierarchical clustering. J Stat Softw 46(11):1-17. -- bicor / maxPOutliers.
• Langfelder P, Luo R, Oldham MC, Horvath S. 2011. Is my network module preserved and reproducible? PLoS Comput Biol 7(1):e1001057.
• Broido AD, Clauset A. 2019. Scale-free networks are rare. Nat Commun 10:1017. -- soft-threshold caveat.
• Khanin R, Wit E. 2006. How scale-free are biological networks? J Comput Biol 13(3):810-818.
• Song WM, Zhang B. 2015. Multiscale embedded gene co-expression network analysis (MEGENA). PLoS Comput Biol 11(11):e1004574.
• Russo PST, et al. 2018. CEMiTool. BMC Bioinformatics 19:56.
• Morabito S, Reese F, Rahimzadeh N, Miyoshi E, Swarup V. 2023. hdWGCNA. Cell Rep Methods 3(6):100498.

Related Skills

• differential-networks - compare co-expression structure between conditions (rewiring)
• scenic-regulons - directed TF regulons from single-cell data (motif-pruned)
• grn-inference - bulk directed GRN inference and TF protein-activity (VIPER)
• differential-expression/batch-correction - remove batch effects before network construction
• single-cell/preprocessing - QC and normalization for single-cell inputs to hdWGCNA
• pathway-analysis/go-enrichment - functional enrichment of co-expression modules