RNA-seq Differential Expression with DESeq2
RNA sequencing (RNA-seq) is the standard method for comparing gene expression between conditions. This guide covers everything in one place: what RNA-seq is and why we use it, how to design the experiment and prepare the materials, how to avoid and correct batch effects, which tools to choose and how they differ, how the DESeq2 design formula works, how to read PCA plots, log2 fold changes and adjusted p-values, and how to troubleshoot the problems you will meet along the way.
1. What is RNA-seq?
RNA-seq is a high-throughput sequencing method that measures the amount of RNA present in a biological sample at a given moment. RNA molecules are converted to complementary DNA (cDNA), fragmented, attached to adapters, and sequenced on a platform such as Illumina. The resulting millions of short reads are then mapped to a reference genome or transcriptome, and the number of reads assigned to each gene is used as a proxy for how strongly that gene is expressed.
The result is a count matrix: one row per gene, one column per sample, and each cell contains the number of reads that came from that gene in that sample. Everything that follows in a differential expression analysis (normalization, statistical testing, visualization) starts from this table.
A few definitions that appear throughout this guide:
- Transcriptome: the complete set of RNA molecules (mRNA, lncRNA, miRNA, rRNA and others) present in a cell or tissue.
- Bulk RNA-seq: the RNA of thousands to millions of cells is extracted together, so you measure the average expression of the whole sample. Single-cell RNA-seq, in contrast, measures each cell separately.
- Differential expression (DE): a statistical statement that a gene has a different expression level between two or more conditions by more than expected from random biological and technical variation.
- Library: the collection of cDNA fragments with adapters that is loaded on the sequencer. One library is prepared per sample.
2. Why and when to use RNA-seq (purpose and applications)
The main purpose of a bulk RNA-seq experiment is to answer the question: which genes change their expression between condition A and condition B, and by how much? Typical applications are:
- Comparing treated and control samples (drug, stress, infection, knockout, growth condition).
- Comparing tissues, developmental stages or disease and healthy individuals.
- Time-course experiments that follow expression after a stimulus.
- Discovery of biomarkers and of pathways affected by a condition.
- Annotation of genes and discovery of new transcripts, isoforms, fusion genes or splice variants (this needs deeper sequencing and different tools than simple DE).
- Studying non-model organisms, for example vector snails, where no microarray exists.
How does it compare with the methods it largely replaced?
| Feature | Microarray | RT-qPCR | RNA-seq |
|---|---|---|---|
| Genes measured | Only probes on the chip | A few to a few dozen | The whole transcriptome |
| Prior knowledge of sequences | Required | Required (primers) | Not required |
| Dynamic range | Limited (saturation) | Wide | Very wide |
| Novel transcripts and isoforms | No | No | Yes |
| Cost per sample | Low to medium | Very low for few genes | Medium |
| Typical use today | Legacy studies | Validation of selected genes | Discovery and genome-wide profiling |
RNA-seq is not always the right tool. If you only care about three genes in many samples, RT-qPCR is cheaper and faster. If you need expression in specific cell types inside a heterogeneous tissue, single-cell or spatial methods are better. Bulk RNA-seq of a mixed tissue can also be misleading when treatment changes the proportion of cell types, because the change looks like differential expression.
3. Experimental design: decide this before you extract any RNA
A poor design cannot be rescued by sophisticated statistics. Spend most of your thinking time here.
1. Define the biological question and the comparisons
Write the comparisons you want in a single sentence: for example, "genes differently expressed in Bulinus snails exposed to the parasite compared with unexposed snails at 24 hours". This tells you the factor of interest (exposure), the levels (exposed, unexposed), the other variables that must be controlled (age, sex, tissue, time) and the number of groups you need.
2. Biological replicates are non-negotiable
- Biological replicates are independent samples from different individuals, cultures or animals. They capture the real variation between individuals, and statistics cannot work without them.
- Technical replicates are the same RNA sample prepared or sequenced more than once. They measure only technical noise. They do not replace biological replicates, and for modern RNA-seq they are usually unnecessary.
How many replicates? Three per condition is the absolute minimum accepted by most reviewers and DE software. A large benchmarking study in yeast (Schurch et al., 2016) concluded that at least six biological replicates per condition are needed for reliable results, and at least twelve if you want to detect genes with small fold changes. When biological variability is high, such as wild-caught animals or human patients, more replicates are far more valuable than deeper sequencing.
3. Avoid confounding
A confounded design is one in which the factor of interest cannot be separated from another factor. The classic mistake: all controls are processed on Monday and all treated samples on Tuesday. Any difference between the groups could then be a processing-day effect, and no software can separate the two. The solution is to randomize and balance: every batch (day, kit, lane, operator) should contain samples from all conditions in similar numbers.
4. Sequencing depth, read length and read type
| Parameter | Typical recommendation for gene-level DE | Comment |
|---|---|---|
| Depth | 20 to 30 million reads per sample for human or mouse | More for isoform analysis or low-expressed genes, less for small genomes |
| Read length | 75 to 150 bp | Longer reads help with isoforms and repeats |
| Single-end vs paired-end | Single-end is enough for gene-level DE, paired-end is better for splicing and novel transcripts | Paired-end improves mapping accuracy |
| Strandedness | Use a stranded protocol | Resolves overlapping genes on opposite strands |
For gene-level DE the number of replicates matters more than depth. Beyond roughly 20 to 30 million reads, adding replicates gives more statistical power than adding reads.
5. Library preparation strategy
| Strategy | How it works | Best for | Limits |
|---|---|---|---|
| Poly(A) selection | Oligo-dT beads capture polyadenylated mRNA | Good-quality RNA from eukaryotes, mRNA profiling | Misses non-polyadenylated RNA, 3' bias if RNA is degraded |
| rRNA depletion | Probes remove ribosomal RNA, everything else stays | Degraded RNA (FFPE), bacteria, non-coding RNA, nascent RNA | More expensive, probes must match the species |
| 3' tag sequencing (for example QuantSeq) | One read per transcript at the 3' end | Large numbers of samples, low-cost DE | No isoform information |
For bacteria, which lack poly(A) tails, rRNA depletion is the only suitable option. For species without commercial probes (many invertebrates), poly(A) selection is usually the practical choice, provided the RNA is intact.
6. Controls, spike-ins and UMIs
- Spike-ins (for example ERCC) are synthetic RNAs added in known amounts. They help in special cases such as global changes in total RNA, but are not required for ordinary DE.
- UMIs (unique molecular identifiers) mark each original molecule so PCR duplicates can be removed. They are most useful with low input and in single-cell protocols.
- Include a negative control (no-template) during library preparation if contamination is a worry.
4. Materials and wet-lab workflow
Plan the experiment around this chain, because every step introduces variation that you can only control if you know about it.
- Sample collection and stabilization. Snap-freeze in liquid nitrogen or place the tissue in an RNA stabilizer (for example RNAlater) immediately. RNA degrades in minutes at room temperature. Record the time, collector, tissue and animal for every sample.
- RNA extraction. Use a guanidinium-phenol reagent (TRIzol or equivalent) or a silica-column kit. Use RNase-free tubes, tips and water, and clean the bench with an RNase decontaminant.
- DNase treatment. Genomic DNA carry-over produces reads in intergenic regions and inflates counts, so treat with DNase.
- RNA quality control. Check concentration with a fluorometer (Qubit), purity with a spectrophotometer (A260/280 close to 2.0 and A260/230 between 2.0 and 2.2) and integrity with a Bioanalyzer or TapeStation. The RIN (RNA integrity number) of 8 or higher is ideal, and 7 or higher is generally acceptable for poly(A) protocols. Some invertebrate samples show a single 28S/18S peak pattern that lowers the RIN artificially, so interpret the electropherogram, not only the number.
- Library preparation. Common kits include the Illumina Stranded mRNA Prep, NEBNext Ultra II Directional RNA and Illumina Stranded Total RNA with Ribo-Zero. Process samples in randomized order and in small groups containing all conditions.
- Library QC and pooling. Check fragment size distribution and quantify by qPCR or fluorometry, then pool at equal molarity.
- Sequencing. Illumina platforms (NovaSeq, NextSeq, MiSeq) are standard. Distribute each library across lanes so lane effects do not align with conditions.
colData (metadata) in DESeq2, and it is what lets you diagnose problems later.5. Batch effects: what they are and what to do about them
A batch effect is systematic variation caused by how, when or by whom samples were handled, rather than by biology. It is one of the most common reasons why an RNA-seq analysis gives disappointing or misleading results.
Common sources
- Different RNA extraction days or operators.
- Different library preparation kits, lots or dates.
- Different sequencing runs or lanes.
- Different RNA quality (RIN) between groups.
- Biological nuisance variables: sex, age, time of day, cage, weight.
Prevention is better than correction
- Randomize sample order at every wet-lab step.
- Make sure every batch contains all conditions.
- Extract and prepare as many samples as possible together.
- Record every variable that could matter.
Detection
- A PCA plot colored by batch, condition and any technical variable. If PC1 or PC2 separates samples by batch rather than condition, you have a batch effect.
- Hierarchical clustering or a sample-to-sample correlation heatmap.
- Check whether PCs correlate with metadata variables such as RIN or library size.
Correction: three strategies
- Include batch in the model (preferred). Add the batch as a covariate in the design formula, for example
~ batch + condition. The test for condition is then adjusted for batch differences, and raw counts are kept intact, which is what the statistical model needs. - Adjust counts with ComBat-seq. The
svapackage can adjust the count matrix before analysis. Use it mainly when you need adjusted counts for a downstream method that cannot take a covariate. - Remove batch for visualization only.
limma::removeBatchEffect()on variance-stabilized data is useful for PCA and heatmaps, but its output must never be used as input for DE tests.
# 1. Model (preferred): batch in the design formula
dds <- DESeqDataSetFromMatrix(counts, meta, design = ~ batch + condition)
# 2. ComBat-seq on counts (when the covariate approach is not possible)
library(sva)
adjusted <- ComBat_seq(as.matrix(counts), batch = meta$batch, group = meta$condition)
# 3. Visualization only
library(limma)
vsd <- vst(dds, blind = FALSE)
assay(vsd) <- removeBatchEffect(assay(vsd), batch = vsd$batch,
design = model.matrix(~ condition, colData(vsd)))
plotPCA(vsd, intgroup = c("condition", "batch"))If the sources of variation are unknown, surrogate variable analysis (sva) or RUVSeq can estimate hidden factors, but they should be used with care because they can also remove real biological signal.
6. From reads to counts: the analysis pipeline step by step
File formats you will meet
| File | Content | Produced by |
|---|---|---|
| FASTQ | Raw reads with base qualities | Sequencer |
| BAM | Reads aligned to the genome | STAR, HISAT2 |
| GTF or GFF | Gene annotation (coordinates of genes and exons) | Ensembl, GENCODE, NCBI |
| counts.txt | Gene counts per sample | featureCounts, STAR, HTSeq |
| quant.sf | Transcript abundances | Salmon |
Step 1: quality control and trimming
Run FastQC on each FASTQ file and combine the reports with MultiQC. Look at per-base quality (should be high, typically Q30 or above for most of the read), adapter content, overrepresented sequences, and GC content. Trim adapters and low-quality ends with fastp, Trimmomatic or cutadapt. Modern aligners handle soft clipping, so aggressive trimming is rarely needed.
Note that RNA-seq data normally shows high sequence duplication because highly expressed genes legitimately produce many identical reads. Do not remove duplicates based on position in standard bulk RNA-seq unless you have UMIs.
fastqc -t 8 raw/*.fastq.gz -o qc/
fastp -i R1.fastq.gz -I R2.fastq.gz -o R1.trim.fq.gz -O R2.trim.fq.gz \
--detect_adapter_for_pe --thread 8 -h fastp.html
multiqc qc/ -o qc/multiqcStep 2: choose a mapping and counting strategy
There are two families of approaches, and they differ in what they do with each read.
| Tool | Approach | Output | Strengths | Limits |
|---|---|---|---|---|
| STAR | Splice-aware alignment to the genome | BAM (and optional gene counts) | Fast, accurate, finds novel splice junctions, BAM useful for visualization and variant work | Needs lots of RAM (about 30 GB for human) |
| HISAT2 | Splice-aware alignment with a small memory footprint | BAM | Runs on modest computers | Slightly less sensitive than STAR for complex cases |
| featureCounts, HTSeq-count | Counts aligned reads overlapping genes | Gene count table | Simple and transparent | Ambiguous multi-mapping and overlapping reads are discarded or ignored |
| Salmon | Quasi-mapping of reads to the transcriptome, with bias correction | Transcript abundance (estimated counts and TPM) | Very fast, handles isoforms and multi-mapping reads, GC and sequence bias models | No BAM unless requested, needs a good transcript set |
| kallisto | Pseudoalignment to the transcriptome | Transcript abundance | Extremely fast and light | Similar limits to Salmon |
| RSEM | Alignment plus expectation-maximization | Gene and isoform estimates | Statistically rigorous isoform quantification | Slower |
How to choose. If you have a well-annotated organism and only need gene-level DE, Salmon followed by tximport is fast and accurate. If you also need BAM files (for visualization in IGV, variant calling, or because your genome annotation is incomplete, as in many non-model organisms), use STAR plus featureCounts. Both routes feed DESeq2 well, and with good annotation they give very similar DE results.
For non-model organisms with a draft genome, align to the genome with STAR or HISAT2 and count with your own GTF. A transcriptome-only approach depends entirely on the quality of the transcript set.
Step 3: check library strandedness
Most modern kits are stranded (dUTP protocols give "reverse stranded" reads, which is -s 2 in featureCounts). Using the wrong setting loses half of your counts or misassigns reads. Infer strandedness with RSeQC's infer_experiment.py, or let Salmon detect it with -l A.
Step 4: example commands
# STAR index (once per genome and annotation)
STAR --runMode genomeGenerate --genomeDir star_index \
--genomeFastaFiles genome.fa --sjdbGTFfile genes.gtf \
--sjdbOverhang 99 --runThreadN 16
# Align one sample
STAR --genomeDir star_index --readFilesIn R1.trim.fq.gz R2.trim.fq.gz \
--readFilesCommand zcat --runThreadN 16 \
--outSAMtype BAM SortedByCoordinate --outFileNamePrefix sample1_
# Count reads per gene for all samples together (reversely stranded paired-end)
featureCounts -T 8 -p --countReadPairs -s 2 -a genes.gtf -o counts.txt sample*_Aligned.sortedByCoord.out.bam
# Alternative: Salmon
salmon index -t gentrome.fa -d decoys.txt -i salmon_index -p 16
salmon quant -i salmon_index -l A -1 R1.trim.fq.gz -2 R2.trim.fq.gz \
-p 8 --gcBias --seqBias -o quant/sample1If you use STAR's --quantMode GeneCounts, the ReadsPerGene.out.tab file has four columns: gene, counts for unstranded data, counts when read 1 matches the gene strand, and counts when read 2 matches the gene strand. Pick the right column for your protocol.
A good outcome looks like this: uniquely mapped reads above 70 to 80 percent for a well-assembled genome, with most reads assigned to exons.
Step 5: build the count matrix in R
# From featureCounts
cnt <- read.delim("counts.txt", comment.char = "#", row.names = 1, check.names = FALSE)
counts <- cnt[, -(1:5)] # drop Chr, Start, End, Strand, Length
colnames(counts) <- sub("_Aligned.*", "", colnames(counts))
# From Salmon with tximport
library(tximport)
files <- file.path("quant", meta$sample, "quant.sf"); names(files) <- meta$sample
txi <- tximport(files, type = "salmon", tx2gene = tx2gene) # tx2gene: transcript ID, gene ID
dds <- DESeqDataSetFromTximport(txi, colData = meta, design = ~ batch + condition)7. Normalization and the many ways to express "expression"
Raw counts depend on sequencing depth and on the composition of the RNA, so samples cannot be compared directly. Different units solve different problems.
| Unit | What it corrects | Use for | Do not use for |
|---|---|---|---|
| Raw counts | Nothing | Input to DESeq2, edgeR, limma-voom | Plots comparing samples |
| CPM | Library size | Quick filtering, edgeR/limma | Comparing different genes |
| FPKM, RPKM | Library size and gene length | Legacy reports | Cross-sample comparison and DE |
| TPM | Library size and gene length, sums to one million per sample | Comparing genes within a sample, reporting abundance | DE tests (use counts) |
| DESeq2 normalized counts | Library size and composition (median of ratios) | Plotting a single gene across samples | Statistical testing again |
| vst or rlog values | Normalization plus variance stabilization | PCA, clustering, heatmaps | DE testing |
DESeq2 estimates a size factor for each sample with the median-of-ratios method, and edgeR uses TMM. Both assume that most genes are not differentially expressed. If a treatment changes the expression of a very large part of the transcriptome, this assumption is violated and you should look at the results with extra care.
8. Differential expression tools compared
| Tool | Statistical model | Normalization | Strengths | When to choose it |
|---|---|---|---|---|
| DESeq2 | Negative binomial GLM, dispersion shrunk toward a fitted trend, Wald or likelihood ratio test | Median of ratios | Robust with small samples, good defaults, built-in LFC shrinkage, outlier detection | Default choice for most projects with 3 to 20 replicates |
| edgeR | Negative binomial GLM, empirical Bayes dispersion, quasi-likelihood F test | TMM | Flexible, fast, very good with complex designs | Equivalent alternative, popular for large designs |
| limma-voom | Linear model on log-CPM with precision weights | TMM | Very fast, excellent for large sample numbers and complex designs | More than about 20 samples or many contrasts |
| NOISeq, others | Non-parametric | Various | Works without replicates | Exploratory only, not for formal inference |
On well-designed data these three main tools usually agree on the strongest signals, and the differences appear for borderline genes. Choosing one and using it correctly matters more than which one you pick.
9. DESeq2 in detail
How DESeq2 works in four sentences
DESeq2 models the counts of each gene as following a negative binomial distribution, which allows the variance to be larger than the mean (overdispersion), as is typical in biological data. For each gene it estimates a dispersion value, then shrinks those estimates toward a trend fitted across all genes so that genes with few replicates borrow information from the rest. It fits a generalized linear model whose terms come from the design formula, and tests whether the coefficient of interest differs from zero (Wald test) or whether dropping a term worsens the fit (likelihood ratio test). P-values are then corrected for multiple testing.
The design formula
The design formula tells DESeq2 which variables explain the counts. It is written with the metadata column names, and the variable of interest goes last.
| Experiment | Design formula | Meaning |
|---|---|---|
| Two groups | ~ condition | Compare treated and control |
| Known batch | ~ batch + condition | Test condition after removing batch differences |
| Paired samples (same individual before and after) | ~ subject + treatment | Each subject is its own control, which removes individual variation |
| Continuous covariate | ~ age + condition | Adjust for a numerical variable |
| Two factors, interaction | ~ genotype + treatment + genotype:treatment | Tests whether the treatment effect depends on genotype |
| Many groups, explicit contrasts | ~ 0 + group | No intercept, then define contrasts for each comparison |
| Time series | ~ batch + time, test with LRT | Genes changing at any time point |
Important rules:
- Variables in the formula must be columns of
colData, and categorical ones should be factors. - The reference level is the denominator of the fold change. Set it explicitly with
relevel(meta$condition, ref = "control"), otherwise R uses alphabetical order and your fold changes may be inverted. - The design must be full rank: no variable can be a perfect combination of others. Complete confounding of batch and condition breaks this.
- Do not add many covariates with few samples; each one costs a degree of freedom.
- Do not include variables that are consequences of the treatment.
# Interaction example: does the treatment effect differ between wild type and knockout?
dds <- DESeqDataSetFromMatrix(counts, meta, design = ~ genotype + treatment + genotype:treatment)
dds <- DESeq(dds)
resultsNames(dds)
res_interaction <- results(dds, name = "genotypeKO.treatmenttreated")
# Time series with the likelihood ratio test
dds <- DESeqDataSetFromMatrix(counts, meta, design = ~ batch + time)
dds <- DESeq(dds, test = "LRT", reduced = ~ batch)
res_time <- results(dds)
# Several groups: choose the comparison explicitly
res_AvsB <- results(dds, contrast = c("group", "A", "B")) # A relative to BThe full workflow
library(DESeq2)
stopifnot(all(colnames(counts) == rownames(meta))) # same samples, same order
meta$condition <- relevel(factor(meta$condition), ref = "control")
meta$batch <- factor(meta$batch)
dds <- DESeqDataSetFromMatrix(countData = round(counts), colData = meta,
design = ~ batch + condition)
# Pre-filter genes with almost no information
keep <- rowSums(counts(dds) >= 10) >= 3 # at least 10 reads in at least 3 samples
dds <- dds[keep, ]
dds <- DESeq(dds) # size factors, dispersion, model fit, Wald test
resultsNames(dds) # coefficient names
res <- results(dds, name = "condition_treated_vs_control", alpha = 0.05)
resLFC <- lfcShrink(dds, coef = "condition_treated_vs_control", type = "apeglm")
summary(res) # numbers of up- and down-regulated genes
out <- as.data.frame(resLFC)
out <- out[order(out$padj), ]
write.csv(out, "DE_treated_vs_control.csv")10. Reading the results table
Each gene gets a row with these columns:
| Column | Meaning |
|---|---|
baseMean | Mean of normalized counts across all samples. Low values mean little information |
log2FoldChange | Effect size: log2 of the ratio treated/control |
lfcSE | Standard error of the log2 fold change |
stat | Wald statistic (log2 fold change divided by its standard error) |
pvalue | Raw p-value for the test |
padj | P-value adjusted for multiple testing (Benjamini-Hochberg) |
What is log2 fold change?
The log2 fold change expresses expression differences on a scale where each unit is a doubling.
| log2FC | Fold change | Interpretation |
|---|---|---|
| 0 | 1 | No change |
| 1 | 2 | Twice as much in treated |
| 2 | 4 | Four times as much |
| -1 | 0.5 | Half as much (down-regulated) |
| -3 | 0.125 | Eight times lower |
The sign depends on the reference level: a positive value means higher in the numerator group (treated) than in the reference (control). A log2FC is only trustworthy when the gene has enough counts, because tiny counts produce enormous but meaningless ratios.
Why shrink fold changes?
For genes with few reads, the raw log2FC estimate is noisy and often exaggerated. lfcShrink() pulls unreliable estimates toward zero while leaving well-measured genes almost unchanged. Use the shrunken values for ranking, plotting (MA plots, volcano plots) and gene set analysis. The choices are:
apeglm: recommended default, works withcoef.ashr: works withcontrast, useful for comparisons that are not a model coefficient.normal: the older method, no longer recommended.
Shrinkage does not change the p-values or padj.
To require a minimum effect size, test it directly instead of filtering after the fact: results(dds, lfcThreshold = 1, altHypothesis = "greaterAbs") tests whether the absolute log2FC is greater than 1.
P-values and multiple testing correction
A p-value is the probability of observing a difference at least this large if the gene were not truly different. A p-value of 0.05 means that, if there were no true effect, about 5 percent of tests would still appear significant by chance.
In RNA-seq you run about 20,000 tests. At a raw threshold of 0.05, you would expect roughly 1,000 false positives even if nothing were changing. Multiple testing correction fixes this:
- Bonferroni multiplies p-values by the number of tests. It controls the chance of even one false positive, and is extremely conservative for genome-wide data.
- Benjamini-Hochberg (BH) controls the false discovery rate (FDR). A
padjof 0.05 means that among the genes called significant at this threshold, about 5 percent are expected to be false discoveries. This is the method DESeq2 uses.
Use padj, not pvalue, to decide which genes are significant. A common reporting rule is padj < 0.05 and |log2FC| > 1, but the fold-change threshold is a biological choice, not a statistical law.
Why are some padj values NA?
- All counts are zero in every sample for that gene.
- The gene has an extreme count outlier according to Cook's distance (the gene is flagged rather than tested).
- Independent filtering: DESeq2 automatically removes genes with very low mean counts, since they have no power, in order to increase the number of discoveries at the chosen
alpha.
Check the p-value histogram. It is a quick diagnostic of the test.
hist(res$pvalue[res$baseMean > 1], breaks = 0:20/20, col = "grey70",
main = "P-value histogram", xlab = "p-value")A healthy histogram is flat (uniform) with a spike near zero. A strongly U-shaped or bumpy histogram suggests an unmodeled batch effect or a problem with the model, and a nearly flat histogram without a spike means little real signal.
11. PCA and sample-level quality control
Before looking at any gene list, ask whether the samples behave as expected.
Principal component analysis (PCA)
PCA reduces thousands of genes to a few axes (principal components, PCs) that capture most of the variation between samples. PC1 explains the largest share, PC2 the next, and so on. Samples that are close on the plot have similar expression profiles overall.
vsd <- vst(dds, blind = TRUE) # blind = TRUE for unbiased QC
plotPCA(vsd, intgroup = c("condition", "batch"))How to read a PCA:
- Replicates of the same condition should cluster together.
- If the condition separates along PC1 or PC2, the biology dominates, which is what you hope for.
- If batch, sex, RIN or any technical variable separates samples, you need to model it.
- A single sample far away from its group is an outlier. Check its QC metrics (mapping rate, RIN, library size) before deciding to remove it, and always report the removal.
- The percentages on the axes show how much variance each PC explains; low percentages mean many small effects and no dominant signal.
Other QC plots
# Sample-to-sample distances
library(pheatmap)
d <- as.matrix(dist(t(assay(vsd))))
pheatmap(d, annotation_col = as.data.frame(colData(dds)[, c("condition", "batch")]))
plotDispEsts(dds) # dispersion should decrease with mean and fit the trend
plotMA(resLFC, ylim = c(-4, 4))In the dispersion plot, the points should scatter around the fitted red curve, and the curve should fall as the mean increases. Genes far above the curve (blue circles) are treated as outliers. In the MA plot, most points should be centered on zero; if the cloud is shifted up or down, a normalization problem is likely.
12. Visualizing differential expression
Volcano plot
A volcano plot shows effect size on the x axis and statistical significance on the y axis. It lets you see at a glance which genes are both large and significant.
library(ggplot2)
df <- as.data.frame(resLFC)
df$sig <- ifelse(df$padj < 0.05 & abs(df$log2FoldChange) > 1,
ifelse(df$log2FoldChange > 0, "Up", "Down"), "NS")
ggplot(df, aes(log2FoldChange, -log10(padj), colour = sig)) +
geom_point(alpha = 0.6) +
scale_colour_manual(values = c(Down = "#2166ac", NS = "grey70", Up = "#c0392b")) +
geom_vline(xintercept = c(-1, 1), linetype = "dashed") +
geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
theme_bw()Heatmap of top genes
library(pheatmap)
top <- rownames(df)[order(df$padj)][1:30]
ann <- as.data.frame(colData(dds)[, c("condition", "batch")]) # rows must match sample names
pheatmap(assay(vsd)[top, ], scale = "row", annotation_col = ann)scale = "row" converts each gene to a z-score so you compare patterns rather than absolute levels.
Plot a single gene
plotCounts(dds, gene = "GENE_ID", intgroup = "condition")13. Functional interpretation: from gene lists to biology
A list of thousands of genes is not a result. Test whether the changed genes concentrate in known pathways or functions.
- Over-representation analysis (ORA): takes the significant genes and asks whether a GO term or pathway contains more of them than expected by chance. The background (universe) must be all genes that could have been detected, not the whole genome.
- Gene set enrichment analysis (GSEA): uses all genes ranked by a statistic, so it does not need an arbitrary cutoff and can detect coordinated small changes.
library(clusterProfiler); library(org.Hs.eg.db) # replace with the organism database
sig <- rownames(subset(df, padj < 0.05 & abs(log2FoldChange) > 1))
universe <- rownames(df)[!is.na(df$padj)]
ego <- enrichGO(gene = sig, universe = universe, OrgDb = org.Hs.eg.db,
keyType = "ENSEMBL", ont = "BP", pAdjustMethod = "BH", qvalueCutoff = 0.05)
dotplot(ego, showCategory = 15)
# GSEA with a ranked list
gl <- res$stat; names(gl) <- rownames(res)
gl <- sort(gl[!is.na(gl)], decreasing = TRUE)
gse <- gseGO(geneList = gl, OrgDb = org.Hs.eg.db, keyType = "ENSEMBL", ont = "BP")For non-model organisms such as Bulinus or Biomphalaria, annotate the genes with orthologs or domain databases (eggNOG-mapper, InterProScan), build a custom TERM2GENE table, and use enricher() or GSEA() from clusterProfiler.
14. Troubleshooting: symptoms, causes and fixes
| Symptom or error | Likely cause | What to do |
|---|---|---|
some values in assay are not integers | Counts come from TPM, FPKM or tximport estimates | Use raw counts, round() them, or use DESeqDataSetFromTximport |
the model matrix is not full rank | A variable is confounded with another (for example batch with condition) | Check table(meta$batch, meta$condition), remove or combine variables, or redesign |
colnames of countData do not match rownames of colData | Different sample order or names | Reorder with counts <- counts[, rownames(meta)] and check with all(colnames(counts) == rownames(meta)) |
| Fold changes have the wrong sign | Wrong reference level | Use relevel() before DESeq() and confirm with resultsNames(dds) |
| Very few reads assigned to genes (featureCounts) | Wrong strandedness (-s), annotation does not match the genome, or wrong chromosome naming | Test -s 0, -s 1, -s 2; use the GTF from the same assembly source as the FASTA |
| Low mapping rate (below 50 to 60 percent) | Contamination, wrong species or genome, adapters, degraded RNA | Check adapters, BLAST a subset of unmapped reads, check rRNA content |
| High rRNA fraction | Poor depletion or poly(A) selection | Run SortMeRNA, consider that samples may need to be resequenced |
| Gene IDs not matching in enrichment | ID versions (ENSG00000123.4) or wrong key type | Strip versions with sub("\\..*", "", ids) and use the correct keyType |
| No or very few DE genes | Too few replicates, high variability, weak biology, or confounding | Check PCA and dispersion, increase replicates, verify sample labels |
| Thousands of DE genes and a flat PCA by condition | Strong global effect, or batch or sample swap | Check size factors, sample labels and composition effects |
All padj are NA for many genes | Low counts, independent filtering, or Cook's outliers | Check baseMean, use results(dds, cooksCutoff = FALSE) carefully |
| Huge log2FC with tiny counts | Noisy estimates at low expression | Use lfcShrink, filter low baseMean, check raw counts of the gene |
| A single sample behaves differently | Poor RNA quality, sample swap, contamination | Inspect QC, check sex or marker genes, and justify any removal |
| Out of memory in STAR or R | Genome too large for the machine | Use --genomeSAindexNbases for small genomes, use HISAT2 or Salmon, increase RAM |
pheatmap annotation error | Annotation rownames do not match matrix colnames | Make sure rownames(ann) == colnames(matrix) |
A practical order for debugging. First check the sample sheet (labels, factors, reference levels). Second check the QC of reads and mapping. Third look at PCA and sample distances. Only then question the statistical model.
15. Reproducibility and reporting
- Record software versions with
sessionInfo()and keep scripts under version control. - Use workflow managers such as Snakemake or Nextflow. The community pipeline nf-core/rnaseq automates QC, alignment, quantification and reports in a standard way.
- Report the number of biological replicates, depth, kit, mapping rates, normalization, the full design formula, thresholds used, and the genome and annotation versions.
- Deposit raw reads and the count matrix in public repositories (GEO and SRA) and cite the accession.
16. Best-practice checklist
- Biological replicates (at least 3, ideally 5 to 6 or more) per condition.
- Samples randomized and balanced across all batches.
- RNA integrity checked and the same extraction method for all samples.
- Stranded library with a protocol suited to the organism (poly(A) or rRNA depletion).
- QC of reads (FastQC, MultiQC) and of mapping rates.
- Raw counts as input; low-count genes filtered.
- Known technical variables included in the design, variable of interest last.
- PCA and sample distances inspected before any DE interpretation.
- Decisions based on
padj, with shrunken log2 fold changes for ranking and plots. - Results validated by an independent method (RT-qPCR on selected genes) when important.
17. Frequently asked questions
Can I run RNA-seq without replicates? You can sequence it, but you cannot assess biological variability, so any p-value is unreliable. Treat such data as exploratory.
TPM or counts for DESeq2? Always counts. TPM is a good unit for describing abundance but violates the assumptions of count-based models.
Should I remove duplicates? Not in standard bulk RNA-seq without UMIs, because duplicates mostly reflect highly expressed genes.
DESeq2, edgeR or limma-voom? All three are reliable. DESeq2 is a good default for small experiments, and limma-voom scales well for large, complex designs.
What padj threshold should I use? 0.05 is conventional, 0.1 is the DESeq2 default for alpha. State your choice and apply it before looking at the genes.
How do I deal with a sample that looks like an outlier? Investigate the cause with QC metrics. Remove it only with a documented technical reason, and show results with and without it.
18. References
- Love MI, Huber W, Anders S (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology 15:550.
- Robinson MD, McCarthy DJ, Smyth GK (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26:139-140.
- Ritchie ME et al. (2015). limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Research 43:e47.
- Dobin A et al. (2013). STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29:15-21.
- Patro R et al. (2017). Salmon provides fast and bias-aware quantification of transcript expression. Nature Methods 14:417-419.
- Zhu A, Ibrahim JG, Love MI (2019). Heavy-tailed prior distributions for sequence count data: removing the noise and preserving large differences. Bioinformatics 35:2084-2092.
- Benjamini Y, Hochberg Y (1995). Controlling the false discovery rate. Journal of the Royal Statistical Society B 57:289-300.
- Schurch NJ et al. (2016). How many biological replicates are needed in an RNA-seq experiment and which differential expression tool should you use? RNA 22:839-851.
- Conesa A et al. (2016). A survey of best practices for RNA-seq data analysis. Genome Biology 17:13.
- Zhang Y, Parmigiani G, Johnson WE (2020). ComBat-seq: batch effect adjustment for RNA-seq count data. NAR Genomics and Bioinformatics 2:lqaa078.
💬 Commentaires