Chapter 6: Exploratory Data Analysis & Dimensionality Reduction: PCA, MDS & t-SNE

1. High-Dimensional Transcriptomic Landscapes

An RNA-Seq experiment typically measures 20,000+ protein-coding genes simultaneously across a cohort of samples. Visualizing this high-dimensional space requires dimensionality reduction techniques that project maximum biological variance onto low-dimensional orthogonal coordinates (PC1, PC2). Exploratory Data Analysis (EDA) answers three critical diagnostic questions:

1. Do biological replicate samples cluster tightly together by condition?
2. Are there extreme sample outliers caused by degradation, sequencing failure, or mislabeling?
3. Is the primary axis of variance driven by the biological variable of interest, or by an uncorrected technical batch effect?

2. Dimensionality Reduction Comparison: PCA vs. MDS vs. t-SNE

3. Production PCA & Scree Plot Script

library(DESeq2)
library(ggplot2)

# Extract top 500 most variable genes from VST matrix
ntop <- 500
rv <- rowVars(assay(vsd))
select <- order(rv, decreasing = TRUE)[seq_len(min(ntop, length(rv)))]
pca <- prcomp(t(assay(vsd)[select, ]))

# Calculate percentage variance explained by each PC
percent_var <- round((pca$sdev^2 / sum(pca$sdev^2)) * 100, 1)

# Construct PCA Dataframe
pca_df <- data.frame(
    PC1 = pca$x[, 1],
    PC2 = pca$x[, 2],
    PC3 = pca$x[, 3],
    Sample = colnames(vsd),
    Condition = colData(vsd)$Condition,
    Batch = colData(vsd)$Extraction_Batch
)

# 1. Publication-Grade PCA Biplot
ggplot(pca_df, aes(x = PC1, y = PC2, color = Condition, shape = Batch)) +
    geom_point(size = 4, alpha = 0.9) +
    stat_ellipse(aes(group = Condition, fill = Condition), geom = "polygon", alpha = 0.1, linetype = "dashed") +
    scale_color_manual(values = c("Control" = "#06b6d4", "Treated" = "#f43f5e")) +
    scale_fill_manual(values = c("Control" = "#06b6d4", "Treated" = "#f43f5e")) +
    labs(
        title = "Principal Component Analysis (PCA) of Transcriptomic Cohort",
        subtitle = "Clear biological separation along PC1; Batch variance checked along PC2",
        x = paste0("PC1: ", percent_var[1], "% variance explained"),
        y = paste0("PC2: ", percent_var[2], "% variance explained")
    ) +
    theme_bw(base_size = 14) +
    theme(panel.grid.minor = element_blank())

# 2. Scree Plot of Eigenvalue Variances
scree_df <- data.frame(
    PC = paste0("PC", 1:10),
    Variance = percent_var[1:10]
)
scree_df$PC <- factor(scree_df$PC, levels = scree_df$PC)

ggplot(scree_df, aes(x = PC, y = Variance)) +
    geom_col(fill = "#3b82f6", width = 0.6) +
    geom_line(aes(group = 1), color = "#f59e0b", linewidth = 1) +
    geom_point(color = "#f59e0b", size = 3) +
    labs(
        title = "PCA Scree Plot (Percentage Variance Explained)",
        x = "Principal Component",
        y = "% Variance Explained"
    ) +
    theme_minimal(base_size = 14)