Differential Expression exercise

differential_expression.md

Introduction

In this exercise, you will analyze RNA-seq data from a p53 activation experiment. Cells were either mock-treated (control) or exposed to ionizing radiation (IR), which activates the p53 tumor suppressor pathway. There are 4 biological replicates per condition (8 samples total).

Data was obtained from https://github.com/pacthoen/BMW2_RNA_clust_vis. To get a copy, drop into a terminal, and download the following files:

wget https://raw.githubusercontent.com/pacthoen/BMW2_RNA_clust_vis/refs/heads/main/data/p53_mock_1.csv
wget https://raw.githubusercontent.com/pacthoen/BMW2_RNA_clust_vis/refs/heads/main/data/p53_mock_2.csv
wget https://raw.githubusercontent.com/pacthoen/BMW2_RNA_clust_vis/refs/heads/main/data/p53_mock_3.csv
wget https://raw.githubusercontent.com/pacthoen/BMW2_RNA_clust_vis/refs/heads/main/data/p53_mock_4.csv
wget https://raw.githubusercontent.com/pacthoen/BMW2_RNA_clust_vis/refs/heads/main/data/p53_IR_1.csv
wget https://raw.githubusercontent.com/pacthoen/BMW2_RNA_clust_vis/refs/heads/main/data/p53_IR_2.csv
wget https://raw.githubusercontent.com/pacthoen/BMW2_RNA_clust_vis/refs/heads/main/data/p53_IR_3.csv
wget https://raw.githubusercontent.com/pacthoen/BMW2_RNA_clust_vis/refs/heads/main/data/p53_IR_4.csv

In this exercise, we'll take this data and:

1. Read in and merge the per-sample count files into a single count matrix
2. Perform quality control by visualizing sample clustering with an MDS plot
3. Identify differentially expressed genes between IR and mock using edgeR
4. Visualize the differentially expressed genes with a heatmap using pheatmap and a volcano plot

The data consist of raw read counts — the number of sequencing reads that mapped to each gene in each sample. This is the starting point for most RNA-seq differential expression workflows.

---

1. Setup

Load the packages you will need. If any are missing, install them first with BiocManager::install("edgeR") and install.packages("pheatmap").

library(edgeR)
library(pheatmap)

Change your working directory to the same directory in which that you saved the gene expression csv files. For example:

setwd("/Users/cmiller/projects/p53data/")

---

2. Reading in the Data

The data are stored as individual CSV files — one per sample. Each file has two columns: the Ensembl gene ID and the read count for that sample. The first step is to read all eight files and merge them into a single count matrix where rows are genes and columns are samples.

files <- Sys.glob("p53*.csv")
files

You should see 8 files listed. Note that they are sorted alphabetically, so the IR samples come before the mock samples. This ordering will matter when we define our experimental groups below.

Now read each file and merge them column-by-column using the Ensembl gene ID as the common key:

#lapply applies the read.csv function to each element of files
counts_list <- lapply(files, function(f) {
    read.csv(f, row.names = 1)
})
#cbind merges the columns
counts <- do.call(cbind, counts_list)

dim(counts)
head(counts)

You should have a matrix with 6,882 rows (genes) and 8 columns (samples). The column names come directly from the headers in each CSV file.

---

3. Building the DGEList

edgeR works with a specialized data structure called a DGEList, which bundles the count matrix together with sample metadata and (later) normalization factors.

First, define the experimental groups. We extract group labels directly from the column names so the grouping is always consistent with the data:

group <- factor(ifelse(grepl("mock", colnames(counts)), "mock", "IR"))
group

Now create the DGEList:

dge <- DGEList(counts = counts, group = group)

Filtering low-count genes

Many genes will have very few counts across all samples — these carry little statistical information and can hurt your ability to detect real differences. The filterByExpr() function uses the library sizes and group structure to automatically determine a sensible count threshold:

keep <- filterByExpr(dge)
dge <- dge[keep, , keep.lib.sizes = FALSE]

dim(dge)

This typically removes a large fraction of genes. The retained genes are those with enough counts to be meaningfully analyzed.

Normalization

Raw counts are affected by differences in sequencing depth (library size) between samples — a gene that appears to have twice as many counts in one sample may simply reflect that sample having twice as many total reads. TMM normalization (Trimmed Mean of M-values) estimates scaling factors to account for this, making counts comparable across samples:

dge <- calcNormFactors(dge)
dge$samples

The norm.factors column shows the TMM scaling factor for each sample. Values close to 1 indicate a sample's composition is close to the overall average.

---

4. Quality Control: MDS Plot

Before running any statistical tests, it is good practice to visualize how the samples relate to one another. An MDS plot (Multidimensional Scaling) is a 2D projection that places samples close together if their expression profiles are similar. In an ideal experiment, you'd want to see:

• Within-group samples clustering together (low biological noise)
• Between-group separation (a signal worth detecting)

sample_colors <- c(IR = "tomato", mock = "steelblue")

plotMDS(
    dge,
    col = sample_colors[group],
    pch = 16,
    cex = 1.5,
    main = "MDS Plot: IR vs Mock"
)
legend(
    "topright",
    legend = levels(group),
    col = sample_colors[levels(group)],
    pch = 16
)

Examine the plot: do IR and mock samples separate cleanly along one of the axes? Are there any obvious outliers within a group? If samples don't cluster by treatment, it may indicate a batch effect or sample mix-up that should be investigated before proceeding.

---

5. Differential Expression with edgeR

Estimating dispersion

