Chapter 5: Expression Count Normalization, Library Sizing & Log Transformation

1. Why Raw Read Counts Cannot Be Directly Compared

Raw integer read counts output by quantification tools (HTSeq, featureCounts) represent the number of sequencing reads that aligned to each gene. However, raw counts are confounded by two massive technical biases:

• Sequencing Depth Bias (Library Size): If Sample A is sequenced to 60 million reads and Sample B is sequenced to 20 million reads, every gene in Sample A will exhibit approximately three-fold higher raw counts purely due to sequencing effort.
• Gene Length Bias: Under uniform expression, an 8,000 bp transcript (such as TTN) will generate 8 times more fragmented cDNA reads than a 1,000 bp transcript, even when both transcripts exist at identical molar concentrations in the cell.

2. Comprehensive Comparison of Normalization Units

3. Step-by-Step R Normalization Implementation

library(DESeq2)
library(vsn) # for meanSdPlot

# Construct DESeqDataSet
dds <- DESeqDataSetFromMatrix(
    countData = mock_counts,
    colData = clinical_clean,
    design = ~ Condition
)

# 1. Filter lowly expressed genes (require at least 10 reads across at least 3 samples)
keep <- rowSums(counts(dds) >= 10) >= 3
dds <- dds[keep, ]
cat("Genes retained after filtering:", sum(keep), "\n")

# 2. Compute DESeq2 Median of Ratios Size Factors
dds <- estimateSizeFactors(dds)
size_factors <- sizeFactors(dds)
cat("Sample Size Factors:\n")
print(round(size_factors, 3))

# Extract normalized counts matrix
normalized_counts <- counts(dds, normalized = TRUE)

# 3. Apply Variance Stabilizing Transformation (VST)
# VST stabilizes the variance across the entire dynamic range,
# preventing high-count genes from dominating distance metrics.
vsd <- vst(dds, blind = FALSE)
vst_matrix <- assay(vsd)

# 4. Diagnostic Mean-Variance Homoskedasticity Plot
# If properly stabilized, the red running trend line should be approximately flat.
meanSdPlot(vst_matrix)