# Bio Applied Network Modules

> Detect PPI/co-expression modules with NetworkX/python-louvain/leidenalg (Louvain, Leiden, modularity Q) and WGCNA eigengenes. Use when clustering a gene network, computing WGCNA modules, or testing DEG/pathway enrichment on network communities.

- Skill: `pavel-kravchenko/bio-applied-network-modules` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/bio-applied-network-modules`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/bio-applied-network-modules/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: Coding & Dev Tools
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/bio-applied-network-modules

---


# Network Modules: Community Detection, WGCNA Eigengenes, Enrichment

## When to Use
- Partitioning a STRING/BioGRID PPI network into functional modules (Louvain/Leiden)
- Building WGCNA-style co-expression modules from an RNA-seq expression matrix
- Testing whether a network module is enriched for a pathway/gene set (hypergeometric test)
- Overlaying DEG results onto network communities (Fisher's exact test)
- Exporting a partitioned network to Cytoscape for figure-quality visualization

## Version Compatibility
networkx ≥ 3.2, python-louvain (`community`) ≥ 0.16, leidenalg ≥ 0.10 + python-igraph ≥ 0.11, scikit-learn ≥ 1.4, scipy ≥ 1.11, pandas ≥ 2.2, Python ≥ 3.10.

## Prerequisites
- `pip install networkx python-louvain leidenalg python-igraph scikit-learn scipy pandas matplotlib seaborn py4cytoscape`
- A network already built (see `bio-database-access-interaction-databases` for STRING/BioGRID fetch) or an expression matrix (samples × genes)
- Familiarity with graph basics (`networkx`) and PCA (`sklearn.decomposition.PCA`)

## Community Detection Algorithms

| Algorithm | Complexity | Notes |
|-----------|------------|-------|
| Louvain | O(n log n) | Greedy modularity; can yield disconnected communities (known bug) |
| Leiden | O(n log n) | Fixes Louvain's disconnection bug; preferred for publication |
| Girvan-Newman | O(m²n) | Edge-betweenness removal; slow, interpretable |
| Spectral clustering | O(n³) | Principled, no resolution limit |
| Infomap | O(m) | Best for directed/flow-based networks |

**Modularity Q**: Q > 0.3 meaningful structure, Q > 0.5 strong modularity, Q ≈ 0 no better than random, Q = 1 perfect (no inter-module edges).

**Goal:** partition a network into communities and score how good the partition is.
**Approach:** try `python-louvain` first, fall back to NetworkX's built-in greedy modularity if it isn't installed; always report Q.

```python
import networkx as nx


def detect_communities(G, seed=42):
    """Partition graph G into communities and return (partition dict, modularity Q).

    partition: {node: community_id}. Falls back to networkx's greedy modularity
    maximization if python-louvain is not installed.
    """
    try:
        import community as community_louvain
        partition = community_louvain.best_partition(G, random_state=seed)
    except ImportError:
        communities = nx.community.greedy_modularity_communities(G)
        partition = {node: i for i, comm in enumerate(communities) for node in comm}

    n_comm = len(set(partition.values()))
    comm_sets = [{n for n, c in partition.items() if c == i} for i in range(n_comm)]
    Q = nx.community.modularity(G, comm_sets)
    return partition, Q


# Leiden (fixes Louvain's disconnected-community bug; needs igraph)
def detect_communities_leiden(G, seed=42):
    """Leiden partitioning via leidenalg + igraph. Returns {node: community_id}."""
    import leidenalg
    import igraph as ig

    ig_G = ig.Graph.from_networkx(G)
    part = leidenalg.find_partition(ig_G, leidenalg.ModularityVertexPartition, seed=seed)
    return {ig_G.vs[i]["_nx_name"]: cid for cid, cluster in enumerate(part) for i in cluster}


# Evaluate against a known ground-truth labeling
from sklearn.metrics import normalized_mutual_info_score, adjusted_rand_score


