Bulk RNA-Seq Data Analysis: From FASTQ to a Trustworthy Differential-Expression Result

Quick answer: Bulk RNA-seq analysis has one core deliverable — a table of genes whose expression differs between conditions, with an effect size and an adjusted p-value you can defend. The pipeline is: QC the reads, align and count (or quantify transcripts), normalize counts, model with a negative-binomial or voom framework (DESeq2, edgeR, or limma-voom), apply multiple-testing correction, and interpret fold changes with shrinkage. The modeling step matters less than the two things beginners skip most: count provenance and experimental design, which silently decide whether the statistics downstream mean anything.

Natural next step after the single-cell vs bulk comparison and the FASTA/FASTQ guides on this blog: bulk RNA-seq averages expression across every cell in the sampled tissue — this article covers what to do with those averaged read counts, FASTQ in, gene table out.

The pipeline at a glance

• Raw data QC — decide whether the reads are usable before anything else.
• Alignment or pseudoalignment + counting — turn reads into a gene-by-sample count matrix.
• Experiment sanity checks — PCA/heatmap of the counts against your design: do samples group by biology, or by batch?
• Normalization + statistical modeling — estimate library sizes, dispersion, and per-gene effects.
• Multiple-testing correction + interpretation — adjusted p-values, log2 fold changes, shrinkage, and honest reporting.

These stages are durable; the tool versions and defaults change. Date your exact commands and versions when you run them.

Step 1: raw-data QC — the cheap decisions

Before counting even begins, tools such as FastQC (or the QC steps built into any pipeline) answer a short list of questions:

• Per-base quality: do reads hold up across their full length, or does the tail drop? Affects whether trimming is required.
• Adapter content: are adapter sequences present in reads? (Trim.)
• Overrepresented sequences: a dominant rRNA, adapter, or anomalous transcript changes what you do next.
• GC distribution: strong skew may indicate contamination or biased transcript coverage.

This step is minutes of compute per sample. Skipping it means problems get discovered two weeks later in Step 4, where you can no longer cleanly tell biology from artifact.

Step 2: reads to counts — the two legitimate routes

There are two broad quantification routes and both are tools-of-the-trade, not competing dogmas:

• Genome alignment + gene counting. Align reads to the genome (STAR, HISAT2) and count reads per gene model (featureCounts). Transparent: you can audit every alignment. Costs more compute.
• Transcript quantification (pseudoalignment). Salmon or kallisto assigns reads to transcripts directly against a transcriptome index, producing transcript-level estimates quickly.

Whichever you choose, keep the same gene annotations (a database version for transcriptome or genome builds) documented with your run, and be aware that count matrices built from different routes are not freely interchangeable down-stream, because a matrix from a different annotation is a slightly different measurement.

The indexing choice within the count step matters more than the aligner: “counts per gene” versus “estimates per transcript.” If you plan to feed a negative-binomial model, deliver it integer counts (gene-level), since count models expect counts, not TPM/FPKM values — this is the most common input formatting error in beginner DE analyses.

Step 3: sample checks — design against the biology

Before any differential test, check the count matrix against your design:

• PCA or MDS plot: samples should cluster by biological condition. Clusters by batch, treatment date, or library prep batch indicate something else is driving variation — add it to the design formula.
• Sample correlation heatmap: a run-away outlier should be stared at, investigated, and documented when removed (the removal should be defensible rather than convenience).
• Replicates at the design level (biological replication, not technical) is the single largest determinant of statistical power — no analysis method rescues an n = 1 design.

This is also where a design formula is chosen: ~ batch + condition when batch overlaps the question, ~ batch + sex + condition when a factor must remain controlled. Design choices drive results far more than the DE tool’s brand.

Step 4: normalization and modeling

Raw counts are not yet comparable: library sizes differ by orders of magnitude, and long or highly expressed genes absorb reads in ways that mislead naive comparisons.

What normalization must do

• Scale counts to account for library size (DESeq2’s median-of-ratios and edgeR’s TMM both handle this without assuming every gene is unchanged; they effectively use the genes not expected to change as anchors).
• Account for RNA composition shifts (a handful of very highly expressed genes can soak up a majority of a lane’s reads).

The statistical reality of RNA-seq counts

Counts are discrete and their variance depends on the mean, so a simple t-test on raw or log counts fails models this structure. The three established frameworks:

• DESeq2 — negative-binomial generalized linear models with empirical-Bayes dispersion estimation, a mean-variance trend, and a design formula carried through; takes raw integer counts directly.
• edgeR — negative-binomial with empirical-Bayes dispersion shrinkage; likewise takes raw counts.
• limma-voom — linear modeling framework with per-observation precision weights derived from the mean-variance relationship; notably fast on larger designs and familiar to users of limma microarray workflows.

These are established alternatives rather than a universal podium. DESeq2 and edgeR differ most in dispersion handling details and defaults; limma-voom wins on speed with large designs; all three, well-used, give concordant sets on moderately powered data. Choose one, learn its assumptions (integer input, design formula, correction method), and record the version you used.

Fold changes: shrinkage is not normalization

