Airway RNA-seq: Differential expression showcase

Load data and metadata

library(airway)
library(DESeq2)
library(EnhancedVolcano)
library(pheatmap)
library(ggplot2)
library(RColorBrewer)

theme_set(theme_minimal(base_size = 13))

data("airway")

meta <- as.data.frame(colData(airway))
meta$dex <- droplevels(meta$dex)
meta$cell <- droplevels(meta$cell)

knitr::kable(meta, caption = "Sample metadata (8 airway RNA-seq libraries)")

Sample metadata (8 airway RNA-seq libraries)
	SampleName	cell	dex	albut	Run	avgLength	Experiment	Sample	BioSample
SRR1039508	GSM1275862	N61311	untrt	untrt	SRR1039508	126	SRX384345	SRS508568	SAMN02422669
SRR1039509	GSM1275863	N61311	trt	untrt	SRR1039509	126	SRX384346	SRS508567	SAMN02422675
SRR1039512	GSM1275866	N052611	untrt	untrt	SRR1039512	126	SRX384349	SRS508571	SAMN02422678
SRR1039513	GSM1275867	N052611	trt	untrt	SRR1039513	87	SRX384350	SRS508572	SAMN02422670
SRR1039516	GSM1275870	N080611	untrt	untrt	SRR1039516	120	SRX384353	SRS508575	SAMN02422682
SRR1039517	GSM1275871	N080611	trt	untrt	SRR1039517	126	SRX384354	SRS508576	SAMN02422673
SRR1039520	GSM1275874	N061011	untrt	untrt	SRR1039520	101	SRX384357	SRS508579	SAMN02422683
SRR1039521	GSM1275875	N061011	trt	untrt	SRR1039521	98	SRX384358	SRS508580	SAMN02422677

Raw counts snapshot and library sizes

counts <- assay(airway)
lib_sizes <- colSums(counts)

ggplot(
  data.frame(
    sample = rownames(meta),
    lib_millions = lib_sizes / 1e6,
    dex = meta$dex,
    cell = meta$cell
  ),
  aes(x = reorder(sample, lib_millions), y = lib_millions, fill = dex)
) +
  geom_col(color = "white", width = 0.8) +
  coord_flip() +
  labs(
    x = "Sample",
    y = "Library size (millions of reads)",
    fill = "Dexamethasone",
    title = "Library size per sample"
  ) +
  geom_hline(
    yintercept = median(lib_sizes / 1e6),
    linetype = "dashed",
    color = "gray30"
  ) +
  theme(plot.title = element_text(face = "bold"))

DESeq2 object and normalization

dds <- DESeqDataSetFromMatrix(
  countData = counts,
  colData = meta,
  design = ~ cell + dex
)

dds <- dds[rowSums(counts(dds)) > 1, ]

dds <- DESeq(dds)
vsd <- vst(dds, blind = FALSE)

Sample-level QC

sampleDists <- dist(t(assay(vsd)))
sampleDistMatrix <- as.matrix(sampleDists)
rownames(sampleDistMatrix) <- colnames(vsd)
colnames(sampleDistMatrix) <- colnames(vsd)

ann <- meta[, c("dex", "cell")]

pheatmap(
  sampleDistMatrix,
  clustering_distance_rows = sampleDists,
  clustering_distance_cols = sampleDists,
  annotation_col = ann,
  main = "Sample-to-sample distances"
)

pca_data <- plotPCA(vsd, intgroup = c("dex", "cell"), returnData = TRUE)
percentVar <- round(100 * attr(pca_data, "percentVar"))

ggplot(pca_data, aes(PC1, PC2, color = dex, shape = cell)) +
  geom_point(size = 4, alpha = 0.9) +
  labs(
    x = paste0("PC1: ", percentVar[1], "%"),
    y = paste0("PC2: ", percentVar[2], "%"),
    title = "PCA on vst-normalized counts"
  ) +
  theme(plot.title = element_text(face = "bold"))

Differential expression results

coef_name <- tail(resultsNames(dds), 1)

res_raw <- results(dds, name = coef_name)
res <- lfcShrink(dds, coef = coef_name, type = "normal")

res_df <- as.data.frame(res[order(res$padj), ])

knitr::kable(
  head(res_df, 10),
  digits = 3,
  caption = "Top 10 differentially expressed genes (shrunken LFC)"
)

Top 10 differentially expressed genes (shrunken LFC)
	baseMean	log2FoldChange	lfcSE	stat
ENSG00000152583	997.440	-4.079	0.162	-24.856
ENSG00000165995	495.093	-3.089	0.124	-24.713
ENSG00000120129	3409.029	-2.796	0.115	-24.274
ENSG00000101347	12703.387	-3.475	0.143	-24.235
ENSG00000189221	2341.767	-3.116	0.132	-23.653
ENSG00000211445	12285.615	-3.389	0.151	-22.495
ENSG00000157214	3009.263	-1.922	0.087	-21.968
ENSG00000162614	5393.102	-1.973	0.091	-21.614
ENSG00000125148	3656.253	-2.125	0.102	-20.929
ENSG00000154734	30315.135	-2.233	0.110	-20.254

summary(res_raw)

