Exercise 8 — Your First Differential Expression Analysis with DESeq2

Welcome to your first real RNA-seq differential expression analysis.

In previous exercises you learned:

  • expression matrices
  • metadata
  • fold change
  • visualization
  • heatmaps

In this exercise you will use DESeq2 to identify genes that respond to treatment.

This workflow is similar to what is performed in many published RNA-seq studies.

Biological Background

Dexamethasone is a glucocorticoid drug that alters gene expression.

The dataset contains RNA-seq profiles from:

  • untreated airway smooth muscle cells
  • dexamethasone-treated airway smooth muscle cells Your goal is to identify genes that respond to treatment.

Part 1 — Load Packages

Task 1.1

Install if necessary:

BiocManager::install("DESeq2")
BiocManager::install("airway")

Task 1.2

Load packages

library(DESeq2)
library(airway)
library(ggplot2)

Part 2 — Load Dataset

Task 2.1

Load the airway dataset

data("airway")

Store as:

airway

Task 2.2

Inspect the object

airway

Questions:

  1. How many genes are present?
  2. How many samples are present?
  3. What type of object is this?

Task 2.3

Inspect metadata

colData(airway)

Questions:

  1. What variables are available?
  2. Which variable describes treatment?

Part 3 — Explore Count Matrix

Task 3.1

Extract counts

counts_matrix <- assay(airway)

Task 3.2

Inspect dimensions

dim(counts_matrix)

Questions:

  1. Number of genes?
  2. Number of samples?

Task 3.3

View first few genes

head(counts_matrix)

Part 4 — Create DESeq Dataset

Task 4.1

Create DESeqDataSet

dds <- DESeqDataSet(
    airway,
    design = ~ dex
)

Task 4.2

Inspect object

dds

Question:

What does:

design = ~ dex

mean biologically?

Part 5 — Pre-filtering

Low count genes often add noise.

Task 5.1

Keep genes with:

rowSums(counts(dds)) >= 10

Store filtered object.

Task 5.2

How many genes were removed?

Part 6 — Run DESeq2

Task 6.1

Run analysis

dds <- DESeq(dds)

This may take a few minutes.

Task 6.2

Question

What steps do you think DESeq2 performs internally?

(You do NOT need exact details.)

Part 7 — Extract Results

Task 7.1

Generate results

res <- results(dds)

Task 7.2

Inspect results

head(res)
summary(res)

Questions:

  1. What is log2FoldChange?
  2. What is pvalue?
  3. What is padj?

Task 7.3

Convert to data frame

res_df <- as.data.frame(res)

Part 8 — Significant Genes

Task 8.1

Find genes with:

padj < 0.05

Store as:

sig_genes

Task 8.2

How many significant genes were identified?

Task 8.3

Find genes with:

padj < 0.05

AND

abs(log2FoldChange) > 1

Task 8.4

Why might we require both criteria?

Part 9 — Top Genes

Task 9.1

Sort genes by adjusted p-value

Task 9.2

Display top 20 genes

Task 9.3

Choose 3 genes.

For each gene answer:

  • Upregulated or downregulated?
  • Strong effect or weak effect?
  • Significant or not?

Part 10 — MA Plot

Task 10.1

Create MA plot

plotMA(res)

Task 10.2

Interpret

Questions:

  1. What does the x-axis represent?
  2. What does the y-axis represent?
  3. Why are low-count genes more variable?

Part 11 — Volcano Plot

Task 11.1

Create plotting dataframe

volcano_df <- as.data.frame(res)

Task 11.2

Create volcano plot

ggplot(volcano_df,
       aes(log2FoldChange,
           -log10(padj))) +
    geom_point()

Task 11.3

Add significance coloring

Significant:

padj < 0.05

Task 11.4

Interpret

Questions:

  1. Which quadrant contains strongly upregulated genes?
  2. Which quadrant contains strongly downregulated genes?
  3. What makes a good candidate biomarker?

Part 12 — Biological Interpretation

Answer:

# What does dexamethasone appear to do?
# Are there many or few responding genes?
# What is the strongest responding gene?
# Why are adjusted p-values necessary?
# Why is fold change alone insufficient?

Part 13 — Mini Analyst Report

Write 5–10 sentences summarizing:

  • dataset
  • comparison performed
  • number of significant genes
  • top findings
  • biological interpretation

Imagine you are presenting to your PI.

Bonus Challenge

B1

Identify:

Top 10 upregulated genes

B2

Identify:

Top 10 downregulated genes

B3

Create a heatmap of:

Top 20 significant genes

Hints:

library(pheatmap)

pheatmap(...)

Key Concepts Learned

  • DESeq2 workflow
  • differential expression
  • significance testing
  • adjusted p-values
  • MA plots
  • volcano plots
  • biological interpretation
  • result prioritization

This workflow forms the foundation of most bulk RNA-seq studies.

⬅️ Previous Exercise 🏠 Home ➡️ Next Exercise