RNA-seq count data are often overdispersed — the variance is larger than you'd expect from a simple Poisson model. edgeR models this using the negative binomial distribution, which has an additional parameter called dispersion. estimateDisp() estimates this from the data, borrowing information across all genes to improve the estimate for any single gene:

dge <- estimateDisp(dge)

You can visualize the estimated dispersions across genes:

plotBCV(dge, main = "Biological Coefficient of Variation")

The BCV plot shows the estimated variability for each gene as a function of its expression level. The red line is the trended dispersion estimate that edgeR uses in testing. Lower BCV values generally mean more statistical power to detect differential expression.

Running the test

The exact test is a classical approach for testing differential expression between two groups in RNA-seq data. The pair argument specifies the comparison direction: c("mock", "IR") means the fold change is calculated as IR relative to mock, so a positive logFC indicates higher expression in IR:

et <- exactTest(dge, pair = c("mock", "IR"))

Extracting results

topTags() returns the results as a ranked table. Setting n = Inf retrieves all genes (not just the top 10):

results <- topTags(et, n = Inf)$table
head(results)

The columns are:

Column	Meaning
logFC	Log2 fold-change (IR vs mock); positive = higher in IR
logCPM	Average log2 counts-per-million across all samples
PValue	Raw p-value from the exact test
FDR	P-value adjusted for multiple testing (Benjamini-Hochberg)

Applying significance thresholds

We define differentially expressed genes as those with FDR < 0.05 and |logFC| > 1 (at least 2-fold change):

sig_genes <- results[results$FDR < 0.05 & abs(results$logFC) > 1, ]

cat("Total DE genes:", nrow(sig_genes), "\n")
cat("Upregulated in IR:", sum(sig_genes$logFC > 0), "\n")
cat("Downregulated in IR:", sum(sig_genes$logFC < 0), "\n")

Look at the top differentially expressed genes:

head(sig_genes, 20)

---

6. Heatmap of Differentially Expressed Genes

A heatmap gives a gene-by-sample view of expression for your DE genes. We'll use log-CPM values (log2 counts-per-million), which are on a continuous scale suitable for visualization. We then scale each row (z-score per gene) so that the color represents relative expression — how each sample compares to the mean for that gene — rather than absolute expression level, which would be dominated by highly expressed genes.

## log-CPM values from the normalized DGEList
logcpm <- cpm(dge, log = TRUE)

## subset to DE genes only
heatmap_mat <- logcpm[rownames(sig_genes), ]

## annotation data frame for the column (sample) color bar
anno_col <- data.frame(
    Treatment = group,
    row.names = colnames(heatmap_mat)
)

anno_colors <- list(
    Treatment = c(mock = "steelblue", IR = "tomato")
)

pheatmap(
    heatmap_mat,
    scale = "row",
    annotation_col = anno_col,
    annotation_colors = anno_colors,
    show_rownames = FALSE,
    cluster_rows = TRUE,
    cluster_cols = TRUE,
    color = colorRampPalette(c("navy", "white", "firebrick3"))(100),
    main = "DE Genes: IR vs Mock (FDR < 0.05, |logFC| > 1)"
)

Interpreting the heatmap:

• Each row is a gene; each column is a sample
• Color represents relative expression (z-scored): red = higher than average, blue = lower than average
• Row clustering groups genes with similar expression patterns
• Column clustering can reveal whether samples separate by treatment
• Look for blocks of genes that are consistently up or down in IR relative to mock — these represent co-regulated gene sets that may share biological function

---

7. Volcano Plot

A volcano plot displays all tested genes at once, with fold-change on the x-axis and statistical significance on the y-axis. Unlike the heatmap — which only shows genes that passed your thresholds — the volcano plot lets you see the full landscape of results and appreciate where your significance cutoffs fall relative to the overall distribution.

We color-code points by their DE status so significant genes stand out clearly:

## classify each gene for coloring
de_status <- ifelse(
    results$FDR < 0.05 & results$logFC > 1,  "Up in IR",
    ifelse(
        results$FDR < 0.05 & results$logFC < -1, "Down in IR",
        "Not significant"
    )
)

point_colors <- c(
    "Up in IR"        = "tomato",
    "Down in IR"      = "steelblue",
    "Not significant" = "grey70"
)

plot(
    results$logFC,
    -log10(results$FDR),
    col = point_colors[de_status],
    pch = 16,
    cex = 0.5,
    xlab = "Log2 Fold Change (IR vs Mock)",
    ylab = "-log10(FDR)",
    main = "Volcano Plot: IR vs Mock"
)

## threshold lines to make the cutoffs visible
abline(h = -log10(0.05), lty = 2, col = "grey40")
abline(v = c(-1, 1),     lty = 2, col = "grey40")

legend(
    "topright",
    legend = names(point_colors),
    col = point_colors,
    pch = 16,
    pt.cex = 1
)

Interpreting the volcano plot:

• Points in the upper corners are the most statistically significant and strongly changed genes — the most compelling DE candidates
• The dashed lines mark your thresholds: FDR = 0.05 (horizontal) and |logFC| = 1 (vertical)
• Genes near the x-axis have large fold changes but poor statistical support — often caused by high variability across replicates
• Genes near the y-axis are highly significant but have small fold changes — statistically real but possibly not biologically meaningful

---

Additional steps beyond this exercise

We could do lots of things as next steps. Some examples:

• use the Ensembl gene IDs to look up gene names and functions (e.g., via org.Mm.eg.db or the Ensembl website)
• Run gene set enrichment or overrepresentation analysis on the DE gene list to identify enriched biological pathways
• add labels for specific genes of interest to our volcano plots