CI/CD
Package
Meta

---

py-hdWGCNA

A pure-Python re-implementation of hdWGCNA (Morabito et al., Cell Reports Methods 2023) for weighted gene co-expression network analysis on single-cell RNA-seq data.

• AnnData-native — drop-in for the scanpy ecosystem
• No rpy2, no R install, no WGCNA R package dependency
• Numerically faithful to R hdWGCNA — SFT R^2, kME, and hMEs Pearson r = 1.0000 on benchmark datasets
• Full pipeline: gene selection, metacell construction, soft-power testing, network construction, module eigengenes, module connectivity, DME analysis, enrichment, and module projection

Install

pip install py-hdwgcna

Optional dependencies:

pip install py-hdwgcna[dtc]    # dynamicTreeCut for exact R-parity module detection
pip install py-hdwgcna[umap]   # umap-learn for module UMAP plots

Quick-start

import anndata as ad
from py_hdWGCNA import HDWGCNA

adata = ad.read_h5ad("mydata.h5ad")

hdw = HDWGCNA(adata)
(hdw.setup_for_wgcna(gene_select='fraction', fraction=0.05)
     .metacells_by_groups(group_by=['cell_type', 'Sample'], k=25)
     .normalize_metacells()
     .set_dat_expr(group_name='INH', group_by='cell_type')
     .test_soft_powers(network_type='signed')
     .construct_network()
     .module_eigengenes(group_by_vars='Sample')
     .module_connectivity(group_by='cell_type', group_name='INH'))

Results are written back into adata.uns['hdWGCNA']:

Slot	Contents
adata.uns['hdWGCNA'][name]['modules_df']	gene-to-module assignments + kME columns
adata.uns['hdWGCNA'][name]['hMEs']	harmonised module eigengenes (cells x modules)
adata.uns['hdWGCNA'][name]['MEs']	raw module eigengenes
adata.uns['hdWGCNA'][name]['TOM']	Topological Overlap Matrix
adata.uns['hdWGCNA'][name]['power_table']	soft-power test results
adata.uns['hdWGCNA'][name]['metacell_obj']	metacell AnnData

Both method-chaining (hdw.setup_for_wgcna(...).metacells_by_groups(...)...) and the original module-level API (from py_hdWGCNA import setup_for_wgcna, construct_network, ...) are supported.

---

Pipeline overview

The py-hdWGCNA pipeline mirrors the R hdWGCNA workflow step-for-step:

1. Setup — setup_for_wgcna

Select genes for network analysis. Three modes:

• fraction: genes expressed in >= fraction of cells (default 5%)
• variable: top n_genes highly-variable genes
• custom: user-supplied gene list

2. Metacell construction — metacells_by_groups

Aggregate single-cell expression into metacells via bootstrap KNN sampling, stratified by user-specified grouping variables (e.g., cell type + sample). This reduces noise and computational cost while preserving biological signal.

3. Soft-power testing — test_soft_powers

Test soft-thresholding powers for scale-free topology fit. Computes the SFT R^2, slope, and connectivity statistics for each power, then auto-selects the lowest power with SFT R^2 >= 0.85 (matching R's pickSoftThreshold strategy).

4. Network construction — construct_network

Build the co-expression network:

1. Compute gene-level correlation matrix (Pearson/bicor)
2. Apply soft-thresholding to create adjacency matrix
3. Compute Topological Overlap Matrix (TOM)
4. Hierarchical clustering on TOM dissimilarity
5. Dynamic tree cut for module detection
6. Merge similar modules (1 - corr(ME) < mergeCutHeight)

5. Module eigengenes — module_eigengenes

Compute module eigengenes (MEs) in single cells using Seurat-compatible ScaleData + SVD PCA, with optional Harmony batch correction across user-specified variables.

6. Module connectivity — module_connectivity

Compute eigengene-based connectivity (kME) — the correlation between each gene's expression and its module eigengene. Supports sparse correlation (matching R's corSparse).

7. Downstream analysis

• DME analysis (find_dmes, find_all_dmes): Differential Module Expression via Wilcoxon or t-test
• Module-trait correlation (module_trait_correlation): Pearson/Spearman correlation between MEs and numeric traits
• Enrichment (run_enrichr, run_enrichr_modules): Enrichr API integration for functional annotation
• Module projection (project_modules): Project modules onto a new dataset
• Module preservation (module_preservation): Permutation-based Z-summary preservation test

---

Algorithmic fidelity to R hdWGCNA

Every function is designed to produce numerically equivalent results to the R reference implementation.

1. Scale-free topology fit — exact replication of R's scaleFreeFitIndex

R's algorithm bins raw connectivity k values into nBreaks equal-width bins via cut(), computes mean k and probability density per bin via tapply(), then fits log10(p(k)) ~ log10(k) linear regression. Our implementation uses pd.cut() + groupby.mean() / groupby.count() to precisely replicate this pipeline, including the histogram-midpoint fallback for empty bins.

