RNA-Seq

RNA-Seq (whole-transcriptome shotgun sequencing) is a high-throughput molecular and bioinformatics methodology that uses massively parallel sequencing or long-read DNA sequencing platforms to catalog and quantify the cellular repertoire of RNA transcripts at a given moment. Introduced in 2008 by Zhong Wang, Mark Gerstein, and Michael Snyder, alongside independent studies by Ali Mortazavi, Barbara Wold et al. and Ugrappa Nagalakshmi et al., RNA-Seq advanced functional transcriptomics by replacing hybridization-based DNA microarrays with direct digital sequence counting.

Unlike hybridization microarrays—which require pre-designed species-specific oligonucleotide probes, suffer from background cross-hybridization noise, and saturate at high fluorescence intensities—RNA-Seq requires no prior genomic sequence annotation, resolves transcript boundaries to single-nucleotide accuracy, detects novel alternative splicing isoforms, allele-specific expression, and oncogenic gene fusions, and quantifies transcript abundance across a >10⁵-fold dynamic range (from <1 transcript copy per cell to >100,000 copies) using 20 to 50 million paired-end reads per sample.

RNA Integrity (RIN), Poly(A)+ Selection, and Ribosomal RNA Depletion

Because cellular RNA is susceptible to degradation by cellular RNases due to the reactive 2′-hydroxyl (2′-OH) group on its ribose sugar ring, RNA-Seq library preparation begins with capillary electrophoretic quality control on an Agilent Bioanalyzer. The RNA Integrity Number (RIN)—scaled from 1 (degraded) to 10 (intact)—evaluates the ratio of 28S to 18S ribosomal RNA (rRNA) peaks (~2.0:1 in intact mammalian RNA); standard poly(A)-selected RNA-Seq requires RIN ≥ 7.0 to 8.0.

In cellular RNA extracts, ribosomal RNAs (28S, 18S, 5.8S, and 5S rRNA) and mitochondrial rRNAs constitute 80% to 90% of cellular RNA mass, whereas protein-coding messenger RNAs (mRNAs) represent only 1% to 5%. To prevent rRNA from consuming 90% of sequencing reads, transcripts are enriched via one of two biochemical pathways:

1. Oligo(dT) Paramagnetic Bead Selection (Poly(A)+ Enrichment): Hybridizes synthetic oligo(dT)₂₅ chains coupled to magnetic beads against the 3′ polyadenylated tail (50–250 adenosine residues) of mature eukaryotic mRNAs and polyadenylated long non-coding RNAs (lncRNAs).
2. Ribosomal RNA Depletion (RNase H or Ribo-Zero): Hybridizes antisense DNA oligonucleotides tiled across rRNA sequences, followed by RNase H enzymatic digestion of the RNA strand in DNA:RNA heteroduplexes. rRNA depletion is required for degraded formalin-fixed paraffin-embedded (FFPE) clinical biopsies (where fragmented 5′ mRNA ends have lost their 3′ poly(A) tails), prokaryotic bacteria (which lack poly(A) tails), histone mRNAs, and circular RNAs (circRNAs).

Comparison of RNA-Seq Library Modalities and Expression Normalization Metrics

The dUTP Second-Strand Marking Method for Stranded RNA-Seq

Early RNA-Seq protocols converted RNA into double-stranded complementary DNA (cDNA) without retaining information on which DNA strand served as the transcriptional template, confounding overlapping sense and antisense genes. Modern strand-specific (stranded) RNA-Seq employs the dUTP second-strand marking method introduced by Dmitry Parkhomchuk and colleagues in 2009:

1. Enriched RNA is chemically fragmented with divalent cations (Mg²⁺ or Zn²⁺) at 94 °C into 200–300 nt fragments and primed with random hexamers.
2. Reverse transcriptase (SuperScript II/IV or M-MLV) synthesizes the first-strand cDNA using canonical dNTPs (dATP, dCTP, dGTP, dTTP), generating an RNA:cDNA hybrid.
3. During second-strand cDNA synthesis (catalyzed by E. coli RNase H, DNA Polymerase I, and E. coli DNA Ligase), deoxyuridine triphosphate (dUTP) is substituted for dTTP, incorporating uracil bases into the second (sense) cDNA strand.
4. After Y-shaped adapter ligation, Uracil-DNA Glycosylase (UDG / USER enzyme) excises uracil bases and degrades the dUTP-containing second strand prior to PCR amplification, preserving >99% strand-of-origin specificity.

Splice-Aware Alignment (STAR, HISAT2) vs. Pseudoalignment (kallisto, Salmon)

Because mature eukaryotic mRNAs have had their non-coding introns excised by the spliceosome, a 150-bp RNA-Seq read spanning an exon-exon junction maps discontinuously across thousands of base pairs of genomic reference DNA. Consequently, standard DNA aligners (BWA, Bowtie2) cannot map junction reads. Instead, splice-aware aligners—such as STAR (Spliced Transcripts Alignment to a Reference, by Alexander Dobin et al., 2013), which uses uncompressed suffix arrays, and HISAT2, which uses a hierarchical graph FM-index—split junction-spanning reads across canonical GT-AG splice sites to assemble transcripts (StringTie) and count gene-level reads (featureCounts, HTSeq).

