I've got a set of raw read counts from a project I'm working on—lung epithelial cells, some treated with a kinase inhibitor and some left as controls. I have the counts.csv (genes in rows, samples in columns) and a metadata.csv that tells me which sample belongs to which group. At first glance, it's just a big table of integers. You might be tempted to jump straight into some basic statistics, and honestly, that's exactly where I started this morning.

Wait, why can't I just use a t-test?

My first instinct was to just loop through the genes and run a t.test() between the treatment and control groups. It seems logical, right? Let's see what happens when we try that on a few rows:

# Quick and dirty t-test on the first gene
gene1_control <- counts[1, 1:3]
gene1_treated <- counts[1, 4:6]
t.test(gene1_control, gene1_treated)

It gives me a p-value, sure. But here's the problem: RNA-Seq data isn't normally distributed. It's count data. It follows a Poisson-like distribution, but with "overdispersion"—meaning the variance is often much higher than the mean. If I use a t-test, I'm going to get a massive amount of false positives because the t-test assumes a normality that just isn't there in raw sequencing reads. I need a model that understands the nature of counts.

Letting DESeq2 handle the heavy lifting

This is where DESeq2 comes in. Instead of me trying to manually normalize the data or guess the distribution, DESeq2 uses a Negative Binomial distribution to model the counts. I'll start by wrapping my data into a DESeqDataSet object. This is the core data structure the package expects.

library(DESeq2)

# Creating the DESeqDataSet
dds <- DESeqDataSetFromMatrix(countData = counts, 
                              colData = metadata, 
                              design = ~ condition)

# Now I run the actual analysis pipeline
dds <- DESeq(dds)

That DESeq() function is actually doing three things under the hood: estimating the size factors (to account for the fact that some libraries were sequenced deeper than others), estimating the dispersion, and fitting the Negative Binomial GLM. It's a lot cleaner than trying to write a custom normalization loop.

The Multiple Testing Trap

Now I'll pull the results. I'm looking for genes that are significantly downregulated by the drug.

res <- results(dds)
head(res)

Looking at the output, I see a pvalue column and a padj column. If I just filter by pvalue < 0.05, I'll probably end up with 2,000 "significant" genes. But remember, I'm testing roughly 20,000 genes simultaneously. By pure chance, 5% of those (1,000 genes!) would show a p-value below 0.05 even if the drug did absolutely nothing. This is the classic multiple testing problem.

I have to use the padj (adjusted p-value), which uses the Benjamini-Hochberg correction. It's much more conservative, but it's the only way to ensure I'm not chasing ghosts.

Making the results actually legible

Staring at a table of 20,000 rows is a great way to get a headache. I need a way to see the "big picture"—which genes have the largest fold change and the highest significance? I'll build a volcano plot. I'll use ggplot2 for this because I want control over the aesthetics.

library(ggplot2)

res_df <- as.data.frame(res)

ggplot(res_df, aes(x = log2FoldChange, y = -log10(padj))) +
  geom_point(aes(color = padj < 0.05), alpha = 0.5) +
  theme_minimal() +
  labs(title = "Kinase Inhibitor Treatment", 
       x = "Log2 Fold Change", 
       y = "-log10 Adjusted P-value") +
  scale_color_manual(values = c("grey", "red"))

Now it clicks. The points in the top-left and top-right corners are my real winners. The ones in the middle are just noise. By moving from a simple t-test to a Negative Binomial model and correcting for multiple testing, I've gone from a list of 2,000 random guesses to a handful of high-confidence targets.

📋 Practical Task

Exercise: Filtering and Extracting High-Confidence Biomarkers

You have been provided with a DESeqResults object called res_kinase from a study on lung cancer cells. Your goal is to extract a specific list of candidate biomarkers for further validation.

Write a script that performs the following:

• Convert the res_kinase object into a standard R data frame.
• Filter the data frame to keep only the genes that meet both of these criteria:
  1. An adjusted p-value (padj) of less than 0.01.
  2. An absolute log2 fold change (abs(log2FoldChange)) of greater than 2.0.
• Sort the resulting list in descending order based on the absolute log2 fold change.
• Create a new data frame containing only the top 10 most significantly changed genes.

Expected Output: A data frame of 10 rows representing the most strongly differentially expressed genes in the dataset.