From Raw Counts to Differential Expression

Where we left off

In the flagship project you took a count matrix and a metadata table, ran two lines of DESeq2, and got a results table plus a volcano plot. That was the speed run. This lesson slows down and answers the questions that actually make you good at this: what is the design formula doing, what happens inside DESeq(), and how do I trust the numbers that come out?

1The design formula: the most important line you'll write

The design formula is how you tell DESeq2 what comparison you care about. In the flagship it was simply:

design = ~ condition

The tilde (~) means "model the counts as a function of...", and condition is the column in your metadata that splits the samples into groups (here, A and B). That one choice frames the whole analysis.

Real experiments are rarely that clean. Suppose your samples were processed in two different batches, and batch differences are leaking into your data. You can tell DESeq2 to account for that:

design = ~ batch + condition

This says: "explain away the differences due to batch first, then test the effect of condition on top of that." A simple rule worth memorizing: the variable you actually want to test goes last, and anything you want to control for (batch, sex, age) goes before it.

2What DESeq2 actually does (three moves)

When you call DESeq(dds), three things happen in sequence. You don't have to do them by hand, but knowing what they are makes the results make sense.

Move 1: Normalization (size factors). Some samples are sequenced more deeply than others, so every gene in them looks "higher" just by accident. DESeq2 computes a size factor per sample and divides it out, so counts become comparable across samples. This corrects for sequencing depth and for differences in library composition between samples (it does not correct for gene length, which is fine because differential expression compares the same gene across samples). Correcting for composition, not just total counts, is exactly why DESeq2 size factors differ from naive total-count scaling.

Move 2: Dispersion estimation. With only a few replicates per group, the per-gene variability is hard to estimate. DESeq2 borrows strength across all genes: it looks at how variable genes are at each expression level and shrinks each gene's noisy estimate toward that trend. This is why DESeq2 works well even with three replicates, where naive methods fall apart.

Move 3: Model and test. DESeq2 fits a negative binomial model (the right shape for count data) and, for each gene, runs a statistical test asking whether the difference between groups is bigger than the noise. Out comes a log2FoldChange, a pvalue, and an adjusted p-value.

3Run it, with two good habits

Here is the workflow again, with two additions that real analyses use: filtering out genes with almost no reads, and checking the coefficient names.

library(DESeq2)

base <- "https://raw.githubusercontent.com/owkin/PyDESeq2/main/datasets/synthetic/"
counts <- as.matrix(read.csv(paste0(base, "test_counts.csv"), row.names = 1))
metadata <- read.csv(paste0(base, "test_metadata.csv"), row.names = 1)
metadata$condition <- factor(metadata$condition)
metadata <- metadata[colnames(counts), , drop = FALSE]

dds <- DESeqDataSetFromMatrix(countData = counts,
                              colData   = metadata,
                              design    = ~ condition)

# Habit 1: drop genes with almost no reads (faster, less noise)
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep, ]

dds <- DESeq(dds)

# Habit 2: see the coefficients DESeq2 estimated (you'll need a name in step 5)
resultsNames(dds)

4Read the results table properly

Pull the comparison you want, then read what came back:

res <- results(dds, contrast = c("condition", "B", "A"))
summary(res)                       # a quick overview: how many up / down
head(res[order(res$padj), ], 10)   # the 10 most convincing genes

Every column earns its place:

• baseMean: the average normalized count for that gene across all samples. A sanity check on how expressed the gene is.
• log2FoldChange: the effect size, on a log2 scale. Plus 1 means doubled in B vs A, minus 1 means halved. Sign is direction, size is magnitude.
• lfcSE: the standard error of that fold change, so how uncertain it is.
• stat: the test statistic (the fold change divided by its error).
• pvalue: the raw probability of seeing a change this big by chance.
• padj: the p-value adjusted for testing thousands of genes at once. This is the column you trust.

To pull a clean list of hits, filter on both significance and effect size:

# significant AND at least a two-fold change
sig <- subset(res, padj < 0.05 & abs(log2FoldChange) > 1)
nrow(sig)

One thing that surprises beginners: some genes have NA in the padj column. That is not a bug. DESeq2 deliberately sets aside genes that are too low-count to ever be detectable, which actually makes your real hits more significant. Treat NA padj as "not testable", not "failed".

5Cleaner effect sizes: LFC shrinkage

A low-count gene can show a huge fold change purely from noise (two reads versus eight reads is "4x", but it's meaningless). Shrinkage fixes this: it pulls unreliable fold changes toward zero while leaving well-measured genes almost untouched. The result is rankings you can trust and much cleaner volcano and MA plots.

# apeglm is the recommended shrinkage method; install it once if needed
# BiocManager::install("apeglm")
res_shrunk <- lfcShrink(dds, coef = "condition_B_vs_A", type = "apeglm")
head(res_shrunk[order(res_shrunk$padj), ], 10)

The coef name comes straight from resultsNames(dds) in step 3. Use the shrunken results whenever you rank genes or draw plots; use the unshrunken res when you just need the raw test outcome.

6Save your work

Write the table to a file so you can open it in Excel, share it, or feed it into the next lesson (volcano plots and enrichment both start from this table):

write.csv(as.data.frame(res), "deseq2_results.csv")

Check your understanding

Next in Track 3

Reading volcano & MA plots in R →