RNA-seq and Transcriptomics

The transcriptome is the complete set of RNA molecules present in a cell, tissue, or organism at a particular time — messenger RNAs, ribosomal RNAs, transfer RNAs, and a large population of non-coding and regulatory RNAs (microRNAs, long non-coding RNAs, small nucleolar RNAs). Unlike the genome, which is essentially fixed, the transcriptome changes continuously with cell type, developmental stage, environment, and disease state. RNA-seq (RNA sequencing) is the method of using high-throughput sequencing to capture and quantify this population directly, rather than inferring it indirectly. Two numbers define what RNA-seq measures for each gene: its presence (is it transcribed at all?) and its abundance (how many transcripts per cell, relative to other genes?). Because sequencing reads are counted, not measured by fluorescence intensity as in microarrays, RNA-seq has a much larger dynamic range — it can detect a transcript present at one copy per cell and one present at ten thousand copies per cell in the same experiment, something array hybridization saturates against. RNA-seq also requires no prior knowledge of what sequences to expect, unlike microarrays, whose probes must be designed in advance; this makes RNA-seq equally suited to well-annotated genomes and to novel transcript discovery, allele-specific expression, fusion transcript detection, and alternative splicing analysis. The output of an RNA-seq experiment is, at its simplest, a matrix: rows are genes, columns are samples, and entries are normalized expression values. Differential expression analysis then asks which rows differ systematically between groups of columns — which genes go up or down between a mutant and wild type, a drug-treated and control sample, or a tumor and matched normal tissue. This matrix and the statistics built on it are the entry point to nearly every downstream functional-genomics analysis, from gene ontology enrichment to network inference to the genome-wide regulatory maps built by projects like ENCODE.

Step 1 — RNA extraction and quality control. Total RNA is extracted from cells or tissue. Its integrity is checked (commonly with a bioanalyzer, yielding an RNA Integrity Number or RIN from 1–10) because degraded RNA biases the final measurement toward the 3′ ends of transcripts.

Step 2 — Enrichment or depletion. Ribosomal RNA (rRNA) makes up roughly 80–90% of total cellular RNA and carries little information about differential gene regulation, so it is usually removed by one of two strategies:

• PolyA selection — oligo-dT beads capture the poly-adenylated tails of mature mRNAs, enriching for protein-coding transcripts but excluding most non-coding RNAs and unprocessed transcripts.
• rRNA depletion — probes hybridize to and remove rRNA specifically, retaining mRNAs, long non-coding RNAs, and immature transcripts; this is essential for degraded samples (e.g. formalin-fixed tissue) or for organisms without strong polyadenylation (bacteria).

Step 3 — Fragmentation and reverse transcription. RNA is chemically or enzymatically fragmented into pieces of roughly 200–500 nucleotides (sequencers cannot read full-length transcripts directly at high throughput) and reverse-transcribed into complementary DNA (cDNA), since sequencing platforms read DNA. Strand-specific protocols mark which original strand each fragment came from, which is essential for correctly assigning reads at loci with overlapping genes on opposite strands.

Step 4 — Adapter ligation and amplification. Sequencing adapters are ligated to the cDNA fragment ends, and the library is PCR-amplified to produce enough material for sequencing. Excessive PCR cycles introduce duplicate reads that must later be flagged, since they reflect amplification bias, not true biological abundance.

Step 5 — Sequencing. The library is loaded onto a sequencer (commonly Illumina), generating millions to hundreds of millions of short reads (typically 50–150 base pairs), either single-end (one read per fragment) or paired-end (reads from both ends of each fragment, which better resolves splice junctions and repetitive regions).

Step 6 — Mapping (alignment). Each read must be assigned back to its genomic origin. This step is genuinely harder than DNA resequencing because mature mRNA has had its introns spliced out — a read may span an exon–exon junction that does not exist as a contiguous stretch in the genome. Two computational strategies exist:

• Splice-aware genome alignment (e.g. STAR, HISAT2) maps reads directly to the genome, allowing reads to be split across known or predicted splice junctions.
• Pseudo-alignment / quasi-mapping to a transcriptome (e.g. Salmon, kallisto) skips genomic coordinates entirely and rapidly assigns reads to transcript sequences based on k-mer compatibility, which is far faster and sufficient when only quantification (not novel splice discovery) is needed.

