Tutorial 13 — WGCNA on COVID-19 PBMCs (GSE152418)

Weighted gene co-expression networks

Author

Single Cell RNA-seq Workshop

Running the code in this tutorial

Every code chunk here is tagged with #| eval: false, so the published page shows the code without running it. To run it yourself, open the downloaded .qmd in RStudio:

Run one chunk: click the green ▶ (Run Current Chunk) at the chunk’s top-right, or press Ctrl/Cmd + Enter. This runs the chunk regardless of its eval setting — the easiest way to work through the tutorial interactively.
Run a chunk on render: change that chunk’s #| eval: false to #| eval: true (or delete the line) so it executes when the document is rendered.
Render the whole tutorial: click the Render button in RStudio (or run quarto render in a terminal) to execute the chunks top-to-bottom and knit a finished HTML. To run everything on render, set eval: true once in the YAML header at the top.

About this tutorial

A worked WGCNA pipeline with workshop-specific commentary and call-out boxes. The dataset is GEO accession GSE152418 — bulk RNA-seq of PBMCs from COVID-19 patients across a severity gradient. We’ll build a signed co-expression network and correlate the resulting modules with disease severity.

Note

Companion book chapter: Chapter 13 — WGCNA — the long-form prose treatment of this tutorial’s material, with cross-references to the prerequisite appendices.

How to use this page

The rendered HTML shows the code but does not execute it. To run it:

Download the .qmd source: Tutorial_13_WGCNA.qmd. If your browser saves the file as Tutorial_13_WGCNA.qmd.txt, drop the trailing .txt so the filename ends in .qmd, then open it in RStudio.
Download the GEO data (block below) into ../data/.
Work through the chunks, flipping eval: true as you go.

Note

Companion lecture: Lecture 13 — WGCNA

Dataset — GSE152418 (COVID-19 PBMC bulk RNA-seq)

What it is. PBMC bulk RNA-seq from a published cross-sectional study of COVID-19 patients (Arunachalam et al. 2020, Science 369: 1210–1220). Roughly 34 samples spanning healthy controls, convalescent, moderate, severe, and ICU patients.

GEO landing page: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE152418
Companion paper: Arunachalam et al. (2020) Science — https://doi.org/10.1126/science.abc6261
Files we’ll use:
- GSE152418_p20047_Study1_RawCounts.txt.gz — gene × sample raw count table
- GEO soft / matrix file for sample metadata (loaded via GEOquery)

Download the data

The tutorial reads the raw count table from ../data/. Download it once from GEO’s supplementary FTP and unzip it in place:

mkdir -p ../data && cd ../data
wget https://ftp.ncbi.nlm.nih.gov/geo/series/GSE152nnn/GSE152418/suppl/GSE152418_p20047_Study1_RawCounts.txt.gz
gunzip GSE152418_p20047_Study1_RawCounts.txt.gz

The per-sample metadata is pulled separately inside R with GEOquery::getGEO() in Step 1 — no extra download needed.

Tip

GEO link rot: if the supplementary file URL changes, get the current one from the Supplementary file table at the bottom of the GSE152418 GEO landing page.

Common errors / things that bite

WGCNA expects samples-as-ROWS — most RNA-seq workflows have samples-as-COLUMNS. Transpose before passing to pickSoftThreshold. The tutorial does this; if you adapt for your own data, watch this.

Scale-free fit R^2 never reaches 0.8 — your data is too noisy or has technical structure. Plot the soft-threshold scan (pickSoftThreshold output) — if it plateaus below 0.8, suspect batch effects or low N. Don’t force a power that gives a low fit; report the issue.

blockwiseModules runs out of memory — for >20k genes, blockwise analysis splits into ~5k-gene blocks. If RAM is the issue, set maxBlockSize = 5000 and accept the loss in cross-block module assignments.

Module-trait correlations are all p > 0.05 — most often you have too few samples (<10). WGCNA’s downstream tests need biological replication. With N = 5–8 the test is underpowered.

Solutions / instructor copy

The Think about it prompts on this page have inline answers (click “Show answer” / “Click for answer” to expand). For the full runnable instructor copy with every chunk pre-evaluated and outputs shown:

