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

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.

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.

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):

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))

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.

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

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")

# 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