def score_against_truth(true_labels, detected_labels):
    """NMI and ARI between detected partition and ground-truth module labels."""
    nmi = normalized_mutual_info_score(true_labels, detected_labels)
    ari = adjusted_rand_score(true_labels, detected_labels)
    return nmi, ari
```

## WGCNA Co-expression Modules

WGCNA builds modules from an expression matrix (samples × genes), not from a PPI database.

1. Pearson correlation matrix of gene expression
2. Soft-threshold power β: raise |correlation| to β so the network approximates scale-free topology (pick β where R² > 0.85 on a log(k) vs log(P(k)) fit)
3. Topological Overlap Measure (TOM): shared-neighbor-based adjacency, more robust than raw correlation
4. Hierarchical clustering of (1 − TOM) distance, dynamic tree cut → modules
5. Module eigengene (ME) = first principal component of the module's expression

**Goal:** derive per-module summary expression (eigengenes) and inspect module-module relationships.
**Approach:** correlate → soft-threshold → PCA per module; for production use `pyWGCNA` or R `WGCNA`, this is the concept-level version.

```python
import numpy as np
from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA


def module_eigengene(expr_matrix, gene_indices):
    """First PC of a module's expression. expr_matrix: (samples x genes) array.

    gene_indices: column indices of genes belonging to one module.
    Returns a (samples,) vector - the module eigengene (ME).
    """
    mod_expr = StandardScaler().fit_transform(expr_matrix[:, gene_indices])
    return PCA(n_components=1).fit_transform(mod_expr).flatten()


def soft_threshold_adjacency(expr_matrix, beta=6):
    """Soft-thresholded adjacency = |Pearson r|^beta.

    beta should be chosen from a scale-free topology fit (R^2 > 0.85);
    this simplified version skips the full TOM (which also folds in
    shared-neighbor connectivity) - use pyWGCNA/R WGCNA TOMsimilarity() for that.
    """
    corr = np.corrcoef(expr_matrix.T)
    return np.abs(corr) ** beta
```

WGCNA in R (production-grade, includes real TOM and dynamic tree cut):

```r
library(WGCNA)
options(stringsAsFactors = FALSE)

# expr: samples x genes matrix, already normalized (e.g. VST/rlog counts)
powers <- c(1:20)
sft <- pickSoftThreshold(expr, powerVector = powers, verbose = 0)
beta <- sft$powerEstimate  # smallest power with fit R^2 > 0.85

net <- blockwiseModules(expr, power = beta, TOMType = "unsigned",
                         minModuleSize = 30, mergeCutHeight = 0.25,
                         numericLabels = TRUE, saveTOMs = FALSE)

MEs <- net$MEs                     # module eigengenes (samples x modules)
moduleColors <- labels2colors(net$colors)

# Module-trait correlation
moduleTraitCor <- cor(MEs, traitData, use = "p")
moduleTraitPvalue <- corPvalueStudent(moduleTraitCor, nrow(expr))
```

| Aspect | WGCNA | STRING/PPI |
|--------|-------|------------|
| Data source | Expression matrix | Protein interactions |
| Edge meaning | Co-expression | Physical/functional |
| Context-specific | Yes (per experiment) | No (static database) |
| Interpretation | Co-regulated genes | Complexes/partners |

## Module Enrichment and DEG Overlay

**Goal:** decide which pathway/gene set each detected module represents, and whether a module is enriched for differentially expressed genes.
**Approach:** hypergeometric test per (module, gene set) pair; Fisher's exact test for module-vs-DEG-list overlap; correct for multiple testing (BH-FDR) before calling anything significant.

```python
import numpy as np
import pandas as pd
from scipy.stats import hypergeom, fisher_exact
from statsmodels.stats.multitest import multipletests


