Bulk RNA-seq: from biological sample to interpretable signal

This guide is intentionally split into two parts:

• Upstream: everything from biological sample and library design through QC, quantification, and the construction of a trustworthy count table.
• Downstream: everything after the count table, starting with core differential-expression analysis and then moving into the extra layers that help you squeeze more biology out of the dataset.

If you want one default path that works for most human or mouse bulk RNA-seq projects, make it this:

1. Plan the biology first: replication, RNA quality, library type, strandedness, and annotation version.
2. Treat FastQC + MultiQC as the universal first stop after FASTQ generation.
3. For expression-first work, quantify with Salmon and import with tximport, ideally keeping provenance with tximeta.
4. Use DESeq2 as the default differential-expression engine.
5. Only after you trust the counts and the design should you add enrichment, pathway activity, TF activity, cell-context inference, and more advanced transcript-aware analyses.

This page is not trying to list every transcriptomics method ever published. It is a practical guideline to the core analysis up to DEGs, and then to the downstream analyses that are most useful when you want to get more biological signal out of bulk RNA-seq without drifting into method tourism. The mainline recommendation is:

Clean mammalian tissue, expression-first question, reasonable RNA quality

poly(A) enrichment -> stranded library -> FASTQ QC -> Salmon or nf-core/rnaseq -> tximport/tximeta -> DESeq2 -> enrichment + pathway activity + TF activity + cell-context inference

Switch away from that default when the question changes:

• Choose rRNA depletion when RNA is degraded, FFPE-like, intronic signal matters, or non-polyadenylated transcripts matter.
• Choose STAR when you need genome alignments for genome-browser coverage, exon-aware counting, or genome-position-aware follow-up.
• Choose paired-end and longer reads when genome-aware follow-up is an important part of the study rather than an afterthought.

The upstream half of the workflow is about getting from biological material to a count table you believe. If this part is weak, no downstream method will save it.

Before sequencing

The best bulk RNA-seq analysis begins before the sequencer sees a library. Bad experimental structure is harder to rescue than bad code.

Planning checklist

Library design choices

Annotation and reference strategy

Use one reference bundle and freeze it:

• GENCODE: strong default for human and mouse when you want well-curated gene and transcript models.
• Ensembl: strong default when your ecosystem or annotation tooling is already Ensembl-centric.
• NCBI RefSeq: use when a collaborator, institution, or downstream clinical context is already standardized on RefSeq.

Do not mix:

• GENCODE transcriptome index with RefSeq gene annotation in downstream joins
• Ensembl IDs from one release with symbols from another release
• human-centric signature resources onto mouse data without explicitly documenting ortholog conversion

Replication, power, and controls

Biological replication matters more than squeezing one more downstream algorithm out of underpowered data.

• Use RNASeqPower or powsimR when you can still influence the experiment.
• If you are already underpowered, be honest about effect-size uncertainty and treat pathway-level summaries as supportive, not rescuing.
• ERCC-style spike-ins can be useful in selected designs, but they do not replace good biological replication and they can create a false sense of security if used as a generic fix.

# Simple planning logic, not a substitute for study-specific assumptions
library(RNASeqPower)
est_power <- rnapower(
  depth = 20,
  n = 4,
  cv = 0.3,
  effect = 1.5,
  alpha = 0.05
)
est_power

From FASTQ to trusted quantification

Keep the data type explicit at each step:

• FASTQ: raw reads and base qualities
• BAM/CRAM: genome alignments
• quant.sf or equivalent: transcript-level abundance estimates
• gene count matrix: integer counts for count-based DE frameworks
• normalized matrix: transformed expression for visualization and sample structure
• ranked vector: gene-wise statistics for GSEA
• splice-event or exon tables: event-aware follow-up
• regional scores or genome-aware summaries: specialized outputs

Universal QC comes first

