Introduction

The dataset used here is from the official DESeq2 tutorial and countless RNA-seq workshops, run here with actual raw counts and DESeq2.

Dataset: Himes et al. 2014, PLOS ONE, “RNA-Seq Transcriptome Profiling Identifies CRISPLD2 as a Glucocorticoid Responsive Gene that Modulates Cytokine Function in Airway Smooth Muscle Cells” — GEO accession GSE52778. Four human airway smooth muscle cell lines, each treated with dexamethasone (a synthetic glucocorticoid used to treat asthma) or left untreated — 8 samples, paired by cell line. This is the same experiment behind Bioconductor’s airway package that DESeq2’s own vignette uses as its primary example.

Where this data came from: GEO’s supplementary files for this series only include FPKM — but the airway package’s raw gene-level counts have also been re-published as plain CSV files by the bioconnector.org teaching materials (airway_rawcounts.csv, airway_metadata.csv), which is what this notebook loads. Same 8 samples, same SRA run accessions, same GEO sample IDs — just as raw integer counts instead of FPKM, which is what DESeq2 actually needs.

(If you already have the airway Bioconductor package installed, there’s an even shorter path — please see the note at the very end.)

Step 1 — Load the data

counts: rows = genes (Ensembl IDs), columns = samples, raw integer read counts — the input DESeq2 expects, and the natural orientation of the CSV, so no transposing needed. meta: which of the 4 cell lines each sample came from, and whether it was treated with dexamethasone.

counts <- as.matrix(read.csv("airway_rawcounts.csv", row.names = 1, check.names = FALSE))
meta   <- read.csv("airway_metadata.csv", row.names = 1)
meta   <- meta[colnames(counts), , drop = FALSE]

cat(sprintf("Count matrix: %d genes x %d samples\n", nrow(counts), ncol(counts)))
## Count matrix: 64102 genes x 8 samples
meta
##                dex celltype     geo_id
## SRR1039508 control   N61311 GSM1275862
## SRR1039509 treated   N61311 GSM1275863
## SRR1039512 control  N052611 GSM1275866
## SRR1039513 treated  N052611 GSM1275867
## SRR1039516 control  N080611 GSM1275870
## SRR1039517 treated  N080611 GSM1275871
## SRR1039520 control  N061011 GSM1275874
## SRR1039521 treated  N061011 GSM1275875
#looking at the top 6 count rows
counts[1:6, ]
##                 SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516
## ENSG00000000003        679        448        873        408       1138
## ENSG00000000005          0          0          0          0          0
## ENSG00000000419        467        515        621        365        587
## ENSG00000000457        260        211        263        164        245
## ENSG00000000460         60         55         40         35         78
## ENSG00000000938          0          0          2          0          1
##                 SRR1039517 SRR1039520 SRR1039521
## ENSG00000000003       1047        770        572
## ENSG00000000005          0          0          0
## ENSG00000000419        799        417        508
## ENSG00000000457        331        233        229
## ENSG00000000460         63         76         60
## ENSG00000000938          0          0          0

Step 2 — QC: library sizes

Same first check as the presentation — how many total reads landed in each sample.

lib_sizes <- colSums(counts)
lib_df <- data.frame(sample = names(lib_sizes), total_counts = lib_sizes, dex = meta$dex)

ggplot(lib_df, aes(x = sample, y = total_counts, fill = dex)) +
  geom_col() +
  scale_fill_manual(values = c("control" = NAVY, "treated" = ORANGE)) +
  labs(y = "Total counts (library size)", x = NULL,
       title = "Library size per sample (orange = Dex-treated, navy = control)") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 30, hjust = 1), legend.position = "none")

Step 3 — Filter low-count genes

keep <- rowSums(counts) >= 10
counts_filt <- counts[keep, ]
cat(sprintf("Genes before filtering: %d\n", nrow(counts)))
## Genes before filtering: 64102
cat(sprintf("Genes after filtering:  %d\n", nrow(counts_filt)))
## Genes after filtering:  22369

Step 4 — Normalize, estimate dispersion, fit the model