For students: the same .qmd you downloaded already contains every solution inline — just expand the collapsed answers
For instructors: render with quarto render --profile solutions from the Exercise_Folder/ directory. The fully-evaluated HTML is written to docs/Exercise_Folder/_solutions/Tutorial_13_WGCNA.html. See Exercise_Folder/_quarto-solutions.yml for the build profile

Setup

# Packages are installed in Tutorial 00 (Setup → bonus modules) — load them here.
library(WGCNA)
library(DESeq2)
library(GEOquery)
library(tidyverse)
library(CorLevelPlot)
library(gridExtra)
library(patchwork)   # combine + annotate multi-panel plots

allowWGCNAThreads()    # multi-threading

set.seed(2026)

# ---------------------------------------------------------------------------
# Output directory for this module's figures and tables. Every figure/table
# chunk below writes a file named Mod13_C<chunk>_<name> into ../output/Mod13/.
#   Mod13 = Module 13 (this tutorial);  C<n> = the nth code chunk.
# ---------------------------------------------------------------------------
out_dir <- "../output/Mod13"
dir.create(out_dir, recursive = TRUE, showWarnings = FALSE)
dir.create("../data",  showWarnings = FALSE)   # shared .rds hand-off folder (sibling of scripts/, created one level up)
message("This module writes its figures & tables to: ", normalizePath(out_dir))

# Prefer dplyr's verbs over identically-named functions from the Bioconductor
# packages loaded above: S4Vectors/IRanges/AnnotationDbi/matrixStats (pulled in by
# org.Hs.eg.db, DESeq2, miloR, SingleCellExperiment, ...) mask select, filter,
# count, desc, slice, rename, first. Re-bind the ones used below so the dplyr
# pipelines work regardless of package load order.
select <- dplyr::select; filter <- dplyr::filter; count <- dplyr::count
desc   <- dplyr::desc;   slice  <- dplyr::slice
rename <- dplyr::rename; first  <- dplyr::first

# --- Figure saving --------------------------------------------------------
# save_fig() writes every figure as a .png (for viewing) AND a .svg (vector,
# editable in Illustrator / Inkscape for a manuscript). Pass the .png path; the
# .svg is written alongside with the same basename. (.svg uses the 'svglite'
# package, installed in Tutorial 00.)
save_fig <- function(filename, plot, width, height, dpi = 300, ...) {
  ggplot2::ggsave(filename, plot, width = width, height = height, dpi = dpi, ...)
  svg_path <- paste0(tools::file_path_sans_ext(filename), ".svg")
  tryCatch(ggplot2::ggsave(svg_path, plot, width = width, height = height, ...),
           error = function(e) message("  (could not write ", basename(svg_path),
                                       " - install 'svglite'? ", conditionMessage(e), ")"))
}
# save_base_fig() does the same for BASE-graphics figures (plot()/heatmap()/etc.):
# pass a no-argument function that draws the figure; it is rendered to .png and .svg.
save_base_fig <- function(filename, draw, width = 7, height = 5, res = 300) {
  grDevices::png(filename, width = width, height = height, units = "in", res = res)
  draw(); grDevices::dev.off()
  svg_path <- paste0(tools::file_path_sans_ext(filename), ".svg")
  tryCatch({ grDevices::svg(svg_path, width = width, height = height); draw(); grDevices::dev.off() },
           error = function(e) message("  (could not write ", basename(svg_path), ")"))
}

Step 1 — Load data and metadata

data <- read.delim("../data/GSE152418_p20047_Study1_RawCounts.txt", header = TRUE)

# Pull GEO metadata
gse       <- getGEO("GSE152418", GSEMatrix = TRUE)
phenoData <- pData(phenoData(gse[[1]]))
phenoData <- phenoData[, c(1, 2, 46:50)]
head(phenoData)

# --- Table out: the trimmed sample metadata (phenoData) ---------------------
phenoData |>
  tibble::rownames_to_column("geo_accession_rowname") |>
  readr::write_csv(file.path(out_dir, "Mod13_C3_pheno_data.csv"))

