Differential Expression: A First Look

Differential expression is the study of which genes change their activity between two situations. A gene is expressed when the cell reads it and makes RNA copies from it, and we can count those RNA copies. Differential expression analysis asks a simple question: for each gene, is that count meaningfully higher or lower in one group of samples than in another?

The two groups being compared are called conditions. A common setup is control samples (untreated) versus treated samples (given a drug, for example). For every gene we have a read count, which is the number of sequencing reads that landed on that gene in a given sample. We want to find the genes whose counts shift between the control and treated conditions.

The DESeq2 Workflow

DESeq2 is a widely used R package (a bundle of ready-made functions) for this job. The whole analysis is just three steps. DESeqDataSetFromMatrix() packages your data into one object, DESeq() does the number-crunching by adjusting for differences in sequencing depth and fitting a small statistical model for each gene, and results() hands back a table of answers. You do not need to understand the math yet, only the shape of the workflow.

# DESeq2 is a Bioconductor package for finding differentially expressed genes.
library(DESeq2)

# counts: a table of raw read counts, one row per gene, one column per sample.
counts <- read.csv("gene_counts.csv", row.names = 1)

# coldata: describes each sample, here which experimental group it belongs to.
coldata <- data.frame(
  condition = c("control", "control", "control", "treated", "treated", "treated"),
  row.names = colnames(counts)
)

# Bundle the counts and the sample info into one object DESeq2 understands.
# design = ~ condition tells DESeq2 to compare genes across the condition groups.
dds <- DESeqDataSetFromMatrix(
  countData = counts,
  colData   = coldata,
  design    = ~ condition
)

# DESeq() normalizes for sequencing depth and fits a model for every gene.
dds <- DESeq(dds)

# results() extracts the per-gene table of effect sizes and p-values.
res <- results(dds)

# head() prints just the first 6 rows so we can read them.
head(res)

log2 fold change (MLE): condition treated vs control
Wald test p-value: condition treated vs control
DataFrame with 6 rows and 6 columns
         baseMean log2FoldChange     lfcSE      stat    pvalue      padj
        <numeric>      <numeric> <numeric> <numeric> <numeric> <numeric>
GENE1    1245.32          0.021     0.128     0.164   0.87012   0.94210
GENE2     532.10          2.415     0.311     7.765  8.16e-15  3.02e-13
GENE3      89.44         -1.874     0.402    -4.662  3.13e-06  4.55e-05
GENE4    2210.87          0.334     0.155     2.155   0.03115   0.08842
GENE5      12.03          0.512     0.889     0.576   0.56480   0.71223
GENE6     745.66         -3.201     0.298   -10.741  6.55e-27  9.10e-25

Two columns carry most of the meaning. The log2FoldChange is the change in expression on a log-base-2 scale: a value of 1 means the gene roughly doubled in the treated group, minus 1 means it roughly halved, and a value near 0 means little change. The padj is the adjusted p-value, a p-value corrected for the fact that we tested thousands of genes at once; a common rule is to trust genes with padj below 0.05. So GENE2 and GENE6 look strongly and reliably changed, while GENE1 does not change at all.

What A Volcano Plot Shows

A volcano plot is the standard way to view the whole results table in one picture. Each dot is one gene. The horizontal axis is the log2FoldChange, so genes drift left or right depending on whether they went down or up. The vertical axis is minus log10 of the adjusted p-value, a transform that pushes the most significant genes (tiny padj) toward the top. The genes worth chasing sit in the upper corners: far from the center because they changed a lot, and high up because that change is statistically solid.