Normalization Strategies: TMM, DESeq2 Size Factors & Quantile Normalization
Normalization Strategies: TMM, DESeq2 Size Factors & Quantile Normalization
Overcoming compositional bias from high-abundance transcripts (e.g., RuBisCO / stress proteins); executing edgeR TMM, DESeq2 Median-of-Ratios, and preprocessCore Quantile Normalization.
Key Learning Deliverables & Competencies
- Understanding compositional distortion caused by massive RuBisCO (rbcS/rbcL) transcripts
- DESeq2 Median-of-Ratios size factor calculation and mathematical derivation
- Trimmed Mean of M-values (TMM) calculation in edgeR for plant transcriptomes
- Quantile normalization implementation using preprocessCore in R
1. Compositional Bias in Plant Leaf Transcriptomes
Plant photosynthetic tissues present an extreme case of compositional library distortion. The enzyme RuBisCO (Ribulose-1,5-bisphosphate carboxylase-oxygenase) can account for up to 50% of soluble leaf protein and 15–35% of total leaf mRNA transcripts (primarily rbcS small subunit nuclear transcripts).
Under severe drought stress, plants rapidly suppress photosynthesis and downregulate rbcS to conserve energy, while sharply inducing hundreds of stress-defensive transcripts (e.g. HSP, LEA, dehydrins). If one relies merely on total sequencing depth normalization (e.g. Reads Per Million), the loss of RuBisCO reads artificially inflates the perceived expression of all other unaffected housekeeping genes!
2. Mathematical Formulations of Transcriptomic Normalization
DESeq2 Median-of-Ratios Method
For each gene $i$, calculate its geometric mean across all $m$ samples: \[ g_i = \left(\prod_{j=1}^m K_{ij}\right)^{1/m} \] Then the size factor $\hat{s}_j$ for sample $j$ is the median of ratios for non-zero genes: \[ \hat{s}_j = \text{median}_{i: g_i > 0} \left( \frac{K_{ij}}{g_i} \right) \] The median protects the size factor from being hijacked by a small handful of super-abundant or massively perturbed genes.
3. Hands-On R: Normalization & Quantile Normalization
# Load required libraries
library(DESeq2)
library(edgeR)
library(preprocessCore)
# Load raw counts matrix
counts_raw <- read.table("potato_gene_counts.txt", header = TRUE, row.names = 1, skip = 1)
counts_matrix <- as.matrix(counts_raw[, 6:ncol(counts_raw)])
colnames(counts_matrix) <- gsub(".sorted.bam", "", colnames(counts_matrix))
# 1. DESeq2 Median-of-Ratios Normalization
dds <- DESeqDataSetFromMatrix(
countData = counts_matrix,
colData = targets,
design = ~ Genotype + Drought_Stage
)
dds <- estimateSizeFactors(dds)
normalized_counts_deseq2 <- counts(dds, normalized = TRUE)
# 2. edgeR TMM (Trimmed Mean of M-values)
dge <- DGEList(counts = counts_matrix)
dge <- calcNormFactors(dge, method = "TMM")
tmm_cpm <- cpm(dge, normalized.lib.sizes = TRUE)
# 3. Quantile Normalization via preprocessCore
# Suitable for microarrays or non-negative continuous values
qn_counts <- normalize.quantiles(as.matrix(counts_matrix))
rownames(qn_counts) <- rownames(counts_matrix)
colnames(qn_counts) <- colnames(counts_matrix)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
chr1or1; divergent naming causes silent alignment drops.