Overview
Curated: · Written: · Reviewed:
RNA-seq measures molecules through an experimental and statistical model
RNA-seq turns RNA molecules into library fragments, sequenced reads and abundance estimates, and every stage of that conversion changes what you can conclude. Tissue composition, collection time, RNA integrity, extraction batch, selection or depletion, strandedness, fragment length, amplification and lane effects can each create systematic differences larger than the biological condition you are studying. The interview-ready framing is this: RNA-seq is a measurement instrument, and the question "what did you actually measure?" is answerable only by walking the pipeline from sample to p-value.
The pipeline and its decision points
Everything starts with the RNA itself. Extraction yield and integrity (RIN, the RNA integrity number, scored 1–10) determine what is even possible: a RIN below roughly 5 on degraded or FFPE material means fragmented RNA, which rules out poly-A selection — the 3′ bias of selection on degraded RNA would measure fragment ends, not transcripts — and pushes you toward rRNA depletion or targeted approaches.
Library strategy is the first big fork:
- Poly-A selection enriches mRNA, gives clean coding-gene coverage at lower depth, and makes non-coding RNA, lncRNAs and degraded fragments invisible.
- rRNA depletion keeps non-polyadenylated species — lncRNAs, intronic and antisense transcription, bacterial RNA in host-microbiome work — at the cost of substantial rRNA contamination you must QC for.
- Total RNA with depletion variants is the fallback for degraded samples and non-coding questions.
Strandedness decides whether a read overlapping two antisense features can be assigned. Unstranded libraries lose that information; stranded libraries (record which protocol — dUTP vs. ligations have opposite read conventions) let you quantify antisense transcription. Getting the strandedness flag wrong in counting silently halves or zeroes your counts for correctly stranded genes.
Read length, depth and pairing trade against replicate number. Paired-end reads improve isoform resolution, junction detection and quantification of overlapping transcripts; single-end is cheaper. Depth beyond roughly 10–30 million reads for gene-level differential expression in bulk tissue buys little — the standard power analysis (RNASeqPower, or the Liu et al. curves) says replicates beat depth once genes are adequately observed. Isoform-level work, fusion detection and rare-cell-type questions need more depth; a gene-level comparison does not.
Quantification: align-and-count vs. pseudo-alignment
Two families of pipelines, with a real trade-off:
- Splice-aware alignment (STAR, HISAT2) maps reads to the genome across introns, giving genomic coordinates, junction evidence and the ability to detect novel splicing and fusions. featureCounts or HTSeq then assigns fragments to genes under explicit rules: feature type (exon), grouping attribute (gene_id), strandedness, whether a pair counts as one fragment, and how multimappers and overlapping features are handled.
- Pseudo-alignment / selective alignment (kallisto, salmon) quantifies against a transcriptome index, distributing ambiguous fragments probabilistically across compatible isoforms. It is much faster and gives transcript-level abundances with bootstrap uncertainty, but it cannot discover unannotated transcription, and its accuracy is bounded by annotation completeness.
The core ambiguity is the same in both: isoforms share exons, so many fragments are compatible with multiple transcripts. Transcript-level estimates are therefore correlated and individually uncertain — a fact interviewers probe with "why are your transcript-level p-values so much weaker than gene-level ones?" Gene-level aggregation is the robust default; go to transcript level only when the biology is isoform-specific and you carry the inferential uncertainty through.
Reference genome and annotation release form one coordinate system. Mixing a GTF from one Ensembl release with a genome build from another silently loses or misassigns reads. Keep stable IDs (ENSG…) with version suffixes for joins, and attach gene symbols as attributes — symbols are renamed, duplicated and retired across releases and are neither stable nor unique.
Normalization: why raw counts go into the model
A worked example, with hypothetical figures, showing a conclusion that reverses:
| Sample A | Sample B | |
|---|---|---|
| Raw counts for GENE1 | 800 | 1,200 |
| Library size (total mapped reads) | 20,000,000 | 60,000,000 |
| CPM | 40.0 | 20.0 |
| Gene length | 2,000 bp | 2,000 bp |
| TPM | 21.7 | 10.8 |
Raw counts say GENE1 is 1.5× higher in B. Normalized for library size, it is 2× higher in A. The raw comparison was measuring sequencing depth, not expression.
Length matters for a different comparison. Within one sample, a 10,000 bp transcript accumulates roughly 5× the reads of a 2,000 bp transcript at equal molar abundance, so CPM ranks long genes above short ones:
Same gene ACROSS samples: library size matters, length cancels -> CPM is enough
Different genes WITHIN one sample: length matters -> TPM or FPKM
Differential expression: neither -- DESeq2/edgeR take raw integer
counts and fit size factors internally
The last line is the one candidates most often get wrong. Feeding TPM into DESeq2 breaks it: the negative binomial likelihood needs integer counts to estimate the mean–variance relationship, and pre-normalized continuous values destroy exactly the information the model uses.
Library-size normalization is also not simple division by totals. If one highly expressed transcript takes 30% of sample B's library, every other gene in B is depressed by composition alone, and CPM will not fix it. DESeq2's median-of-ratios computes, per gene, the ratio of each sample's count to the geometric mean across samples, then takes the median ratio per sample as its size factor — robust because it assumes most genes do not shift in one direction. edgeR's TMM (trimmed mean of M-values) trims the log-fold-changes and absolute intensities with the largest extremes before averaging, under the same majority-stable assumption. Both fail when that assumption fails: a genuine global transcriptional shift (mitochondrial, dosage, global activation) is absorbed into the size factor and erased. Spike-ins can support absolute or global-shift questions only when their addition and recovery are technically controlled.
TPM is the right unit for reporting within-sample composition and for some visualizations; it is not the input to a count-based DE test.
The statistical model: negative binomial, dispersion, design
RNA-seq counts are overdispersed relative to Poisson — biological variability between replicates makes the variance grow faster than the mean. The negative binomial adds a dispersion parameter α so that Var = μ + αμ². With n = 3 replicates per group, a per-gene variance estimate from the data alone is nearly meaningless, so DESeq2 and edgeR share information across genes: DESeq2 fits a gene-wise estimate, a fitted trend against mean expression, and shrinks each gene toward the trend (empirical Bayes moderation). edgeR does the same with a common/trended/tagwise dispersion hierarchy. This is why these tools work at all at small n — and why n below 3 gives dispersion estimates that are mostly borrowed assumption, with power to detect only large effects.
The GLM design formula connects samples to coefficients: ~ batch + condition models log-scale expression as an intercept, a batch coefficient and a condition coefficient. The contrast is the specific coefficient combination that answers your question — condition_B_vs_A in the simplest case, or an interaction term in a time-course or paired-plus-treatment design. Write the estimand in words before extracting it; a reversed reference level flips every sign silently.
Fold-change shrinkage (DESeq2's lfcShrink, apeglm) is separate from dispersion shrinkage. Low-count genes produce unstable extreme ratios — 4 reads vs. 1 read is a 4× fold change with no support — so shrinkage pulls those toward zero, trading a little bias for much better ranking and visualization. Use shrunken LFCs for ranking and plots; the test statistics for significance come from the unshrunken model.
Report effect estimate, uncertainty and normalized counts, not a thresholded gene list. A tiny effect can be significant in a large study; an important effect can be uncertain in a small one.
Multiple testing: BH, FDR, and what an adjusted p-value means
You tested ~20,000 genes, so raw p-values are meaningless as evidence — at p < 0.05 you expect ~1,000 false positives under the null. Benjamini–Hochberg controls the false discovery rate: sort p-values, find the largest i where p(i) ≤ i·q/N, and everything at or above it is significant at FDR q. The guarantee is about the expected proportion of false discoveries among the set under repeated use — not the probability that any particular gene is false, and not family-wise error (the probability of even one false positive, which Bonferroni controls far more conservatively).
Interpret an FDR cutoff alongside effect size: a significant adjusted p-value with a log2FC of 0.1 is usually noise-with-power, and a log2FC of 3 at FDR 0.2 in a pilot may be the finding worth pursuing. Independent filtering (removing very low-count genes before testing) improves power when the filter statistic is independent of the test under the null — low average expression is such a statistic, which is why DESeq2 does it automatically.
Failure modes, confounding, and what interviewers probe
The canonical failure, as a worked example: a study sequenced 8 treated libraries from plate A and 8 controls from plate B. A ~condition fit reported 1,842 genes at FDR 0.05. Adding batch to the design surfaces the problem — DESeq2 refuses a rank-deficient design at construction, because condition and plate are the same column and no software can separate them. The 1,842 genes were never a treatment result; they were a treatment-or-plate result, and the experiment cannot say which. The lesson generalizes: if condition is perfectly confounded with batch, no downstream correction can identify which caused the difference. Batch correction (including it in the design, or methods like ComBat for visualization) works only when batch and condition are identifiable — each batch contains both conditions.
Other failure modes worth naming:
- Bulk averages cells. A changed gene count can reflect regulation within cells, altered cell composition, or both. Check markers, histology, deconvolution or single-cell evidence before claiming regulation.
- Outcome-driven outlier removal invalidates inference. Inspect outliers against identity, QC and influence; exclude only under predefined rules and report sensitivity with and without.
- PCA on raw counts is dominated by library size. Use a variance-stabilizing transform; label known covariates on the plot; do not regress away a component that may contain the biology.
- Gene-set analysis inherits biases: identifier mapping, universe selection (the universe is the genes that could have been tested, not every gene in a database), gene-length and correlation biases. A pathway label is a hypothesis, not mechanism.
What interviewers probe, and what a weak answer sounds like: expect "why not just use TPM/FPKM for DE?" (weak answer: "TPM is normalized" — strong answer names the integer-count likelihood and composition assumptions); "what does an FDR of 0.05 mean?" (weak: "5% chance this gene is false" — that is the per-gene misreading); "you have 3 vs. 3 replicates and 40 genes at FDR 0.05 — what do you trust?" (probe for power, dispersion borrowing, and independent validation); "how would you handle a global transcriptional shift?" (probe for spike-ins and the failure of median-of-ratios assumptions). Likely follow-ups: paired designs and subject terms, interaction contrasts in time courses, and what you would change at scale — thousands of samples bring batch-aware distributed processing (e.g. SNRNA-seq-style workflows), identity and consent metadata protection (expression data can re-identify individuals), and immutable provenance: checksums, reference and annotation versions, counting rules, design formula, contrasts, session info and every exclusion.
The production invariant is design-valid inference: every reported expression change traceable to exact samples, molecules, reference and counting rules, a full-rank model and a declared contrast, with technical and biological uncertainty, multiplicity and alternative explanations kept visible.