Step 7 — Quantification. For each gene or transcript, the number of reads (or fragments, for paired-end data) mapping to it is counted. This raw count must be corrected for two confounds before genes or samples can be compared: 1. Gene length — a long gene accumulates more reads than a short gene expressed at the same molar concentration, simply because it offers more positions for fragments to land. 2. Sequencing depth — a sample sequenced to 40 million reads will show higher raw counts for every gene than the same sample sequenced to 20 million reads, with no biological difference. This is corrected using FPKM (Fragments Per Kilobase of transcript per Million mapped reads, for paired-end data; RPKM is the single-end equivalent) or, preferably, TPM (Transcripts Per Million):

The critical distinction: TPM normalizes for gene length first, then for sequencing depth, so that TPM values across all genes in a sample always sum to exactly one million — this makes the proportion each gene represents directly comparable between samples. FPKM normalizes depth first, so its per-sample sums vary, making cross-sample comparison technically invalid even though it is common practice. This subtlety is one of the most frequently misapplied points in transcriptomics.

Step 8 — Differential expression statistics. Comparing TPM or FPKM values directly between two conditions using a simple fold-change threshold is statistically naive, because RNA-seq counts are subject to both technical (sequencing) noise and biological variability between replicates, and this noise is not uniform — low-count genes are proportionally noisier than high-count genes. Modern tools (DESeq2, edgeR) instead work from raw counts and model them with a negative binomial distribution, which captures overdispersion (variance greater than the mean) better than the Poisson distribution a purely random sampling process would predict. Each gene's variance is estimated by "borrowing strength" across all genes in the dataset — shrinking noisy per-gene variance estimates toward a shared trend — which stabilizes statistics when only a few replicates (often 3) are available. The output for each gene is a fold-change, a p-value, and a false discovery rate (FDR)-adjusted p-value (correcting for the tens of thousands of simultaneous tests performed across the genome).

RNA-seq did not arrive labeled "accurate" — it had to be validated against pre-existing, trusted methods, and its properties (dynamic range, reproducibility, detection of novel transcripts) had to be demonstrated experimentally.

Mortazavi et al. (2008), Nature Methods — mouse tissue transcriptome mapping. This paper, from the lab of Barbara Wold, is one of the founding demonstrations of RNA-seq quantification and introduced the RPKM normalization scheme still in wide use. The authors sequenced polyadenylated RNA from mouse brain, liver, and skeletal muscle using the Illumina platform, generating tens of millions of short reads per tissue, and mapped them to the mouse genome. They showed that RPKM values were highly reproducible between technical replicates (correlation coefficients above 0.95) and, critically, that RNA-seq quantification correlated strongly with independent measurements from quantitative real-time PCR (qPCR) — the pre-existing "gold standard" for measuring individual transcript abundance. This cross-validation was essential: it proved RNA-seq wasn't just precise (reproducible with itself) but also accurate (agreeing with an independently validated method). The study also directly detected tissue-specific alternative splicing and previously unannotated exons, demonstrating a capability microarrays, limited to pre-designed probes, structurally could not offer.

Marioni et al. (2008), Genome Research — RNA-seq versus microarrays, and technical reproducibility. This study directly compared RNA-seq to microarray platforms on the same liver and kidney RNA samples and found strong concordance in detecting differentially expressed genes, while showing that RNA-seq had a substantially wider dynamic range — it did not saturate at high expression levels the way array fluorescence does, and it could detect low-abundance transcripts that fell below microarray background noise. The paper also directly tested technical reproducibility by sequencing the same RNA sample on multiple lanes of the same flow cell (a technical, not biological, replicate), showing that read counts followed a distribution close to what a Poisson sampling process predicts — an important early result that shaped how the field would later model (and depart from, with the negative binomial) the statistics of count-based expression data.

The MAQC/SEQC consortium studies (2006, 2014). The MicroArray/Sequencing Quality Control consortium, organized in part by the U.S. FDA, ran large multi-site benchmarking studies using reference RNA samples (including the Universal Human Reference RNA) sequenced or hybridized independently at many laboratories on multiple platforms. These studies established the reproducibility of RNA-seq across different sequencing centers and instruments, benchmarked it against TaqMan qPCR measurements on hundreds of genes, and directly informed the field's confidence that RNA-seq-derived fold-changes and expression calls were robust enough for clinical and regulatory use, including in the FDA's own evaluation of biomarker assays. Together these studies established the three claims RNA-seq is now built on: it is reproducible (technical replicates agree), it is accurate (agrees with an independent, orthogonal method — qPCR), and it has superior dynamic range and resolution compared to the hybridization arrays it replaced, including the ability to discover novel transcripts and splice junctions that arrays, by design, could never detect.

Example 1 — Computing TPM by hand. Suppose an RNA-seq experiment on a small three-gene organism yields the following, after paired-end sequencing and mapping:

Step 1 — divide fragment counts by gene length in kb (this gives "reads per kb", RPK):

Step 2 — sum the RPK values across all genes:

Step 3 — divide each RPK by this sum, and multiply by one million:

Check: 227,273 + 681,818 + 90,909 = 1,000,000 (up to rounding) — TPM values for all genes in a sample always sum to one million, which is exactly what makes them comparable across samples with different total read depths. Note the biological point this reveals: gene B has by far the fewest raw fragments (3,000, less than A's 4,000) yet is the most highly expressed gene per molecule once its short length (0.5 kb) is accounted for — raw counts alone would have ranked it lowest. This is precisely the length bias TPM is designed to correct.

Example 2 — Deciding whether a gene is truly differentially expressed. A researcher compares a drug-treated cell line (3 replicates) against untreated controls (3 replicates) and finds:

A naive read of "fold-change" alone would flag Reg-X as the more dramatic and exciting hit, with 8-fold induction versus Hsp70's 2-fold. But raw count magnitude matters enormously for statistical confidence. Hsp70 is measured from thousands of reads per replicate; a doubling from ~4,000 to ~8,000 reads is a large, well-sampled signal, and the negative binomial model built into DESeq2 or edgeR will assign it a very small, confident p-value once the modest between-replicate variance is accounted for. Reg-X, in contrast, is measured from single-digit read counts — the kind of gene where one replicate happening to get 3 reads instead of 9 by chance alone can swing the fold-change by several-fold, since Poisson-like sampling noise at low counts is proportionally huge (a change of even 2–3 reads is a large fraction of the total). A count-based statistical test will typically assign Reg-X a large, non-significant p-value despite its impressive-looking fold-change, correctly flagging it as unreliable rather than biologically important. This is exactly why differential expression tools work from raw counts and model count-dependent variance directly, rather than filtering on fold-change computed from already-normalized TPM values — a mistake that would systematically favor exactly this kind of false positive among lowly expressed genes.

RNA-seq is a technology, not a single scientist's discovery, but Barbara Wold's laboratory at Caltech produced one of the two founding papers (2008) that established RNA-seq as a quantitative, genome-wide expression method and introduced RPKM, the normalization scheme from which TPM was later derived. Wold's background was in developmental biology and gene regulation — she had spent much of her earlier career studying muscle differentiation and transcription factor networks using the pre-genomic tools of Northern blots and in situ hybridization, which made her acutely aware of how limited single-gene assays were for understanding genome-wide regulatory programs. When massively parallel sequencing became available, her lab was among the first to recognize it could be repurposed not just for genome sequencing but for direct, unbiased transcript counting, and to rigorously benchmark the new method against the trusted quantitative PCR standard rather than simply asserting its superiority.

Cancer genomics and diagnostics. Tumor transcriptome sequencing identifies fusion transcripts (such as BCR-ABL1 in chronic myeloid leukemia), classifies tumors into molecular subtypes invisible to histology alone, and reveals which signaling pathways are aberrantly active, guiding targeted therapy selection.

Drug development and mechanism-of-action studies. Pharmaceutical screens routinely use RNA-seq to characterize how a candidate compound reprograms cellular gene expression, revealing on-target effects, off-target toxicity signatures, and resistance mechanisms long before a drug reaches clinical trials.

Rare disease diagnosis. When genome or exome sequencing finds a variant of uncertain significance, RNA-seq from the patient's own tissue can show directly whether that variant disrupts splicing or silences the gene, often solving diagnostic cases that DNA sequencing alone left ambiguous.

Agriculture and crop improvement. Comparing the transcriptomes of drought-tolerant and drought-sensitive crop varieties under stress identifies candidate genes for breeding programs aimed at climate resilience, without needing to know in advance which genes to look at.

Conservation and ecological genomics. RNA-seq applied to non-model, non-reference-genome species (via de novo transcriptome assembly) allows researchers to study stress responses, immune function, or population differences in endangered or ecologically important organisms that lack a sequenced genome.

Single-cell and spatial transcriptomics. The core RNA-seq logic — fragment, sequence, map, count, normalize — now extends to single-cell resolution (revealing cellular heterogeneity within a "bulk" tissue that averages it away) and to spatial transcriptomics, which preserves the physical location of each transcript within a tissue section, both direct technical descendants of the bulk RNA-seq pipeline described here.

Ask the tutor: How would you redesign an RNA-seq experiment to reliably detect a rare, lowly-expressed transcription factor that shows only a modest 1.5-fold induction?