Bulk RNA-Seq: reads to differential expression

A complete bulk RNA-Seq workflow — quality control, adapter trimming, kallisto pseudoalignment, and differential expression with tximport and DESeq2.

Worked example Mouse — Mus musculus GRCm38

Design 6 samples, 2 conditions, single-end

Time to run A working day

Tools FastQC · Trim Galore · kallisto · tximport · DESeq2

This is a full bulk RNA-Seq workflow, split into three parts. It starts with the fastq.gz files your sequencing core hands you and finishes with a table of differentially expressed genes annotated with gene symbols and descriptions.

The worked example throughout is a real experiment: six mouse samples, three treated with the Eg5 inhibitor STLC and three vehicle controls (DMSO), sequenced single-end. Every output block you see is the actual result from that run, not an invented example.

What you will be able to do

• Judge whether your reads need trimming, and trim them if they do
• Pseudoalign reads to a transcriptome with kallisto, and explain what “pseudoalignment” means
• Summarise transcript-level abundance to gene level with tximport
• Fit a negative binomial model with DESeq2 and read its output
• Export an annotated results table your collaborators can open in Excel

Before you start

NotePrerequisites

You need a computer you can install software on, roughly 30 GB of free disk space, and 8 GB of RAM (16 GB is more comfortable). You do not need a compute cluster — kallisto is fast enough that this entire workflow runs on a laptop.

You do not need prior programming experience. You do need to type commands and read error messages.

The walkthroughs

Set up your machine Install everything in one sitting — WSL or Terminal, conda, R and RStudio, and every package the other two parts depend on. Do this first; it is the part that goes wrong. Raw reads to kallisto Quality control with FastQC, adapter and quality trimming with Trim Galore, then pseudoalignment against a kallisto index to produce per-sample abundance estimates. kallisto to DESeq2 Build a transcript-to-gene map, summarise abundances with tximport, fit the DESeq2 model, inspect the fit, and write out an annotated differential expression table.

Why kallisto rather than a genome aligner

A traditional workflow aligns every read to the genome with STAR or HISAT2, then counts reads per gene with featureCounts or HTSeq. For some questions — novel splice junctions, allele-specific expression, and anything where you need the actual alignment coordinates — it is the right tool.

kallisto pseudoaligns reads against a transcriptome For a straightforward differential expression analysis. Asking which transcripts a read is compatible with is far less work, so the whole quantification step runs in minutes on a laptop instead of hours on a cluster. Salmon works the same way and is an equally good choice.

A trade-off worth knowing: because you never produce a BAM file, if your analysis might need to go back and look at where a read landed, use a genome aligner instead.

Citing the tools

If you publish work using this pipeline, cite the tools:

• kallisto — Bray NL, Pimentel H, Melsted P, Pachter L. Near-optimal probabilistic RNA-seq quantification. Nature Biotechnology 34, 525–527 (2016).
• tximport — Soneson C, Love MI, Robinson MD. Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences. F1000Research 4:1521 (2015).
• DESeq2 — Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology 15:550 (2014).
• Trim Galore — Krueger F. Trim Galore. Babraham Bioinformatics.
• FastQC — Andrews S. FastQC: a quality control tool for high throughput sequence data. Babraham Bioinformatics.

Record the exact versions you used. conda list --export > environment.txt in your active environment and sessionInfo() in R will both do this for you, and your future self will thank you.