The design formula is ~ celltype + dex rather than just ~ dex — the celltype term tells DESeq2 to account for baseline differences between the 4 donor cell lines before estimating the treatment effect.

dds <- DESeqDataSetFromMatrix(
  countData = counts_filt,
  colData   = meta,
  design    = ~ celltype + dex
)
dds$dex <- relevel(factor(dds$dex), ref = "control")
dds <- DESeq(dds)
plotDispEsts(dds)

Step 5 — Test each gene: dexamethasone vs. control

res <- results(dds, contrast = c("dex", "treated", "control"))
res <- res[order(res$padj), ]

cat(sprintf("Genes tested: %d\n", nrow(res)))
## Genes tested: 22369
cat(sprintf("padj < 0.05 (genome-wide significant): %d\n", sum(res$padj < 0.05, na.rm = TRUE)))
## padj < 0.05 (genome-wide significant): 4000
head(as.data.frame(res), 10)
##                   baseMean log2FoldChange      lfcSE     stat        pvalue
## ENSG00000152583   997.4447       4.574967 0.18424142 24.83137 4.110667e-136
## ENSG00000165995   495.0957       3.291099 0.13305274 24.73530 4.463384e-135
## ENSG00000120129  3409.0384       2.947850 0.12187644 24.18720 3.033839e-129
## ENSG00000101347 12703.4128       3.767022 0.15599195 24.14882 7.682657e-129
## ENSG00000189221  2341.7807       3.353655 0.14218099 23.58722 5.212706e-123
## ENSG00000211445 12285.7001       3.730439 0.16639840 22.41872 2.585059e-111
## ENSG00000157214  3009.2729       1.976796 0.09042512 21.86113 6.091348e-106
## ENSG00000162614  5393.1145       2.035697 0.09468470 21.49974 1.565462e-102
## ENSG00000125148  3656.2674       2.211006 0.10609056 20.84074  1.849644e-96
## ENSG00000154734 30315.1132       2.345635 0.11635963 20.15849  2.266757e-90
##                          padj
## ENSG00000152583 7.412355e-132
## ENSG00000165995 4.024187e-131
## ENSG00000120129 1.823540e-125
## ENSG00000101347 3.463342e-125
## ENSG00000189221 1.879910e-119
## ENSG00000211445 7.768964e-108
## ENSG00000157214 1.569131e-102
## ENSG00000162614  3.528552e-99
## ENSG00000125148  3.705865e-93
## ENSG00000154734  4.087415e-87

A gene-symbol sanity check, the same idea as before: these are well-established glucocorticoid receptor target genes. If a new pipeline doesn’t recover them here, something is wrong with the pipeline, not the biology.

known_glucocorticoid_genes <- c(
  ENSG00000096060 = "FKBP5",
  ENSG00000101347 = "SAMHD1",
  ENSG00000163884 = "KLF15",
  ENSG00000157214 = "STEAP2"
)

known_ids <- intersect(names(known_glucocorticoid_genes), rownames(res))
known_df <- as.data.frame(res[known_ids, ])
rownames(known_df) <- known_glucocorticoid_genes[known_ids]
known_df
##         baseMean log2FoldChange      lfcSE     stat        pvalue          padj
## FKBP5   2564.384       4.046656 0.36977841 10.94346  7.141983e-28  9.005891e-26
## SAMHD1 12703.413       3.767022 0.15599195 24.14882 7.682657e-129 3.463342e-125
## KLF15    561.111       4.459169 0.23479342 18.99188  1.990619e-80  2.243428e-77
## STEAP2  3009.273       1.976796 0.09042512 21.86113 6.091348e-106 1.569131e-102

Step 6 — PCA

vsd <- vst(dds, blind = FALSE)
pca <- prcomp(t(assay(vsd)))
pct_var <- round(100 * (pca$sdev^2 / sum(pca$sdev^2)), 1)

pca_df <- data.frame(PC1 = pca$x[, 1], PC2 = pca$x[, 2], meta[rownames(pca$x), ])

