RNA Seq Differential Expression Analysis

RNA sequencing differential expression (RNA seq DE) analysis is the computational process of identifying genes whose expression levels change significantly between experimental conditions, such as treated versus control or disease versus healthy. This guide is for bench biologists, bioinformatics newcomers, and experienced researchers transitioning to RNA seq, who need a practical, source anchored framework to design, execute, and interpret a differential expression study without getting lost in statistical jargon.

RNA seq DE analysis starts with raw sequencing reads and ends with a list of differentially expressed genes (DEGs) and their biological context. The core steps include read alignment or pseudo alignment, quantification, normalization, statistical modeling, and multiple testing correction. Every decision point, from choice of aligner to statistical method, affects the final gene list. Reliable results depend on careful experimental design, robust quality control, and awareness of the method's limitations NCBI Bookshelf.

This guide assumes you have access to raw FASTQ files from a sequencing facility or public repository NCBI Sequence Read Archive. You will also need a computational environment with R/Bioconductor or a Galaxy instance Galaxy Training Network. The framework below is tool agnostic but highlights two widely used Bioconductor packages, edgeR and DESeq2, which are compared in a recent benchmark Tool choice matters: evaluating edgeR vs. DESeq2 for sensitivity, robustness, and cross-study performance.

At a Glance

Decision Criteria

Before running any command, you must decide on a statistical framework. The choice depends on your sample size, expected effect sizes, and tolerance for false positives.

Sample size and replicates. For at least three biological replicates per condition, DESeq2 and edgeR perform similarly. With two replicates, use edgeR with a robust dispersion estimate or DESeq2 with the default method. A recent evaluation found that no single tool dominates across all scenarios, edgeR gives higher sensitivity in some settings, while DESeq2 offers better control of false positives with low counts Tool choice matters: evaluating edgeR vs. DESeq2 for sensitivity, robustness, and cross study performance.

Normalization method. Both packages apply internal normalization. DESeq2 uses median of ratios. edgeR uses TMM (trimmed mean of M values). For samples with highly variable library sizes or strong contamination, TMM is often more robust EMBL EBI Training.

Gene level counts vs. transcript level. If you want isoform specific expression, use a transcript quantifier like Salmon or kallisto and then import counts into a DE tool that supports tximport. For gene level analysis, align reads with STAR or HISAT2 and summarize with featureCounts.

Batch effects. Always check whether samples cluster by batch (e.g., sequencing run, extraction date). DESeq2 includes a ~ batch + condition design formula. edgeR uses model.matrix with batch as a covariate. If batch effects are severe, consider limma which handles complex designs.

Small numbers of DEGs. For experiments expecting few DE genes (e.g., nearly isogenic organisms), use a moderated test in limma trend or edgeR's exact test.

Practical Workflow or Implementation Sequence

Step 1: Raw Data Quality Control

Run FastQC on all FASTQ files. Inspect per base quality, GC content, adapters, and duplication levels. Use MultiQC to aggregate reports. Remove low quality bases and adapters with Trimmomatic, Cutadapt, or fastp.

Step 2: Alignment or Quantification

For alignment based analysis, map reads to a reference genome with STAR or HISAT2. STAR is fast and accurate for mammalian genomes. Use default parameters except for outFilterMismatchNmax (set to 10 for 100 bp reads). For alignment free quantification, use Salmon or kallisto in mapping based mode. These tools are faster and directly output transcript abundance Galaxy Training Network.

Step 3: Generating the Count Matrix

From alignment output, use featureCounts (part of the subread package) to assign reads to genes. For alignment free methods, use tximport to convert transcript level estimates to gene level counts.

Step 4: Quality Assessment of Counts

Filter lowly expressed genes. A common rule: keep genes with at least 10 counts in at least two samples. Then inspect sample clustering with principal component analysis (PCA) or multidimensional scaling (MDS). Look for outliers. If one sample separates entirely, consider removing it or checking for contamination. A study on Senegalese sole used PCA to confirm distinct expression patterns in olfactory rosettes Specialization of the upper and lower olfactory rosettes reveals sex and origin dependent plasticity in gene expression in Senegalese sole.

Step 5: Normalization and Dispersion Estimation

In DESeq2, the estimateSizeFactors function computes median of ratios. In edgeR, calcNormFactors computes TMM. Then estimate gene wise dispersions. DESeq2 shrinks dispersion estimates toward a fit, which is beneficial for small sample sizes. edgeR uses a similar empirical Bayes approach.

Step 6: Differential Expression Testing

Fit a model using the design matrix. For DESeq2: DESeq(dds). For edgeR: estimateGLMCommonDisp, estimateGLMTagwiseDisp, then glmQLFit and glmQLFTest. Extract results with a log2 fold change threshold (e.g., 1 or 1.5) and adjusted p value < 0.05.