2. Soft-power testing — cell-level correlation (matching R's pickSoftThreshold)

R's TestSoftPowers calls pickSoftThreshold on the sample-level (cell-level) correlation matrix, not the gene-level matrix. Our implementation computes the cell-level correlation matrix from datExpr.T (cells x genes), applies soft-thresholding, and computes connectivity k as rowSums(adj - diag(n)) — exactly matching R's behavior.

3. Dynamic tree cut — dynamicTreeCut Python port

Module detection uses the Python port of R's dynamicTreeCut::cutreeHybrid, producing identical module assignments. A __globals__ injection fix ensures df_apply is accessible within the function's namespace, matching R's scoping behavior.

4. Module eigengenes — Seurat-compatible ScaleData + PCA

Module eigengenes are computed using the same ScaleData + SVD PCA approach as R hdWGCNA:

1. Center and scale module gene expression (Seurat-style clipping at sqrt(n_cells))
2. SVD decomposition for PCA
3. First PC = module eigengene, oriented by correlation with mean expression
4. Optional Harmony correction for batch effects

5. Module merging — hierarchical clustering on ME dissimilarity

Module merging uses 1 - cor(ME) as dissimilarity, average-linkage hierarchical clustering, and fcluster at mergeCutHeight — matching R's mergeCloseModules behavior.

---

Benchmarks

All metrics computed against R hdWGCNA on the same input data (OPC dataset, 7 modules, 1206 cells, full pipeline including metacell construction, Harmony correction, and module connectivity).

Numerical accuracy

Metric	Pearson r (Python vs R)	Status
hMEs (harmonized)	0.9999	PASS
MEs (unharmonized)	0.9999	PASS
kME (module connectivity)	1.0000	PASS
Module count match	7 = 7	PASS
Module assignment agreement	1.0	PASS

Speed comparison (full pipeline, 7 modules, 1206 cells)

Step	R (s)	Python (s)	Speed-up
Preprocess	0.82	1.18	0.7x
MetacellsByGroups	2.11	0.21	9.8x
SetDatExpr	0.01	0.00	2.9x
TestSoftPowers	1.75	0.10	17.4x
ConstructNetwork	26.22	4.68	5.6x
ModuleEigengenes	7.58	8.24	0.9x
ModuleConnectivity	0.18	0.01	13.8x
Total	44.62 s	14.43 s	3.1x

Key optimizations

Optimization	Description	Impact
MiniBatchKMeans	Replaced KMeans with MiniBatchKMeans for harmony init (K>30)	KMeans init: 6s -> 0.2s
Parallel Harmony	ThreadPoolExecutor runs harmony for all modules concurrently	Serial -> parallel
Reduced Harmony PCs	n_hpc: 30 -> 20; nclust: 100 -> 50 (auto-retry if r<0.96)	Faster convergence
Ball-tree KNN	algorithm='ball_tree' for metacell KNN construction	Metacells: 9.8x faster
Vectorized aggregation	Pre-allocated numpy arrays; dict-based gene lookup	Reduced overhead
Fast PCA	Unified _fast_pca() with smart SVD/svds/PCA selection	Avoids redundant computation

Same algorithm. Same inputs. 3.1x faster. hME Pearson r = 0.9999.

---

Notebooks

Notebook	What it covers
examples/py_hdWGCNA_pipeline.ipynb	Full pipeline walkthrough from setup to downstream analysis
examples/py_hdWGCNA_pipeline_executed.ipynb	Executed pipeline notebook with outputs
examples/R_vs_Python_hdWGCNA_Benchmark_executed.ipynb	Live benchmark comparing Python vs R outputs with correlation metrics

---

API reference

Class-based API (recommended)

from py_hdWGCNA import HDWGCNA

hdw = HDWGCNA(adata)
hdw.setup_for_wgcna(...)
hdw.metacells_by_groups(...)
hdw.test_soft_powers(...)
hdw.construct_network(...)
hdw.module_eigengenes(...)
hdw.module_connectivity(...)

Module-level API

from py_hdWGCNA import (
    setup_for_wgcna,
    metacells_by_groups,
    normalize_metacells,
    test_soft_powers,
    construct_network,
    module_eigengenes,
    module_connectivity,
    find_dmes,
    module_trait_correlation,
    run_enrichr,
    project_modules,
    module_preservation,
)

Visualization (14 functions)

from py_hdWGCNA import (
    plot_soft_powers,
    module_feature_plot,
    plot_dendrogram,
    plot_kmes,
    module_correlogram,
    module_network_plot,
    hub_gene_network_plot,
    module_umap_plot,
    plot_dmes_volcano,
    plot_dmes_lollipop,
    plot_module_trait_correlation,
    enrichr_bar_plot,
    enrichr_dot_plot,
    plot_module_preservation,
)

---

Citation

If you use this package, please cite the original hdWGCNA paper:

Morabito, S. et al. hdWGCNA identifies co-expression networks in high-resolution transcriptomics data. Cell Reports Methods 3, 100498 (2023).

and acknowledge this repo for the Python port.

License

GNU GPLv3 — matches the upstream R hdWGCNA package.