ggplot(pca_df, aes(x = PC1, y = PC2, color = dex, label = celltype)) +
  geom_point(size = 4) +
  geom_text_repel(size = 3, show.legend = FALSE) +
  scale_color_manual(values = c("control" = NAVY, "treated" = ORANGE)) +
  labs(x = sprintf("PC1 (%.1f%%)", pct_var[1]), y = sprintf("PC2 (%.1f%%)", pct_var[2]),
       title = "PCA (VST counts) — cell line labeled") +
  theme_minimal()

PC1 separates treatment; PC2 largely separates cell line (you can see each cell line’s treated/control pair sitting at a similar position on PC2) — an illustration of exactly what the ~ celltype + dex design formula is modeling explicitly.

Step 7 — Volcano plot

res_plot <- as.data.frame(res)
res_plot$gene_id <- rownames(res_plot)
res_plot$neg_log10_padj <- -log10(pmax(res_plot$padj, 1e-300))
res_plot$group <- "Not significant"
res_plot$group[!is.na(res_plot$padj) & res_plot$padj < 0.05 & res_plot$log2FoldChange > 1]  <- "Up with Dex"
res_plot$group[!is.na(res_plot$padj) & res_plot$padj < 0.05 & res_plot$log2FoldChange < -1] <- "Down with Dex"

label_df <- res_plot[res_plot$gene_id %in% names(known_glucocorticoid_genes), ]
label_df$symbol <- known_glucocorticoid_genes[label_df$gene_id]

ggplot(res_plot, aes(x = log2FoldChange, y = neg_log10_padj, color = group)) +
  geom_point(alpha = 0.6, size = 1.2) +
  geom_text_repel(data = label_df, aes(label = symbol), color = "black",
                   fontface = "bold", size = 3.2, show.legend = FALSE) +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed", color = "grey40") +
  scale_color_manual(values = c("Not significant" = GREY, "Up with Dex" = ORANGE, "Down with Dex" = TEAL)) +
  labs(x = "log2 fold change (Dex vs. control)", y = "-log10(adjusted p-value)",
       title = "Volcano plot — airway smooth muscle cells (real DESeq2)", color = NULL) +
  theme_minimal()

Step 8 — Heatmap of the top differentially expressed genes

top_genes <- rownames(res)[1:20]
vst_mat <- assay(vsd)[top_genes, ]
z_mat <- t(scale(t(vst_mat)))  # row z-score

sample_order <- rownames(meta)[order(meta$dex, meta$celltype)]
z_mat <- z_mat[, sample_order]

annotation_col <- data.frame(dex = meta[sample_order, "dex"], row.names = sample_order)
ann_colors <- list(dex = c(control = NAVY, treated = ORANGE))

pheatmap(
  z_mat,
  cluster_rows = TRUE,
  cluster_cols = FALSE,
  annotation_col = annotation_col,
  annotation_colors = ann_colors,
  color = colorRampPalette(c(TEAL, "white", ORANGE))(100),
  main = "Top 20 genes by adjusted p-value"
)

Save results

write.csv(as.data.frame(res), "airway_DESeq2_DexVsControl_results.csv")
cat("Saved: airway_DESeq2_DexVsControl_results.csv\n")
## Saved: airway_DESeq2_DexVsControl_results.csv

Note: an even shorter path if you have the airway package

Since this is literally the dataset the airway Bioconductor package ships, you can skip the CSVs and manual data wrangling entirely by running the code below:

# BiocManager::install("airway")
library(airway)
data(airway)

dds <- DESeqDataSet(airway, design = ~ cell + dex)
dds$dex <- relevel(dds$dex, ref = "untrt")
dds <- DESeq(dds)
res <- results(dds, contrast = c("dex", "trt", "untrt"))
res <- res[order(res$padj), ]

vsd <- vst(dds, blind = FALSE)
plotPCA(vsd, intgroup = c("dex", "cell"))

airway’s own column names differ slightly (cell instead of celltype, trt/untrt instead of treated/control) but everything else — design formula, workflow, expected results — is identical to what this document just ran from the CSVs.