R-Package Presentation: DESeq2

Laura Henze

What is DESeq2?

  • DESeq2 is an R package for analyzing RNA-seq count data
  • It identifies genes that are expressed differently between conditions (ex. control vs treated)
  • Widely used in biology, genomics, and medicine

Why does this matter?

  • Helps answer:
    • Which genes are active in diseased vs healthy cells?
    • How do treatments affect gene expression?

Background: RNA-seq

  • RNA sequencing measures gene expression
  • Output = counts of how many RNA fragments come from each gene
  • More reads = higher gene activity

Challenges with RNA-seq Data

  • Data are counts (not continuous)
  • Different samples have different sequencing depth
  • High variability in genes with low expresssion

What DESeq2 Does

  • Normalizes count data
  • Models counts using a negative binomial distribution
  • Tests for differential expression

Load Packages

# Install Packages
install.packages("BiocManager")
BiocManager::install("DESeq2")

library(DESeq2)
# Add dplyr for data manipulation
library(dplyr)

Create Dataset

set.seed(123)

# Base counts
count_data <- matrix(
  rnbinom(1000, mu = 50, size = 1),
  nrow = 100,
  ncol = 10
)

rownames(count_data) <- paste0("Gene", 1:100) # rows = genes
colnames(count_data) <- paste0("Sample", 1:10) # columns = samples

# Metadata
col_data <- data.frame(
  condition = rep(c("Control", "Treated"), each = 5) # 2 conditions
)
rownames(col_data) <- colnames(count_data)

# First 20 genes MUCH higher in treated group
count_data[1:20, 6:10] <- count_data[1:20, 6:10] * 5

Prepare Data

# convert to DESeq2 dataset object
dds <- DESeqDataSetFromMatrix(
  countData = count_data, 
  colData = col_data,
  design = ~ condition # specify experimental design
)

Run DESeq2

dds <- DESeq(dds) # performs normalization and runs statistical tests

Get Results

log2FoldChange → size of change

res <- results(dds)

head(res)
log2 fold change (MLE): condition Treated vs Control 
Wald test p-value: condition Treated vs Control 
DataFrame with 6 rows and 6 columns
       baseMean log2FoldChange     lfcSE      stat      pvalue      padj
      <numeric>      <numeric> <numeric> <numeric>   <numeric> <numeric>
Gene1   57.3574        2.32134  1.054971   2.20038 0.027779895 0.1432858
Gene2  127.4924        2.71200  0.892098   3.04003 0.002365571 0.0352274
Gene3   99.9888        3.07118  0.902736   3.40208 0.000668748 0.0336437
Gene4  109.2556        1.95675  0.671892   2.91231 0.003587714 0.0352274
Gene5  147.4220        1.65342  0.732936   2.25589 0.024077461 0.1387995
Gene6  218.1978        2.79957  0.993186   2.81878 0.004820723 0.0393692

Significant Genes

sig_genes <- res %>%
  as.data.frame() %>%
  filter(!is.na(padj) & padj < 0.05) # filter for genes under adjusted p-value threshold

head(sig_genes)
        baseMean log2FoldChange     lfcSE     stat       pvalue       padj
Gene2  127.49239       2.712000 0.8920975 3.040027 0.0023655712 0.03522736
Gene3   99.98881       3.071181 0.9027360 3.402081 0.0006687476 0.03364372
Gene4  109.25562       1.956754 0.6718917 2.912306 0.0035877140 0.03522736
Gene6  218.19777       2.799569 0.9931863 2.818775 0.0048207229 0.03936924
Gene9  124.07012       2.257087 0.7462938 3.024395 0.0024913076 0.03522736
Gene11 208.48812       2.863022 0.9708600 2.948954 0.0031885121 0.03522736

MA Plot

  • each point = gene
  • x = average expression
  • y = how much the gene changes between conditions
  • far from zero, genes change the most
  • highlighted blue = statistically significant

MA Plot

plotMA(res, main = "MA Plot")

Top Genes

top_genes <- sig_genes |>
  arrange(padj) |> # top genes based on adjusted p-value
  head(10)
top_genes # genes with strong evidence of differential expression
        baseMean log2FoldChange     lfcSE      stat       pvalue       padj
Gene3   99.98881       3.071181 0.9027360  3.402081 0.0006687476 0.03364372
Gene12  94.64533       2.639304 0.8041200  3.282226 0.0010299099 0.03364372
Gene19 156.71391       2.468916 0.7514423  3.285570 0.0010177624 0.03364372
Gene2  127.49239       2.712000 0.8920975  3.040027 0.0023655712 0.03522736
Gene4  109.25562       1.956754 0.6718917  2.912306 0.0035877140 0.03522736
Gene9  124.07012       2.257087 0.7462938  3.024395 0.0024913076 0.03522736
Gene11 208.48812       2.863022 0.9708600  2.948954 0.0031885121 0.03522736
Gene14 181.02964       2.816623 0.9351417  3.011975 0.0025955427 0.03522736
Gene16 172.24764       2.218786 0.7699304  2.881801 0.0039540916 0.03522736
Gene75  48.09083      -2.276942 0.7752239 -2.937141 0.0033125366 0.03522736

Why DESeq2 is Powerful

  • Works with small sample sizes
  • Controls false positives
  • Widely used in research

Limitations

  • Requires raw count data
  • Assumes a specific statistical model
  • Can be slow for large datasets

Summary

  • DESeq2 analyzes RNA-seq data
  • Finds differentially expressed genes
  • Uses strong statistical methods
  • Essential tool in modern biology

References

  • Love et al. (2014), Genome Biology
  • Bioconductor DESeq2 documentation