A raw log2 fold change for a low-count gene swings wildly on one or two reads. Methods like DESeq2’s lfcShrink (with the apeglm or similar options) pull noisy fold changes toward a central tendency, which mainly helps ranking and visualization of low-expression genes. Shrinkage is a stabilization of effect-size estimates, not an additional normalization step, and it does not change which genes pass the adjusted-p threshold.

Step 5: multiple testing and honest interpretation

A genome-wide screen tests 20,000+ genes simultaneously, so raw p-values are meaningless without correction (Benjamini–Hochberg FDR is standard). Then interpret:

• Adjusted p-value < 0.05 with a practical fold-change cutoff like |log2FC| > 1 is a typical (never universal) threshold pair; the choice should be justified per study.
• “Not significant” is not “no change”: absent power, small true effects will not clear FDR; conclusions about “no difference” require the design to have had that power.
• Look at MA/volcano plots and the top hits in a genome browser rather than trusting sorted tables alone: a handful of high-fold-change hits with tiny counts behave differently than a broad set of moderate, well-supported hits.

Bulk vs single-cell, in one paragraph

Bulk averages every cell in a tissue: strong signal per gene, but mixed-population signals average away cell-type-specific changes, and “which cell type” becomes a deconvolution question. Single-cell resolves this at the price of sparser counts and different modeling (as covered in the single-cell explainer on this blog). Choose bulk when the question is “does the sampled tissue change overall”; choose single-cell when subpopulation behavior is the question.

Common failure modes

• Feeding TPM or normalized values into DESeq2/edgeR. Count models need integers; conversion through TPM deletes the count-level variance structure the model exists to model.
• Forgetting to check for batch in the PCA. Everything after a batch-confounded design is decoration.
• Confusing biological and technical replication. Three lanes of one sample are one biological replicate.
• Treating DESeq2/edgeR/limma rankings as interchangeable rankings. They usually share most findings, but borderline genes can differ; report the tool and version rather than asserting a single truth.
• Interpreting small fold changes as false just because p-value barely passed, or large fold changes as real because the p-value looks good at low counts. Shrunken fold changes plus plotted counts give a far safer picture.
• Overwriting raw counts with intermediate normalizations during the analysis — keep the count matrix, the model input, and each intermediate distinct and documented.

Key takeaways

• The pipeline is durable: QC, quantify (alignment+counts or pseudoalignment), design checks, normalize/model, adjust, interpret.
• The data quality question — count provenance and experimental design — outranks the tool choice; DESeq2, edgeR, and limma-voom are three well-regarded legitimate framings, not a ranked podium.
• Variance modeling is why simple tests fail on RNA-seq counts; use a count model, feed it integers, and choose a design formula mapped to the experiment.
• Fold-change shrinkage stabilizes noisy effect estimates for ranking; it is not a normalization step and does not refilter significance.
• Every claim should survive the checks of adjusted p-values, plotted counts, and a framed design statement.

Challenge

Pick any public bulk RNA-seq dataset with two conditions and ≥3 biological replicates per condition (e.g. a GEO series with raw FASTQ or a featureCounts matrix). Rebuild the analysis from Step 3 onward: PCA first, then model with one of DESeq2/edgeR/limma-voom, and compare your significant gene count and top-20 hit list against the paper’s figure. A single completed rebuild from counts is worth more than ten tutorials.

Need help planning or executing an RNA-seq analysis for a project or course? Ampersand Academy offers one-to-one mentoring in bioinformatics and analysis — part of Mahadhi, a development and digital marketing company.

—

FAQ

Bulk measures each gene’s total expression averaged across all cells in a sample and usually focuses on differential expression between conditions. Single-cell measures expression in individual cells and then focuses on cell populations, clustering, and annotation. Costs, count depth, and statistical models differ correspondingly.

• What is the difference between bulk RNA-seq and single-cell RNA-seq analysis?

Both are legitimate. Alignment plus counting is transparent and ties counts directly to gene models; pseudoquantification is fast and produces transcript-level outputs. Pick one, record the annotation version, and do not mix count matrices built from different routes.

• Should I align to the genome or use Salmon/kallisto-style quantification?

No. Both expect raw integer counts, because their variance models are built on count noise structure. Use TPM only for cross-gene visual comparisons, not for count-model input.

• Can I use TPM or FPKM for DESeq2 or edgeR?

It adjusts per-sample library size and composition effects so counts are comparable, typically through size factors (DESeq2) or TMM (edgeR). It deliberately avoids assuming most genes change.

• What does “normalization” actually do in RNA-seq?

Shrinkage stabilizes the effect-size estimates without refiltering significance, which mainly affects ranking and visualization, especially for low-count genes. Significance still comes from the model’s adjusted p-values.

• Why doesn’t fold-change shrinkage change my gene list?

At minimum three biological replicates per condition for reasonable FDR behavior; more is needed for complex designs or noisy tissue. Technical replicates (re-running the same sample) do not substitute.

• How many replicates do I need?

They are all established options; choice is shaped by design complexity, speed needs, and lab precedent rather than a universal winner. The decisive factors are count provenance, design formula, and documented assumptions, not brand.

• Which DE tool should I pick — DESeq2, edgeR, or limma-voom?