# Reshape to a tidy long form, attach phenoData (matching `samples` to `title`),
# then spread back to a counts matrix indexed by Ensembl ID.
data <- data |>
  gather(key = "samples", value = "counts", -ENSEMBLID) |>
  mutate(samples = gsub("\\.", "-", samples)) |>
  inner_join(phenoData, by = c("samples" = "title")) |>
  select(1, 3, 4) |>
  spread(key = "geo_accession", value = "counts") |>
  column_to_rownames(var = "ENSEMBLID")

dim(data)
data[1:5, 1:5]

# --- Table out: counts-matrix dimensions (genes x samples) ------------------
tibble::tibble(
  stage   = "Step 1 — tidy counts matrix (Ensembl x GEO accession)",
  genes   = nrow(data),
  samples = ncol(data)
) |>
  readr::write_csv(file.path(out_dir, "Mod13_C4_counts_dimensions.csv"))

Step 2 — Outlier detection

Note

Why before normalization? Outlying samples bias the variance-stabilizing transform and inflate the dispersion estimates downstream. Drop them first.

gsg <- goodSamplesGenes(t(data))
summary(gsg)
gsg$allOK

# --- Table out: goodSamplesGenes outcome (genes/samples flagged) ------------
tibble::tibble(
  all_ok        = gsg$allOK,
  n_good_genes   = sum(gsg$goodGenes),
  n_bad_genes    = sum(!gsg$goodGenes),
  n_good_samples = sum(gsg$goodSamples),
  n_bad_samples  = sum(!gsg$goodSamples)
) |>
  readr::write_csv(file.path(out_dir, "Mod13_C5_good_samples_genes.csv"))

# Drop outlier genes
data <- data[gsg$goodGenes == TRUE, ]

htree <- hclust(dist(t(data)), method = "average")
plot(htree,
     main = "Sample clustering for outlier detection — GSE152418 PBMCs",
     sub  = "Average linkage, Euclidean distance on samples",
     xlab = "Sample")

# --- Figure out (base graphics) --------------------------------------------
save_base_fig(file.path(out_dir, "Mod13_C6_sample_dendrogram.png"), width = 10, height = 6, draw = function() {
  plot(htree,
       main = "Sample clustering for outlier detection — GSE152418 PBMCs",
       sub  = "Average linkage, Euclidean distance on samples",
       xlab = "Sample")
})

Reading the output

The x-axis lists every sample (GEO accession), and the y-axis is the Euclidean distance at which two branches merge. Samples that fuse at a very long branch length on the right side of the tree are outliers — they are far from every other sample in raw-count space. Anything that merges last, at a height visibly above the main cluster of samples, is a candidate to drop. If all samples fuse at roughly similar heights, the dataset is clean; a single isolated branch at 2–3× the height of the rest is a reliable warning sign.

pca       <- prcomp(t(data))
pca.dat   <- pca$x
pca.var   <- pca$sdev^2
pca.var.p <- round(pca.var / sum(pca.var) * 100, digits = 2)

p_pca <- ggplot(as.data.frame(pca.dat), aes(PC1, PC2)) +
  geom_point() +
  geom_text(label = rownames(pca.dat), nudge_y = 5) +
  labs(title    = "Sample PCA for outlier detection — GSE152418 PBMCs",
       subtitle = "Each point is a sample; outliers sit far from the main cloud",
       x = paste0("PC1: ", pca.var.p[1], "%"),
       y = paste0("PC2: ", pca.var.p[2], "%"),
       caption  = "Module 13 · WGCNA") +
  theme_minimal()
p_pca

# --- Figure out ------------------------------------------------------------
save_fig(file.path(out_dir, "Mod13_C7_sample_pca.png"), p_pca,
       width = 7, height = 6, dpi = 300)

Reading the output

Each point is one sample, labelled by its GEO accession. The axes show PC1 and PC2 with their percent-variance contributions; together they capture the largest two dimensions of sample-to-sample variation. Samples that sit far from the main cloud — visibly isolated along either axis — are outliers. Because PCA captures magnitude of distance (not just merge order), a sample can look borderline in the dendrogram yet be clearly outlying here; the PCA is generally more informative on a small N. Most points should form a loose ellipse; anything beyond ~3 standard deviations in either PC is a removal candidate.

