Chapter 3: Count Filtering, Variance-Stabilizing Transformation (VST) & Sample Outlier QC

1 min read • Verified October 2026 • Peer-Reviewed Learning Resource

Tested Stacks: GATK 4.6.1 Samtools 1.21 Nextflow 24.10 Seurat v5.1 Python 3.12

Prerequisites: Basic Linux CLI & Sequencing Formats

📚 Chapter 3 of 10 0% Course Completed

Module 03 • Preprocessing & Homoscedasticity Track 23 • 10 Modular Chapters WGCNA • DESeq2 • clusterProfiler • Cytoscape

Chapter 3: Count Filtering, Variance-Stabilizing Transformation (VST) & Sample Outlier QC

Biological count filtering across sample replicates, applying DESeq2 variance-stabilizing transformation (vst) with blind=FALSE, transposing to sample-by-gene format, and running goodSamplesGenes with hierarchical outlier clustering.

⏱️ Estimated Study: 45–60 mins 🧪 Workflows: R / Bioconductor (WGCNA 1.73+, DESeq2, clusterProfiler, Cytoscape 3.10+) 📋 Level: Intermediate to Advanced Systems Transcriptomics

1. The Heteroscedasticity Challenge: Why Raw Counts Distort Networks

Standard Pearson or Spearman correlation assumes continuous, homoscedastic data where variance remains roughly constant across the dynamic range. High-throughput RNA-seq count data violates this assumption fundamentally:

Mean-Variance Dependency: In un-transformed count matrices, variance scales quadratically with mean expression ($\text{Var}(Y) = \mu + \alpha \mu^2$). Highly expressed housekeeping genes exhibit massive absolute variances that artificially dominate pairwise covariance calculations.
The Flaw of $\log_2(\text{count} + 1)$: While common, naive log transformation severely inflates the variance of low-count transcripts near zero, producing spurious correlation spikes driven by Poisson sampling noise rather than biology.
Variance-Stabilizing Transformation (VST): DESeq2's vst() calculates a parametric dispersion-mean trend, stabilizing variance across all expression tiers so that genes across all abundance classes contribute equitably to the co-expression topology.

Design Parameter: blind = FALSE vs blind = TRUE When preparing expression data for WGCNA, supply the full experimental design formula (design = ~ treatment * timepoint) and set blind = FALSE during VST. Blind dispersion estimation calculates global variance across all samples, treating genuine treatment differences as uncontrolled noise. Specifying blind = FALSE isolates technical dispersion while preserving true biological group shifts.

2. Low-Count Filtering & Homoscedastic Transformation Workflow

Before transformation, filter out unexpressed and ultra-low-count transcripts. Genes with negligible expression across all conditions add thousands of uninformative dimensions to the Topological Overlap Matrix (TOM) without biological value:

# Build DESeqDataSet with multi-factor stress design
dds <- DESeqDataSetFromMatrix(
    countData = counts,
    colData = sample_meta,
    design = ~ treatment * timepoint
)

# Rule: Retain genes with count >= 10 in at least 4 biological samples
# (Adapt threshold to your minimum replicate group size)
min_reps <- 4
keep <- rowSums(counts(dds) >= 10) >= min_reps
dds <- dds[keep, ]

cat("Genes retained after filtering:", nrow(dds), "\\n")
if (nrow(dds) < 100) {
    warning("Extremely few genes remain after filtering. Review sequencing depth and filter criteria.")
}

# Perform Variance-Stabilizing Transformation (homoscedastic scale)
vsd <- vst(dds, blind = FALSE)

# WGCNA Standard Geometry: Samples in Rows, Genes in Columns
expr_mat <- t(assay(vsd))

# Remove genes with zero or non-finite variance across retained samples
gene_var <- apply(expr_mat, 2, var, na.rm = TRUE)
expr_mat <- expr_mat[ , is.finite(gene_var) & gene_var > 0, drop = FALSE ]

# Verify sample alignment
stopifnot(identical(rownames(expr_mat), rownames(sample_meta)))
cat("Transformed expression matrix:", nrow(expr_mat), "samples ×", ncol(expr_mat), "genes.\\n")

3. Automated QC with `goodSamplesGenes` & Outlier Dendrograms

WGCNA provides an automated iterative diagnostic, goodSamplesGenes(), to scan the expression matrix for genes or biological samples with excessive missing values or invariant expression trajectories:

# Execute iterative quality check
gsg <- goodSamplesGenes(expr_mat, verbose = 3)

if (!gsg$allOK) {
    cat("Trimming problematic samples and invariant genes...\\n")
    expr_mat <- expr_mat[ gsg$goodSamples, gsg$goodGenes, drop = FALSE ]
    sample_meta <- sample_meta[rownames(expr_mat), , drop = FALSE]
}
stopifnot(identical(rownames(expr_mat), rownames(sample_meta)))

# Hierarchical Sample Clustering to Detect Batch/Degradation Outliers
sample_tree <- hclust(dist(expr_mat), method = "average")

par(mar = c(0, 4, 2, 0))
plot(
    sample_tree,
    main = "Sample Clustering to Detect Biological Outliers",
    sub = "",
    xlab = "",
    cex.lab = 1.2,
    cex.axis = 1.0,
    cex.main = 1.3
)
# Optional: Add cut line if an extreme outlier sample is visually separated
# abline(h = 120, col = "red", lty = 2)

Following sample verification, prepare the numeric trait matrix. Continuous physiological variables (e.g., Relative Water Content, stomatal conductance) or binary indicator vectors (e.g., is_drought) must be aligned with rownames(expr_mat):

# Construct numeric trait matrix for correlation testing
traits <- sample_meta

# Encode binary stress indicator
traits$is_drought <- as.numeric(traits$treatment == "Drought")

# Incorporate real physiological metrics (e.g., Relative Water Content %)
# if ( "RWC" %in% colnames(sample_meta) ) traits$RWC <- as.numeric(sample_meta$RWC)

trait_numeric <- data.frame(
    is_drought = traits$is_drought,
    row.names = rownames(traits)
)
stopifnot(identical(rownames(expr_mat), rownames(trait_numeric)))
cat("✓ Trait matrix successfully synchronized with expression matrix.\\n")

← Previous: Ch 2 📚 All 10 Chapters Next: Ch 4 →

⚠️

Common Pitfalls & Troubleshooting Checklist

Java Heap Exhaustion (GATK / Picard): If processes terminate with java.lang.OutOfMemoryError, set your environment variable explicitly: export _JAVA_OPTIONS="-Xmx16g -Djava.io.tmpdir=/tmp".
Index Timestamp Mismatch: If tools complain that a BAM index is missing or out of date, regenerate with samtools index -@ 4 sample.bam.
Conda Channel Conflicts: Avoid ABI incompatibilities between bioconda and conda-forge by setting strict priority: conda config --set channel_priority strict.
Chromosome Naming Contamination: Always verify whether reference FASTA headers use chr1 or 1; divergent naming causes silent alignment drops.