Your first differential expression analysis, step by step

This is a walk through a standard bulk RNA-seq differential expression analysis, written for someone who generated the data and wants to understand what happens to it. It names the tools a core facility would use, explains what each step is for, and says what to look at before moving on.

What is differential expression analysis?

Differential expression analysis asks, for every gene, whether its expression differs between conditions by more than you would expect from ordinary variation between samples. The input is a count of sequencing reads per gene for each sample; the output is a table with one row per gene, giving the size of the change (a log2 fold change) and how confident the call is (an adjusted p-value).

The confidence comes from replicates. Two samples will always differ a little, so a model needs several independent samples per condition to estimate how much genes vary on their own. That is why the number of biological replicates matters more than almost any setting in the software.

Before you start: replicates and design

Plan for at least three biological replicates per condition, and more if small changes matter: Schurch and colleagues found that three replicates mainly detect the larger changes, and recommend at least six where finding most of the changed genes is the goal. Technical replicates of the same library do not count, because they do not capture the variation between independent samples.

Write down anything that differs between samples apart from the condition: extraction date, operator, sequencing run, library kit. These are your batches, and they belong in the model. A sample sheet with a column for each is the most useful thing you can bring to the analysis.

1. Check and trim the reads

Start with read quality. fastp trims adapter sequence and low-quality bases and writes a report for every sample. Read it. You are looking for a sharp drop in quality along the read, adapter content that was not removed, or a sample whose numbers look unlike the others.

High duplication is normal in RNA-seq, because abundant transcripts are sequenced many times over; it is not, on its own, a sign of a bad library.

Worth changing: almost nothing, at first. The defaults are sensible, and adapter detection for paired-end data finds the adapter from the read overlap, so you do not need to know the kit.

2. Count reads per gene

There are two established routes. STAR aligns reads to the genome, and reads are then counted per gene. Salmon skips the full alignment and estimates transcript abundance directly, then sums transcripts to genes. For gene-level differential expression both work well. Salmon is faster and needs far less memory, which matters on a laptop; STAR is the choice if you also want to look at reads in a genome browser, find new splice junctions or call variants from RNA.

Check: the mapping rate for each sample. For a well-annotated organism, most samples map above about 70%. One sample well below the rest usually points to contamination or degraded RNA; every sample being low usually means the wrong reference or annotation.

3. Test each gene with DESeq2

DESeq2 takes the counts and your sample sheet and fits a model to each gene. It estimates how much each gene varies between replicates, borrowing information across genes so that the estimate is stable even with few samples, and then tests whether the difference between conditions is larger than that variation.

Because it tests every gene at once, some genes will look significant by chance. DESeq2 adjusts the p-values with the Benjamini–Hochberg method, so calling genes significant at an adjusted p-value of 0.05 means about 5% of those calls are expected to be false discoveries.

Worth changing: the design formula. If samples were processed in two batches, write ~ batch + condition, so the model accounts for the batch before it tests the condition. Leave the batch out and a batch difference can be read as biology. Also set the contrast explicitly: treated vs control means a positive fold change is higher in the treated samples.

4. Look at the results

Make the PCA plot first. It places each sample by its overall expression profile, and replicates of the same condition should sit together. If the main split is by sequencing run rather than condition, go back and add the run to the design. If one replicate sits far from its group, find out why before reading anything into the gene list.

Then the volcano plot: log2 fold change on one axis, significance on the other, so the genes that changed a lot and confidently sit in the two upper corners. A heatmap of the top genes across samples shows whether the change is consistent across replicates or driven by one sample.

A useful sanity check: genes you expect not to change, such as common housekeeping genes, should not be at the top of the list.

What to take to lab meeting

• The QC report, with anything unusual marked.
• The PCA plot, and a sentence on how the samples group.
• The volcano plot and the top genes by adjusted p-value.
• The design formula and the settings you used, so someone can ask about them.

Running it without code

In OpenBionode this is the RNA-seq starter pipeline: the steps above, already connected, with each setting explained in the step’s panel. The bulk RNA-seq protocol goes through it setting by setting, with a troubleshooting table, and the RNA-seq use case shows the pipeline and what it produces.