Think about it

The hierarchical-clustering tree (code) and the PCA (code) both expose outliers. Which would you trust more on a small dataset?

Show answers

PCA — it gives you a continuous view of how far each sample is from the rest, in the same space the model will use. hclust only tells you the order of merging, which can hide a sample that’s mildly aberrant but never gets clustered alone.

# This example drops a couple of obvious outliers; modify to fit your eyeballed picks.
samples.to.exclude <- c("GSM4614993", "GSM4614994", "GSM4614995")
data.subset <- data[, !(colnames(data) %in% samples.to.exclude)]

# Subset metadata accordingly
colData <- phenoData |>
  filter(!(geo_accession %in% samples.to.exclude)) |>
  column_to_rownames(var = "geo_accession")

# Make sure rows of colData match cols of counts
all(rownames(colData) %in% colnames(data.subset))
all(rownames(colData) == colnames(data.subset))

Step 3 — Normalize with DESeq2 VST

WGCNA expects expression values that are roughly homoscedastic across the mean range. Variance-stabilizing transformation does that for count data:

dds <- DESeqDataSetFromMatrix(countData = data.subset,
                              colData   = colData,
                              design    = ~ 1)        # WGCNA is unsupervised

# Drop genes with < 15 reads in 3/4 of samples
dds75 <- dds[ rowSums(counts(dds) >= 15) >= 24, ]
nrow(dds75)

dds_norm <- vst(dds75)
norm.counts <- assay(dds_norm) |> t()      # WGCNA wants samples × genes

Tip

Why ~ 1? WGCNA’s network is unsupervised — the design should not encode the trait you’re going to test against, otherwise you’re “double dipping”. Use ~ 1 for an unmodelled VST.

Step 4 — Pick the soft-thresholding power

power <- c(1:10, seq(12, 30, by = 2))
sft <- pickSoftThreshold(norm.counts,
                         powerVector = power,
                         networkType = "signed",
                         verbose     = 5)
sft.data <- sft$fitIndices

a1 <- ggplot(sft.data, aes(Power, SFT.R.sq, label = Power)) +
        geom_point() + geom_text(nudge_y = 0.05) +
        geom_hline(yintercept = 0.8, color = "red") +
        labs(title = "Scale-free topology fit",
             x = "Power", y = "Scale-free topology fit, signed R²") +
        theme_classic()
a2 <- ggplot(sft.data, aes(Power, mean.k., label = Power)) +
        geom_point() + geom_text(nudge_y = 0.1) +
        labs(title = "Mean connectivity",
             x = "Power", y = "Mean connectivity") +
        theme_classic()
grid.arrange(a1, a2, nrow = 2)

# --- Figure out (multi-panel via patchwork) --------------------------------
p_sft <- (a1 / a2) +
  patchwork::plot_annotation(
    title    = "Soft-thresholding power selection — GSE152418 PBMCs",
    subtitle = "Smallest power crossing signed R² ≈ 0.8 with reasonable connectivity",
    caption  = "Module 13 · WGCNA"
  )
save_fig(file.path(out_dir, "Mod13_C10_soft_threshold.png"), p_sft,
       width = 7, height = 8, dpi = 300)

Reading the output

The two-panel plot (code) shows candidate powers on the x-axis for both panels. Top panel: the scale-free topology fit R² — the red horizontal line at 0.8 is your target; choose the smallest power that crosses it, because higher powers reduce sensitivity. Bottom panel: mean network connectivity — it should decrease steadily as power increases; if connectivity is already very low (near zero) at your chosen power, the network is too sparse. Pick the power where the top panel just clears 0.8 and the bottom panel still shows meaningful connectivity. If R² never reaches 0.8, suspect batch effects or too few samples rather than trying higher powers.

Warning

Pick the smallest power that crosses R² ≈ 0.8 and gives a reasonable mean connectivity. This worked example lands at β = 18 for this dataset; yours will vary.

Step 5 — Build the network

soft_power     <- 18
temp_cor       <- cor                    # save base::cor before WGCNA shadows it
cor            <- WGCNA::cor

