Differential expression analysis can feel like a black box if you just follow a tutorial blindly. Today, we're going to walk through a real scenario: imagine we have RNA-seq data from a set of lung cancer cells. Half were treated with a new kinase inhibitor, and half were left as a control. We want to know which genes were significantly "turned down" or "turned up" because of that drug.

Getting our counts and sample info ready

For DESeq2 to work, you need two things: a matrix of raw integer counts (not normalized values!) and a data frame describing your samples. I've put together a small synthetic dataset here to keep things focused.

library(DESeq2)

# Raw counts: rows are genes, columns are samples
countData <- matrix(rnbinom(1000, mu=100, size=1), ncol=6)
colnames(countData) <- c("Ctrl1", "Ctrl2", "Ctrl3", "Trt1", "Trt2", "Trt3")
rownames(countData) <- paste0("Gene", 1:167) # A few genes for the example

# Sample metadata
colData <- data.frame(
  condition = factor(c("Control", "Control", "Control", "Treated", "Treated", "Treated")),
  row.names = colnames(countData)
)

Notice I wrapped condition in factor(). DESeq2 relies heavily on factors to understand the experimental design. If you leave them as strings, you're asking for a headache later.

Building the DESeqDataSet object

Now we wrap these into a DESeqDataSet. This is the core object where all the magic happens. We tell it the counts, the metadata, and the design formula—which basically tells R, "I want to see how the expression changes based on the condition."

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

The "Reference Level" trap

Here is where I almost messed up this analysis. I ran DESeq(dds) immediately, but I forgot that R assigns factor levels alphabetically by default. In our case, "Control" comes before "Treated," so it worked out. But if my groups were "Placebo" and "Control," R would have treated "Control" as the baseline and "Placebo" as the treatment, flipping the sign of all my Fold Change results.

I always explicitly set the reference level to avoid this ambiguity. It's a small step that saves hours of debugging later.

# Force "Control" to be the baseline
dds$condition <- relevel(dds$condition, ref = "Control")

# Now we run the actual differential expression pipeline
dds <- DESeq(dds)

Extracting the meaningful genes

The DESeq() function does a lot under the hood: it estimates size factors for normalization, estimates dispersion, and fits a generalized linear model. To actually see the results, we use results(). I usually filter for an adjusted p-value (padj) of less than 0.05 to account for the thousands of tests we're running simultaneously.

res <- results(dds)

# Let's look at the most significant genes
resOrdered <- res[order(res$padj), ]
head(resOrdered)

If you see a positive log2FoldChange, the gene is upregulated in the Treated group compared to the Control. A negative value means it was suppressed. Simple, but incredibly powerful.

📋 Practical Task

Identify Drug-Responsive Genes from a New Dataset

You have been handed a new count matrix lung_counts and a metadata table lung_meta. The metadata contains a column called treatment with two levels: "DMSO" (the vehicle control) and "Drug_A".

Write a script to:

• Create a DESeqDataSet object using treatment as the design.
• Explicitly set "DMSO" as the reference level for the treatment factor.
• Run the DESeq pipeline and extract the results.
• Create a new data frame called significant_genes that contains only the genes with an adjusted p-value (padj) less than 0.01.