Alternatively, alignment-free transcript quantifiers—kallisto (Nicolas Bray, Lior Pachter et al., 2016) and Salmon (Rob Patro et al., 2017)—bypass base-by-base genomic alignment. By building a colored transcriptome De Bruijn graph of annotated transcript k-mers, these algorithms perform pseudoalignment (or quasi-mapping) to determine the set of transcript isoforms compatible with each read, using an Expectation-Maximization (EM) likelihood algorithm to assign multimapping isoform counts in minutes.

Why TPM Replaced RPKM/FPKM and Negative Binomial Modeling (DESeq2 / edgeR)

Because a 5,000-bp transcript produces five times as many fragmented sequencing reads as a 1,000-bp transcript expressed at the same molar concentration, raw read counts (C_i) must be normalized by effective transcript length (L_i). Mortazavi et al. (2008) introduced RPKM (Reads Per Kilobase of transcript per Million mapped reads; FPKM for paired-end fragments), which divides by library size first and gene length second. However, because the sum of RPKM/FPKM values across all genes varies from sample to sample whenever the mean expressed transcript length changes, FPKM values cannot be compared directly across samples. TPM (Transcripts Per Million) (proposed by Bo Li, Colin Dewey, and Günter Wagner) normalizes by transcript length C_i / L_i first and then scales by the sum of length-normalized rates ∑ (C_j / L_j) / 10⁶, ensuring that ∑_i TPM_i = 1,000,000 in every sample.

For statistical differential gene expression (DGE) testing between experimental conditions, raw integer count matrices are analyzed in DESeq2 (Michael Love, Wolfgang Huber, and Simon Anders, 2014) or edgeR (Mark Robinson, Davis McCarthy, and Gordon Smyth, 2010). Because biological replicates exhibit overdispersion (variance σ² exceeds the mean μ predicted by a Poisson model), DESeq2 models read counts K_ij for gene i in sample j using a Negative Binomial generalized linear model:

Var(K_ij) = μ_ij + α_i μ_ij²

where size factors normalize library composition (median-of-ratios method), gene-wise dispersion α_i is stabilized via empirical Bayes shrinkage, and Wald or likelihood-ratio p-values are adjusted via the Benjamini-Hochberg False Discovery Rate (FDR).

Why can RNA-Seq raw read counts not be modeled with a Poisson distribution across biological replicates?

In a Poisson distribution, the variance is constrained to equal the mean (σ² = μ). While technical replicates (sequencing the same cDNA library across multiple flow cell lanes) follow Poisson sampling noise, biological replicates (independent cell cultures, mice, or human patients) introduce genuine inter-individual biological variability in transcription rates, resulting in overdispersion where σ² ≫ μ at moderate-to-high expression levels. Applying a Poisson test underestimates variance and produces thousands of false-positive differentially expressed genes; the Negative Binomial distribution adds a quadratic dispersion parameter (α_iμ²) that captures both technical shot noise (μ) and biological variance (α_iμ²).

Why does DESeq2 use the median-of-ratios (or edgeR use TMM) normalization instead of simple library read depth?

Read-count (or TPM) normalization assumes that cellular mRNA output is equal across conditions. However, if a small subset of abundantly expressed genes (such as hemoglobin in blood or tissue-specific secretory proteins) is upregulated in condition B, those transcripts consume a large fraction of the fixed sequencing read budget, causing all other non-changing housekeeping genes in condition B to appear artificially downregulated (composition bias). DESeq2's median-of-ratios and edgeR's Trimmed Mean of M-values (TMM) calculate sample scaling factors from the median expression ratio of non-differentially expressed majority genes, neutralizing composition bias.

How does long-read RNA-Seq (PacBio Iso-Seq and Oxford Nanopore direct RNA sequencing) overcome the isoform reconstruction limits of short-read RNA-Seq?

Short-read Illumina RNA-Seq breaks 2,000–10,000 nt mRNAs into 200–300 bp fragments; when a gene contains multiple distant alternatively spliced exons, short reads cannot link which 5′ exon combined with which 3′ exon on the same molecule, forcing statistical inference. PacBio Iso-Seq sequences full-length 1–10 kb cDNA molecules end-to-end via Circular Consensus (HiFi) reads, capturing full splice chains in a single read. Furthermore, Oxford Nanopore Direct RNA Sequencing threads native polyadenylated mRNA molecules directly through protein nanopores without reverse transcription or PCR, measuring poly(A) tail length and mapping epitranscriptomic modifications such as N⁶-methyladenosine (m⁶A) and pseudouridine (Ψ) on native transcripts.

1. Wang, Z., Gerstein, M., & Snyder, M. (2009). "RNA-Seq: A Revolutionary Tool for Transcriptomics." Nature Reviews Genetics, 10(1), 57–63. https://doi.org/10.1038/nrg2484
2. Mortazavi, A., Williams, B. A., McCue, K., Schaeffer, L., & Wold, B. (2008). "Mapping and Quantifying Mammalian Transcriptomes by RNA-Seq." Nature Methods, 5(7), 621–628. https://doi.org/10.1038/nmeth.1226
3. Dobin, A., Davis, C. A., Schlesinger, F., Drenkow, J., Zaleski, C., Jha, S., Batut, P., Chaisson, M., & Gingeras, T. R. (2013). "STAR: Ultrafast Universal RNA-seq Aligner." Bioinformatics, 29(1), 15–21. https://doi.org/10.1093/bioinformatics/bts635
4. Love, M. I., Huber, W., & Anders, S. (2014). "Moderated Estimation of Fold Change and Dispersion for RNA-seq Data with DESeq2." Genome Biology, 15(12), 550. https://doi.org/10.1186/s13059-014-0550-8