def module_enrichment(partition, gene_sets, background_size):
    """Hypergeometric enrichment of each community against each reference gene set.

    partition: {gene: community_id}. gene_sets: {name: set(genes)}.
    Returns a DataFrame with BH-adjusted p-values (q_value).
    """
    rows = []
    n_comm = len(set(partition.values()))
    for comm_id in range(n_comm):
        comm_genes = {g for g, c in partition.items() if c == comm_id}
        n = len(comm_genes)
        for gs_name, gs_genes in gene_sets.items():
            k = len(comm_genes & gs_genes)
            if k == 0:
                continue
            K = len(gs_genes)
            pval = hypergeom.sf(k - 1, background_size, K, n)
            fold = (k / n) / (K / background_size)
            rows.append({"community": comm_id, "gene_set": gs_name,
                         "overlap": k, "fold_enrichment": fold, "p_value": pval})
    df = pd.DataFrame(rows)
    if len(df):
        df["q_value"] = multipletests(df["p_value"], method="fdr_bh")[1]
    return df


def deg_overlap_test(partition, degs, all_genes):
    """Fisher's exact test: is community enriched for DEGs vs. rest of the network?"""
    n_total = len(all_genes)
    n_comm = len(set(partition.values()))
    rows = []
    for comm_id in range(n_comm):
        comm_genes = {g for g, c in partition.items() if c == comm_id}
        a = len(degs & comm_genes)                       # DEG in module
        b = len(comm_genes) - a                           # not-DEG in module
        c = len(degs) - a                                 # DEG outside module
        d = n_total - len(comm_genes) - c                  # not-DEG outside module
        odds, pval = fisher_exact([[a, b], [c, d]], alternative="greater")
        rows.append({"community": comm_id, "module_size": len(comm_genes),
                      "n_deg_in_module": a, "odds_ratio": odds, "p_value": pval})
    return pd.DataFrame(rows)
```

## Visualization and Cytoscape Export

```python
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import networkx as nx

pos = nx.spring_layout(G, seed=42, k=3.0)  # larger k visually separates modules
colors = [cm.tab10(partition[n] / max(partition.values())) for n in G.nodes()]
nx.draw_networkx(G, pos=pos, node_color=colors, node_size=150,
                  edge_color="gray", width=0.5, font_size=5)
plt.axis("off")
```

```python
# Export for Cytoscape (requires Cytoscape Desktop running with CyREST)
import py4cytoscape as p4c
p4c.create_network_from_networkx(G)

# Or, without Cytoscape: dump edgelist + module labels for later import
import networkx as nx
nx.write_edgelist(G, "network.edgelist", data=["weight"])
pd.Series(partition, name="module").to_csv("modules.tsv", sep="\t")
```

## Pitfalls
- **Resolution limit**: Louvain/Leiden cannot detect modules smaller than ~√(2m) nodes; use hierarchical methods for small sub-communities.
- **Disconnected communities (Louvain bug)**: prefer Leiden, or verify with `nx.is_connected(G.subgraph(comm))` after partitioning.
- **WGCNA soft-threshold selection**: always check scale-free fit (R² of log(k) vs log(P(k)) > 0.85) before trusting a chosen β.
- **Module eigengene sign is arbitrary**: PCA sign can flip between runs/datasets; use absolute correlation for trait associations, not raw sign.
- **Dense networks hide modules**: pre-filter STRING edges to confidence ≥ 700 before community detection; near-complete graphs give Q ≈ 0.
- **Batch effects**: check for batch confounding before interpreting module-trait or module-DEG correlations.
- **Multiple testing**: always BH-FDR correct enrichment/overlap p-values across all (module × gene set) pairs before reporting hits.

## See Also
- `bio-database-access-interaction-databases` - fetching STRING/BioGRID PPI edges before module detection
- `bio-gene-regulatory-networks-coexpression-networks` - building the co-expression graph WGCNA modules come from
- `bio-gene-regulatory-networks-differential-networks` - comparing module structure across conditions
- `bio-pathway-analysis-go-enrichment` - richer GO/KEGG enrichment beyond the hypergeometric snippet here