bio-gene-regulatory-networks-coexpression-networks

What this skill does

## 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.

```r
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.

```r
# 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 d