mkdir -p qc/fastqc
fastqc data/fastq/*.fastq.gz --outdir qc/fastqc
multiqc qc -o qc/multiqc

For RNA-specific QC after a genome-alignment step, check:

• RSeQC: infer_experiment.py, read_distribution.py, geneBody_coverage.py
• RNA-SeQC: broad process metrics, coverage and expression summaries
• Picard CollectRnaSeqMetrics: coding, UTR, intronic, intergenic fractions, ribosomal bases

infer_experiment.py -r refs/hg38.bed -i align/sample.bam
geneBody_coverage.py -r refs/hg38_housekeeping.bed -i align/sample.bam -o qc/rseqc/sample
picard CollectRnaSeqMetrics \
  I=align/sample.bam \
  O=qc/picard/sample.rna_metrics.txt \
  REF_FLAT=refs/gencode.vM35.refFlat.txt \
  STRAND_SPECIFICITY=SECOND_READ_TRANSCRIPTION_STRAND

Trimming: useful, but do not do it reflexively

trim_galore --paired \
  --cores 4 \
  sample_R1.fastq.gz sample_R2.fastq.gz \
  --output_dir trimmed/

Quantification strategy

salmon quant \
  -i refs/gencode.v44.salmon_index \
  -l A \
  -1 trimmed/sample_R1.fq.gz \
  -2 trimmed/sample_R2.fq.gz \
  -p 8 \
  --validateMappings \
  -o quants/sample

library(tximport)
library(tximeta)

coldata <- read.csv("metadata/samples.csv")
se <- tximeta(coldata)
gse <- summarizeToGene(se)

The downstream half starts once you have a count table and credible metadata. The first job is still boring and essential: import, QC, design, normalization, and differential expression. Only then does it make sense to start asking more ambitious biological questions.

Core analysis

Once you have a count-like object and trustworthy metadata, the work becomes statistical and interpretive rather than purely computational.

Import, normalization, and sample structure

For most users, the best default is DESeq2 on imported gene-level counts. If you quantified with Salmon, do not throw away transcript-level work by flattening it carelessly. Use tximport, and use tximeta when you can.

library(DESeq2)
library(tximeta)

se <- tximeta(read.csv("metadata/samples.csv"))
gse <- summarizeToGene(se)

dds <- DESeqDataSet(gse, design = ~ batch + condition)
dds <- DESeq(dds)

vsd <- vst(dds, blind = FALSE)
plotPCA(vsd, intgroup = c("condition", "batch"))

Differential expression

res <- results(dds, contrast = c("condition", "treated", "control"))
res_shrunk <- lfcShrink(dds, coef = "condition_treated_vs_control", type = "apeglm")

sig <- subset(as.data.frame(res_shrunk), padj < 0.05 & abs(log2FoldChange) > 1)
head(sig[order(sig$padj), ])

Batch effects and unwanted variation

library(sva)

mod  <- model.matrix(~ condition, data = colData(dds))
mod0 <- model.matrix(~ 1, data = colData(dds))
sv   <- svaseq(counts(dds), mod = mod, mod0 = mod0)$sv

colData(dds)$sv1 <- sv[, 1]
design(dds) <- ~ sv1 + batch + condition
dds <- DESeq(dds)

What to inspect before believing any DE table

• PCA or MDS on transformed counts
• Sample-sample distances
• Library-size spread and composition bias
• Gene-body coverage and strandedness sanity checks
• Outlier samples in MultiQC and RNA-specific QC
• Whether low-count filtering was reasonable
• Whether the design formula matches the actual experiment

Slow down before interpretation

Downstream layers after DE

Once the DE result is credible, move upward in interpretation one layer at a time. Start with individual genes, then ask whether coherent gene sets shift, then whether individual samples carry programs, then whether those programs imply signaling or regulatory activity, and finally whether the bulk profile may reflect changing tissue composition.

Gene-level

This is the anchor layer. Differential expression tells you which genes move and by how much. Everything that follows should be read as an interpretation layered on top of that gene-level result, not as a substitute for it.

Stay here until you trust:

• the contrast and design formula
• the effect sizes and shrinkage behavior
• the sample structure and batch handling
• a few manually inspected genes that fit the biology you think you are seeing

Once you have a ranked gene-level result you believe, the next shift is from individual genes to coordinated programs.

Gene set-level

This is the first abstraction step. Instead of asking whether one gene is up or down, ask whether a pathway or process moves as a coherent set.

library(fgsea)

stats <- res_shrunk$stat
names(stats) <- rownames(res_shrunk)
stats <- sort(stats, decreasing = TRUE)

hallmark <- gmtPathways("refs/h.all.v2024.1.Hs.symbols.gmt")
fg <- fgsea(pathways = hallmark, stats = stats, minSize = 15, maxSize = 500)
head(fg[order(fg$padj), c("pathway", "NES", "padj")])

library(clusterProfiler)

ora <- enrichGO(
  gene = rownames(subset(sig, log2FoldChange > 1)),
  OrgDb = org.Hs.eg.db,
  keyType = "ENSEMBL",
  ont = "BP"
)

ORA and GSEA still summarize the group-level contrast. If the next question is whether individual samples carry a program strongly or weakly, move from gene sets to per-sample scores.

Sample-level pathway activity

Here the unit of interpretation shifts from the contrast to the sample. GSVA and ssGSEA turn a normalized expression matrix into per-sample pathway scores, which is useful when heterogeneity, subtype structure, responder patterns, or outlier biology matter.

Use this same layer for cell state or functional state scoring. In practice that means scoring signatures such as:

• cell cycle
• proliferation
• hypoxia
• EMT
• interferon response
• stress or injury programs

These are still sample-level summaries. They tell you which samples carry a program, not yet which upstream regulator or pathway driver caused it.

library(GSVA)

ssgsea_scores <- gsva(
  as.matrix(vst_mat),
  hallmark,
  method = "ssgsea"
)

Once you have sample-level programs, the next shift is from descriptive scoring to inferred upstream drivers.

Regulatory / signaling inference

This layer asks what might be driving the expression pattern rather than just describing it. PROGENy estimates signaling pathway activity from downstream responsive genes, while DoRothEA + VIPER estimates transcription factor activity from coordinated target-gene behavior.

library(progeny)

pathway_scores <- progeny(vst_mat, scale = TRUE, organism = "Human", top = 500)

library(dorothea)
library(viper)

data(dorothea_hs, package = "dorothea")
regulon <- subset(dorothea_hs, confidence %in% c("A", "B", "C"))
tf_scores <- viper(vst_mat, regulon, method = "scale")

These methods still assume the bulk profile is mostly reflecting state within the sampled mixture. If the mixture itself is changing, apparent pathways or TFs may partly be composition effects rather than cell-intrinsic regulation.

Tissue composition

This is the mixture layer. Here the question shifts from what programs are active to who is present in the sample. In bulk RNA-seq, that distinction matters because apparent inflammation, stromal biology, or differentiation can reflect changing proportions rather than changed state within one cell type.

Recommended practical default:

• start with xCell2 for a conservative context readout
• use EPIC or CIBERSORTx when the reference assumptions fit the biological system
• use MuSiC when you have a truly matched single-cell reference and want proportion-aware estimation

At this point you should know whether your story is mainly gene-level, program-level, sample-state-level, regulatory, or compositional. If the remaining question is specifically about isoforms, splicing, or genomic position, move into a more structure-aware layer.

Genome-aware / transcript-level analysis

This is the advanced or optional layer. Once you collapse signal to gene-level counts, you lose exon structure, isoform usage, and genomic position. These methods are for cases where that lost structure is the actual biological question.

Use this layer when the unresolved question is about transcript structure or genomic organization, not when the real question is pathway membership, sample state, or tissue composition.

library(DEXSeq)
library(DRIMSeq)

rmats.py --b1 treated_bams.txt --b2 control_bams.txt \
  --gtf refs/genes.gtf \
  --od rmats_out --tmp rmats_tmp \
  --readLength 101 --nthread 8

library(PREDA)

Targeted panels and ortholog sanity checks

Do not let an enrichment package decide your biology for you. Always return to a few hand-curated gene panels:

• interferon and inflammatory genes
• ECM and adhesion genes
• mitochondrial and stress genes
• lineage markers relevant to your tissue
• known perturbation targets

For mouse studies pushed through human-centric resources, create an explicit ortholog table and keep both original and mapped identifiers side by side.

Common pitfalls

This section is intentionally blunt. Most bulk RNA-seq pain comes from a small set of repeated mistakes.

Technical and analytical pitfalls

1. No real replication. Statistical sophistication does not replace biological replication.
2. Annotation drift. If the index, tx2gene mapping, and visualization symbols came from different annotation releases, joins will break or silently lie.
3. Wrong strandedness assumption. Use infer_experiment.py or equivalent checks. Counting with the wrong strand setting can quietly invert confidence in some loci.
4. Blind trimming. Trim because there is a problem, not because it feels hygienic.
5. Ignoring multi-mapping. Repetitive genes, paralogs, pseudogenes, and short-read ambiguity are biological and computational realities.
6. Over-correcting batch. Batch correction can remove signal of interest if design and biology are entangled.
7. Short reads for coordinate-aware questions. If the goal is genome-aware interpretation or positional follow-up, cheap short single-end data are often the real limitation.
8. rRNA contamination denial. High rRNA fractions are not just annoying; they waste depth and can indicate library problems.
9. Assuming weak signal means no biology. It can also mean wrong tissue compartment, weak perturbation, bad timing, or overwhelming composition changes.
10. Treating every downstream method as mandatory. Most projects do not need every specialized scoring, regional, and network layer all at once.

What to do when signal is weak or libraries look bad

• Re-check sample labels and metadata first.
• Re-check strandedness and annotation consistency.
• Inspect gene-body coverage and read distribution.
• Ask whether the perturbation timing or tissue compartment was wrong.
• Look at known positive-control genes or pathways before declaring failure.
• Consider that composition shifts may dominate weak within-cell-state changes.
• If everything depends on rare populations or spatial structure, bulk may not be the right assay.

Interpretation pitfalls

This final list focuses on interpretation traps that still matter even after you have already checked for technical structure and composition effects.

1. Compensatory versus functional change. Upregulation does not imply more function. It may reflect compensation after stress, damage, or failed homeostasis.
2. Structure versus expression mismatch. Higher mRNA does not guarantee correct protein structure, localization, or tissue architecture.
3. Stress and injury responses are nonspecific. Hypoxia, inflammation, interferon, EMT-like programs, and immediate-early responses are reused across many contexts.
4. Pathway redundancy. Overlapping gene sets can make several pathways look independently significant when they really reflect one underlying program.
5. Directionality ambiguity. Up or down does not automatically define activation, inhibition, or biological outcome.
6. mRNA is not protein and not function. Translation, post-translational regulation, localization, and turnover matter.
7. Temporal ambiguity. A single time point cannot cleanly separate cause from consequence.
8. Reverse causality. Especially in tumor and inflammatory settings, observed changes may be downstream responses rather than drivers.
9. Signatures are not mechanisms. Enrichment indicates similarity to a known program, not mechanistic activation.
10. Cross-species limitations. Mouse-to-human ortholog mapping loses nuance and sometimes meaning.

If you want the shortest honest recommendation:

1. Design the experiment around the biological question, not around whichever downstream package you hope to use.
2. For most human and mouse bulk RNA-seq projects, use stranded poly(A) if RNA quality is good and the question is coding-gene expression.
3. Start every dataset with FastQC + MultiQC, then add RSeQC, RNA-SeQC, or Picard when RNA-specific sanity checks matter.
4. Use Salmon + tximport + tximeta for expression-first work.
5. Use STAR when you genuinely need genome-aware analyses.
6. Use DESeq2 for gene-level DE, then add lfc shrinkage, the right flavor of gene-set analysis, one sensible regulatory inference layer, and a conservative cell-context layer.
7. Treat genome-aware regional analysis as a targeted follow-up when the biological question justifies it, not as a default decoration.
8. Pick downstream methods based on the biological question they answer, not because they are popular on other people's RNA-seq figures.