Step 7: Multiple Testing Correction

The output contains raw p values. Apply Benjamini Hochberg FDR correction. Report the number of DEGs at a given FDR threshold.

Step 8: Downstream Interpretation

Annotate DEGs with gene symbols and functional databases (GO, KEGG, Reactome). Validate a subset with qPCR or independent cohorts. For example, a diabetes study validated GYS2 downregulation by integrating RNA seq with metabolomics Glycogen metabolic dysfunction in T2DM with MASLD: linking alpha hydroxybutyrate to GYS2 downregulation.

Common Mistakes

Insufficient replication. Using only two biological replicates per condition leads to unreliable dispersion estimates and high false positive rates. At least three replicates are minimum, six are better.

Ignoring batch effects. Failing to include batch in the design matrix can confound biological results. Always plot PCA colored by batch before testing.

Using raw counts in downstream tools. Do not feed raw counts into clustering or machine learning without proper normalization (e.g., variance stabilizing transformation, regularized log).

Overinterpreting small fold changes. A statistically significant DEG may have a log2 fold change of 0.5. Decide a biological meaningful threshold beforehand.

Not filtering low counts. Including genes with very low counts increases multiple testing burden and can cause false positives from technical noise.

Forgetting to check alignment artifacts. rRNA contamination, adapter dimers, or mapping to mitochondrial genes can distort counts. Remove mitochondrial reads if they dominate.

Limits of Interpretation

Statistical significance does not equal biological importance. A gene can be highly significant with a tiny fold change that has no functional effect. Conversely, a gene with large fold change may fail significance due to high variance.

DE analysis is correlational. It identifies changes associated with a condition, not necessarily causal. Validation in independent samples or with perturbation experiments is required.

Small sample sizes inflate false discovery. The Benjamini Hochberg procedure assumes independent tests, which is reasonable for genes, but with only two replicates per group, the actual FDR may be higher than reported [EMBL EBI Training](https://www.ebi.ac.uk/training/online/courses/functional genomics/chapter-5-rna-seq-experiment/).

Results depend on the reference genome and annotation. Different versions can produce different gene lists. Always document the reference genome build and annotation version.

Interpreting fold change direction depends on the reference condition. If you treat control as the baseline, a log2 fold change of 2 means the gene is upregulated in the treatment group. Reversing the contrast flips signs.

DE methods assume a negative binomial model. This model works well for most RNA seq data, but it can fail for highly expressed outliers or for data with excessive zeros (e.g., single cell RNA seq requires different methods).

Frequently Asked Questions

1. What is the minimum number of replicates for a differential expression experiment?

At least three biological replicates per condition are recommended. With two replicates, you must use a method that shares information across genes (e.g., edgeR with a robust dispersion estimate) but power is low. Two technical replicates are not a substitute for biological replication.

2. Should I use fold change or p value to filter DEGs?

Both are important. A common threshold is adjusted p value < 0.05 and absolute log2 fold change > 1 (two fold). However, for exploratory analysis you might use a relaxed fold change with a stringent p value. Always report both.

3. How do I deal with samples that have very low total read counts?

Low total counts ( < 5 million reads) indicate poor sequencing depth or degradation. Remove such samples if they are outliers. If all samples have low depth, you may need more sequencing. Normalization cannot fix missing data.

4. Can I reuse the same control group for multiple treatment groups?

Technically yes, but it violates the assumption of independent samples. This practice increases the chance of false positives because the control group is used multiple times. Instead, run all samples in one experiment or use a limma contrast approach with a proper design matrix.

References and Further Reading

• NCBI Bookshelf Free biomedical books covering RNA seq methodology.
• EMBL EBI Training Official training resources for RNA seq data analysis.
• Galaxy Training Network Hands on tutorials for RNA seq with Galaxy.
• Bioconductor Open source software for genomics, including DESeq2 and edgeR manuals.
• NCBI Sequence Read Archive Public repository for RNA seq data.
• Tool choice matters: Evaluating edgeR vs. DESeq2 for sensitivity, robustness, and cross-study performance PLoS One 2025 benchmark.
• Specialization of the upper and lower olfactory rosettes reveals sex and origin dependent plasticity in gene expression in Senegalese sole BMC Biol 2025 example of PCA usage.
• Construction and multi omics validation of a five gene ferroptosis based model for predicting prognosis and therapy response in pancreatic ductal adenocarcinoma Discov Oncol 2025 case study.
• Glycogen metabolic dysfunction in T2DM with MASLD linking alpha hydroxybutyrate to GYS2 downregulation Front Nutr 2025 validation example.

Related Articles

• Somatic Cell
• Phases In Cell Cycle
• Endoplasmic Reticulum Cell Function
• Stages Of Cell Cycles
• Pcr Test