bwnet <- blockwiseModules(norm.counts,
                          maxBlockSize     = 14000,
                          TOMType          = "signed",
                          power            = soft_power,
                          mergeCutHeight   = 0.25,
                          numericLabels    = FALSE,
                          randomSeed       = 1234,
                          verbose          = 3)

cor <- temp_cor                          # restore base::cor

Warning

cor shadowing. WGCNA replaces base R’s cor(). If you forget to restore it (cor <- temp_cor), unrelated downstream code can break in confusing ways.

plotDendroAndColors(bwnet$dendrograms[[1]],
                    cbind(bwnet$unmergedColors, bwnet$colors),
                    c("unmerged", "merged"),
                    dendroLabels = FALSE,
                    addGuide     = TRUE,
                    hang         = 0.03,
                    guideHang    = 0.05,
                    main         = "Gene clustering dendrogram & module colours")

# --- Figure out (base graphics) --------------------------------------------
save_base_fig(file.path(out_dir, "Mod13_C12_module_dendrogram.png"), width = 10, height = 6, draw = function() {
  plotDendroAndColors(bwnet$dendrograms[[1]],
                      cbind(bwnet$unmergedColors, bwnet$colors),
                      c("unmerged", "merged"),
                      dendroLabels = FALSE,
                      addGuide     = TRUE,
                      hang         = 0.03,
                      guideHang    = 0.05,
                      main         = "Gene clustering dendrogram & module colours")
})

Reading the output

The dendrogram (code) clusters all genes by their topological overlap; genes that are highly co-expressed sit on adjacent leaves. Below the tree, two colour bars show the module assignments: the top bar is unmerged (many fine-grained modules) and the bottom bar is merged (modules with similar eigengenes fused by mergeCutHeight = 0.25). Each distinct colour is one module; grey genes were not assigned to any module. A good result shows a modest number of discrete-coloured branches — typically 5–20 modules. Very many tiny modules (especially if unmerged and merged look similar) suggest mergeCutHeight could be raised; a single giant turquoise module may indicate the power was set too low.

Step 6 — Module eigengenes & module–trait correlations

module_eigengenes <- bwnet$MEs
head(module_eigengenes)
table(bwnet$colors)        # genes per module

# --- Tables out: module eigengenes and genes-per-module --------------------
module_eigengenes |>
  tibble::rownames_to_column("sample") |>
  readr::write_csv(file.path(out_dir, "Mod13_C13_module_eigengenes.csv"))

tibble::enframe(table(bwnet$colors),
                name = "module", value = "n_genes") |>
  dplyr::mutate(n_genes = as.integer(n_genes)) |>
  dplyr::arrange(dplyr::desc(n_genes)) |>
  readr::write_csv(file.path(out_dir, "Mod13_C13_genes_per_module.csv"))

# Build a numeric trait matrix: severity (Healthy / Convalescent / ICU / Moderate / Severe).
traits <- colData |>
  mutate(disease_state_bin = ifelse(grepl("severe", `disease state:ch1`,
                                          ignore.case = TRUE), 1, 0)) |>
  select(disease_state_bin)

severity.out <- binarizeCategoricalColumns(colData$`severity:ch1`,
                                           includePairwise = FALSE,
                                           includeLevelVsAll = TRUE,
                                           minCount = 1)
traits <- cbind(traits, severity.out)

nSamples <- nrow(norm.counts)
nGenes   <- ncol(norm.counts)

module.trait.corr   <- cor(module_eigengenes, traits, use = "p")
module.trait.corr.p <- corPvalueStudent(module.trait.corr, nSamples)

heatmap.data <- merge(module_eigengenes, traits, by = "row.names")
heatmap.data <- column_to_rownames(heatmap.data, var = "Row.names")

CorLevelPlot(heatmap.data,
             x      = names(traits),
             y      = names(module_eigengenes),
             col    = c("blue1", "skyblue", "white", "pink", "red"),
             main   = "Module–trait correlations — GSE152418 severity")