## 
## out of 29391 with nonzero total read count
## adjusted p-value < 0.1
## LFC > 0 (up)       : 2218, 7.5%
## LFC < 0 (down)     : 2607, 8.9%
## outliers [1]       : 0, 0%
## low counts [2]     : 11397, 39%
## (mean count < 5)
## [1] see 'cooksCutoff' argument of ?results
## [2] see 'independentFiltering' argument of ?results

cat("Genes with padj < 0.05:", sum(res_raw$padj < 0.05, na.rm = TRUE), "\n")

## Genes with padj < 0.05: 3994

Global views of differential signal

plotMA(res_raw, ylim = c(-5, 5), main = "MA plot (unshrunken results)")

EnhancedVolcano(
  res,
  lab = rownames(res),
  x = "log2FoldChange",
  y = "padj",
  xlab = "log2 fold change",
  ylab = "-log10(padj)",
  title = "Volcano: dex treatment effect",
  subtitle = coef_name,
  pCutoff = 0.05,
  FCcutoff = 1.0,
  pointSize = 2.5,
  labSize = 3.5,
  colAlpha = 0.7,
  legendPosition = "right",
  legendLabSize = 10
)

Heatmap of strongest changes

top_genes <- rownames(res_df)[1:30]
top_genes <- top_genes[!is.na(top_genes)]

mat <- assay(vsd)[top_genes, ]
mat <- mat - rowMeans(mat)

pheatmap(
  mat,
  annotation_col = ann,
  cluster_rows = TRUE,
  cluster_cols = TRUE,
  show_rownames = TRUE,
  color = colorRampPalette(rev(brewer.pal(9, "RdBu")))(50),
  main = "Top 30 DE genes (z-scored per gene)"
)

Spotlight on the top gene

top_gene <- top_genes[1]

gene_df <- data.frame(
  expression = assay(vsd)[top_gene, ],
  dex = meta$dex,
  cell = meta$cell,
  sample = rownames(meta)
)

ggplot(gene_df, aes(dex, expression, fill = dex)) +
  geom_boxplot(alpha = 0.7, outlier.color = NA) +
  geom_jitter(aes(color = cell), width = 0.12, size = 3, alpha = 0.9) +
  labs(
    title = paste("Expression of", top_gene),
    x = "Condition",
    y = "vst expression"
  ) +
  scale_fill_brewer(palette = "Set2") +
  scale_color_brewer(palette = "Dark2") +
  theme(plot.title = element_text(face = "bold"))

Session info

sessionInfo()

## R version 4.3.2 (2023-10-31 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 11 x64 (build 26100)
## 
## Matrix products: default
## 
## 
## locale:
## [1] LC_COLLATE=English_United States.utf8 
## [2] LC_CTYPE=English_United States.utf8   
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.utf8    
## 
## time zone: Asia/Singapore
## tzcode source: internal
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] RColorBrewer_1.1-3          pheatmap_1.0.13            
##  [3] EnhancedVolcano_1.20.0      ggrepel_0.9.6              
##  [5] ggplot2_4.0.1               DESeq2_1.42.1              
##  [7] airway_1.22.0               SummarizedExperiment_1.32.0
##  [9] Biobase_2.62.0              GenomicRanges_1.54.1       
## [11] GenomeInfoDb_1.38.8         IRanges_2.36.0             
## [13] S4Vectors_0.40.2            BiocGenerics_0.48.1        
## [15] MatrixGenerics_1.14.0       matrixStats_1.5.0          
## 
## loaded via a namespace (and not attached):
##  [1] gtable_0.3.6            xfun_0.52               bslib_0.9.0            
##  [4] lattice_0.21-9          vctrs_0.6.5             tools_4.3.2            
##  [7] bitops_1.0-9            generics_0.1.3          parallel_4.3.2         
## [10] tibble_3.2.1            fansi_1.0.6             pkgconfig_2.0.3        
## [13] Matrix_1.6-1.1          S7_0.2.0                lifecycle_1.0.4        
## [16] GenomeInfoDbData_1.2.11 compiler_4.3.2          farver_2.1.1           
## [19] codetools_0.2-19        htmltools_0.5.8.1       sass_0.4.9             
## [22] RCurl_1.98-1.17         yaml_2.3.10             pillar_1.9.0           
## [25] crayon_1.5.2            jquerylib_0.1.4         BiocParallel_1.36.0    
## [28] DelayedArray_0.28.0     cachem_1.1.0            abind_1.4-8            
## [31] tidyselect_1.2.1        locfit_1.5-9.12         digest_0.6.37          
## [34] dplyr_1.1.4             labeling_0.4.3          fastmap_1.2.0          
## [37] grid_4.3.2              colorspace_2.1-0        cli_3.6.2              
## [40] SparseArray_1.2.4       magrittr_2.0.3          S4Arrays_1.2.1         
## [43] utf8_1.2.4              withr_3.0.2             scales_1.4.0           
## [46] rmarkdown_2.30          XVector_0.42.0          evaluate_1.0.3         
## [49] knitr_1.50              rlang_1.1.3             Rcpp_1.0.12            
## [52] glue_1.7.0              rstudioapi_0.17.1       jsonlite_2.0.0         
## [55] R6_2.5.1                zlibbioc_1.48.2

Airway RNA-seq: Differential expression showcase

Mikhail Bovt

2025-11-28

Load data and metadata

Raw counts snapshot and library sizes

DESeq2 object and normalization

Sample-level QC

Differential expression results

Global views of differential signal

Heatmap of strongest changes

Spotlight on the top gene

Session info