# --- Figure out (base graphics) --------------------------------------------
save_base_fig(file.path(out_dir, "Mod13_C15_module_trait_heatmap.png"), width = 8, height = 8, draw = function() {
  CorLevelPlot(heatmap.data,
               x      = names(traits),
               y      = names(module_eigengenes),
               col    = c("blue1", "skyblue", "white", "pink", "red"),
               main   = "Module–trait correlations — GSE152418 severity")
})

# --- Table out: module-trait correlations and p-values ---------------------
as.data.frame(module.trait.corr) |>
  tibble::rownames_to_column("module") |>
  readr::write_csv(file.path(out_dir, "Mod13_C15_module_trait_corr.csv"))

as.data.frame(module.trait.corr.p) |>
  tibble::rownames_to_column("module") |>
  readr::write_csv(file.path(out_dir, "Mod13_C15_module_trait_pvalues.csv"))

Reading the output

The heatmap (code) is a modules × traits matrix. Each cell shows the Pearson r between a module eigengene and a severity trait, with the colour scale running from blue (strong negative correlation) through white (no correlation) to red (strong positive correlation). Each cell also displays the correlation value and asterisks for significance (* p < 0.05, ** p < 0.01). Modules coloured strongly red or blue against the “severe” or “ICU” columns are your candidates for follow-up — their member genes co-vary with disease severity. A module that is red for “severe” and blue for “healthy” is behaving as expected for an immunity/inflammation module; any module showing conflicting associations (red for both ends of severity) may reflect a technical or confounding signal.

Think about it

Modules with strong |r| and small p against the severity column are candidates for follow-up. What’s the next thing you’d do with such a module?

Show answers

Two things in parallel:

GO / pathway enrichment of the module’s gene list (e.g. clusterProfiler::enrichGO) to give the module a biological label.
Hub gene identification: rank module members by module membership (kME) = correlation of each gene’s expression with the module eigengene. The top-kME genes are your most representative — and the most worth following up at the bench.

Step 7 — Hub genes within a module

module.gene.mapping <- as.data.frame(bwnet$colors)
module.gene.mapping |>
  rownames_to_column("gene") |>
  filter(`bwnet$colors` == "turquoise") |>     # change to your module of interest
  head(20)

# kME (gene–module correlation) — high values = hubs
mm <- as.data.frame(cor(norm.counts, module_eigengenes, use = "p"))
mm |>
  rownames_to_column("gene") |>
  arrange(desc(MEturquoise)) |>
  head(15)

# --- Tables out: gene-module assignments and top kME hub genes -------------
module.gene.mapping |>
  rownames_to_column("gene") |>
  readr::write_csv(file.path(out_dir, "Mod13_C16_gene_module_assignments.csv"))

mm |>
  rownames_to_column("gene") |>
  arrange(desc(MEturquoise)) |>
  head(15) |>
  readr::write_csv(file.path(out_dir, "Mod13_C16_turquoise_hub_genes.csv"))

Reading the output

The first table (from [code](#lst-M13-kme)) lists every gene assigned to the module you chose (here turquoise); the second table ranks all module members by their kME value — the Pearson correlation between each gene’s expression profile and the module eigengene. A kME close to 1.0 means the gene tracks the module’s average activity almost perfectly, making it a strong hub candidate. The top 5–15 hub genes are typically the most biologically interpretable members of the module and the best starting points for literature validation or follow-up wet-lab work. clusterProfiler::enrichGO() on the full gene list gives the module a pathway label.

Stretch goals

GO / KEGG enrichment per module via clusterProfiler to label modules biologically.
Validation in a held-out cohort — re-project the module eigengenes onto a second study (e.g. another COVID-19 PBMC GSE) and check whether the trait correlation reproduces.
Single-cell extension — try hdWGCNA on a scRNA-seq object (metacells per cell type, signed network).

What you’ve practiced

Outlier detection with goodSamplesGenes / hclust / PCA
VST normalization for downstream WGCNA
Picking a soft-thresholding power so the network is approximately scale-free
blockwiseModules end-to-end and reading the dendrogram
Module eigengenes, module–trait correlation, and hub genes via kME

Credits

Data: GEO GSE152418 (Arunachalam et al. 2020, Science)