Modern Transcriptomics: From Bulk RNA-Seq to Single-Cell and Spatial Resolution

Introduction

Transcriptomics encompasses the comprehensive study of RNA transcripts, including their identity, abundance, and regulatory dynamics within a biological system. The field has evolved from low-throughput methods such as northern blotting and quantitative PCR to high-throughput sequencing-based approaches that capture the entire transcriptome. This review details the progression from bulk RNA sequencing (RNA-seq) through single-cell RNA-seq (scRNA-seq) to spatial transcriptomics, with a focus on the computational and biophysical principles underlying each technique. The discussion emphasizes applications in veterinary medicine and comparative biology, where transcriptomic profiling informs host-pathogen interactions, tissue-specific gene regulation, and diagnostic biomarker discovery.

Bulk RNA-Seq: Expression Quantification and Normalization

Bulk RNA-seq measures the average gene expression across a population of cells from a tissue sample. The core workflow begins with RNA extraction, poly-A selection or ribosomal RNA depletion, cDNA synthesis, and library preparation followed by high-throughput sequencing. The resulting reads are aligned to a reference genome or transcriptome using splice-aware aligners.

Quantification Strategies

Gene expression quantification from aligned reads proceeds at the gene, transcript, or exon level. The simplest metric is read count per gene, which sums all reads mapping to exonic regions. However, gene length and sequencing depth introduce biases that necessitate normalization. The three most common normalization methods for bulk RNA-seq are fragments per kilobase of transcript per million mapped reads (FPKM), transcripts per million (TPM), and reads per kilobase per million (RPKM). FPKM and RPKM normalize for gene length and library size but are not directly comparable across samples, whereas TPM normalizes for gene length after correcting for library size, providing better cross-sample comparability [1]. For differential expression analysis, count-based methods such as DESeq2 and edgeR model the raw count data using a negative binomial distribution, which accounts for the mean-variance relationship inherent in RNA-seq data [2]. These tools estimate dispersion parameters and apply statistical tests to identify differentially expressed genes.

Technical and Biological Considerations

Bulk RNA-seq provides a population average that obscures heterogeneity within the sample. This limitation is critical in veterinary contexts such as tumor biopsy analysis or immune cell infiltrate characterization, where distinct cell subpopulations may respond differently to infection or treatment. Deconvolution algorithms can infer cell type proportions from bulk data using reference signatures, but these methods lack single-cell resolution. The dynamic range of transcript detection in bulk RNA-seq spans approximately four to five orders of magnitude, with high-abundance transcripts dominating the signal. This dynamic range is constrained by the sequencing depth and the efficiency of library preparation.

Differential Expression Analysis

Differential expression (DE) analysis identifies genes whose expression levels change significantly between experimental conditions. The statistical framework must account for the discrete nature of count data and the presence of biological variability.

Statistical Models

The negative binomial distribution is the standard model for RNA-seq counts. DESeq2 estimates a dispersion parameter for each gene as a function of the mean expression level, using a shrinkage estimator to stabilize dispersion estimates for genes with low counts [2]. edgeR employs a similar approach with empirical Bayes methods to moderate the dispersion estimates across genes [3]. Both tools output fold change estimates and adjusted p-values, typically corrected for multiple testing using the Benjamini-Hochberg procedure. The selection of a significance threshold (commonly adjusted p-value < 0.05 and log2 fold change > 1) defines the set of differentially expressed genes. For an in-depth discussion of these algorithms, see RNA-Seq Differential Expression.

Confounding Factors

Batch effects, library preparation variability, and differences in sequencing depth are common confounders in DE analysis. Principal component analysis (PCA) of normalized count data is a standard diagnostic for identifying major sources of variation. ComBat-seq and RUVseq are computational methods that adjust for unwanted variation when batch information is known or unknown. In veterinary studies, factors such as animal age, sex, breed, and sample collection site can introduce systematic biases that must be modeled in the design matrix.

Single-Cell RNA-Seq: Preprocessing and Quality Control

Single-cell RNA-seq captures the transcriptome of individual cells, revealing cellular heterogeneity invisible to bulk methods. The technology relies on isolating single cells into nanoliter-volume reaction chambers, typically using microfluidic devices or droplet-based systems. Each cell is lysed, and its mRNA is reverse transcribed with a unique molecular identifier (UMI) and a cell-specific barcode.

Preprocessing Pipeline

The preprocessing pipeline for scRNA-seq includes read alignment, UMI deduplication, cell barcode assignment, and quality filtering. UMIs are random oligonucleotide sequences that label each mRNA molecule before amplification, enabling the removal of PCR duplicates by collapsing reads with identical UMIs and cell barcodes [4]. This step is critical for accurate quantification because amplification bias can skew transcript counts. After counting, a gene-cell expression matrix is constructed with rows representing genes and columns representing cells.

Quality Control Metrics

Low-quality cells are identified using three metrics: library size (total UMI counts per cell), number of expressed genes, and the fraction of reads mapping to mitochondrial genes. Cells with low library size or few genes likely represent empty droplets or dead cells. A high mitochondrial fraction (typically > 20%) indicates cellular stress or damage, as mitochondrial transcripts are released from damaged mitochondria. Doublets, which are two cells captured in a single droplet, show abnormally high gene counts and are removed using computational doublet detection tools.

Normalization and Batch Correction

Normalization in scRNA-seq must address the high dropout rate (zero inflation) and the large dynamic range of expression across cells. The standard approach is to compute size factors by pooling cells and using deconvolution to normalize library sizes [5]. Global scaling normalization, where each cell's counts are divided by its total count and multiplied by a scaling factor, is also used. For detailed methods, see Single-Cell RNA-Seq Normalization.

Batch correction is essential when integrating data from multiple sequencing runs or samples. Methods such as Harmony, Seurat's CCA (canonical correlation analysis), and Scanorama align cells in a shared low-dimensional space while preserving biological variation [6]. These tools identify mutual nearest neighbors across batches to correct for technical differences.

Dimensionality Reduction and Clustering

The high-dimensional gene expression matrix (tens of thousands of genes per cell) is reduced to a low-dimensional representation for visualization and clustering. PCA is applied first to capture the dominant axes of variation. The top 20-50 principal components are typically retained. For visualization, t-distributed stochastic neighbor embedding (t-SNE) and uniform manifold approximation and projection (UMAP) generate two-dimensional embeddings that preserve local and global structure [7]. UMAP is generally preferred for its faster runtime and better preservation of global topology.

Clustering partitions cells into groups with similar expression profiles. Graph-based clustering, implemented in Seurat and Scanpy, constructs a k-nearest neighbor graph on the PCA-reduced data and applies the Louvain or Leiden algorithm for community detection [8]. Marker genes for each cluster are identified by differential expression testing between clusters, enabling cell type annotation.

Trajectory Inference and Lineage Tracing

Trajectory inference algorithms order cells along a continuous path representing a biological process such as differentiation, cell cycle progression, or response to stimulus. Tools like Monocle, Slingshot, and fateID construct a minimum spanning tree or principal graph through the low-dimensional representation [9]. Pseudotime is a scalar value assigned to each cell that represents its position along the trajectory. These methods require highly expressed genes that change monotonically across the process. For a comprehensive review, see Single-cell RNA-seq Trajectory Inference and Cell Lineage Tracing.

Spatial Transcriptomics

Spatial transcriptomics preserves the spatial coordinates of RNA transcripts within a tissue section, providing context for cell-cell interactions and tissue architecture. Two major approaches exist: sequencing-based methods (e.g., 10x Visium, Slide-seq) and imaging-based methods (e.g., MERFISH, seqFISH+).

Sequencing-Based Spatial Transcriptomics

In sequencing-based methods, a tissue section is placed on a slide with spatially barcoded oligonucleotides. The tissue is permeabilized, allowing mRNA to diffuse and bind to the barcoded spots. After reverse transcription and library preparation, sequencing reads are assigned to spatial coordinates via the barcodes [10]. Each spot contains 10-100 cells, depending on the spot size. The resulting data is a spatially resolved gene expression matrix with lower resolution than single-cell methods but with intact spatial information.

Imaging-Based Spatial Transcriptomics

Imaging-based methods use fluorescence in situ hybridization (FISH) with sequential rounds of hybridization to detect individual RNA molecules. MERFISH uses error-robust barcodes and achieves subcellular resolution [11]. seqFISH+ extends the multiplexing capacity to several thousand genes through sequential hybridization and combinatorial barcoding. These methods provide single-cell resolution but are limited in transcript detection efficiency and throughput.

Computational Analysis of Spatial Data

Spatial transcriptomics data analysis involves detecting spatial patterns of gene expression, mapping cell types to spatial coordinates, and identifying ligand-receptor interactions. Non-negative matrix factorization can identify spatially coherent gene expression programs. Cell type deconvolution assigns cell type proportions to each spatial spot using reference scRNA-seq data. For alignment of multiple tissue sections and neighborhood analysis, see Spatial Transcriptomics Alignment and Cellular Neighborhood Analysis.

Workflow Summary: From Bulk to Spatial

The following diagram summarizes the key computational steps across the three transcriptomic modalities.

flowchart TD
    A[Start: RNA Sample], > B[Bulk RNA-Seq]
    A, > C[Single-Cell RNA-Seq]
    A, > D[Spatial Transcriptomics]
    
    B, > B1[Read Alignment & Quantification]
    B1, > B2[Count Normalization: TPM/FPKM]
    B2, > B3[Differential Expression: DESeq2/edgeR]
    B3, > B4[Pathway Enrichment & Interpretation]
    
    C, > C1[Cell Barcode & UMI Processing]
    C1, > C2[Quality Control: Library Size, MT Fraction]
    C2, > C3[Normalization & Batch Correction]
    C3, > C4[Dimensionality Reduction: PCA + UMAP]
    C4, > C5[Clustering & Cell Type Annotation]
    C5, > C6[Trajectory Inference & Differential Testing]
    
    D, > D1[Spatial Barcode Assignment]
    D1, > D2[Gene Expression Matrix with Coordinates]
    D2, > D3[Cell Type Deconvolution]
    D3, > D4[Spatial Pattern Detection]
    D4, > D5[Ligand-Receptor Interaction Analysis]

Integration of Multi-Omics Data

Transcriptomic data is increasingly integrated with other omic layers such as genomics, epigenomics (e.g., scATAC-seq), and proteomics. Single-cell multi-omic technologies enable the simultaneous measurement of gene expression and chromatin accessibility from the same cell [12]. Computational integration methods, including weighted nearest neighbor analysis and matrix factorization, combine these data modalities to identify regulatory relationships. For details on chromatin accessibility profiling, see Single-Cell ATAC-Seq Bioinformatics.

Applications in Veterinary Medicine

In veterinary research, transcriptomics has been applied to understand immune responses to viral and bacterial pathogens, characterize tumor microenvironments in canine and feline cancers, and map tissue-specific gene expression in livestock species. Spatial transcriptomics is particularly valuable for studying infectious disease histopathology, where the spatial organization of infected cells and immune infiltrates influences disease progression. For host-pathogen interaction studies, see Single-Cell Transcriptomics of Host-Pathogen Interactions During Viral Infection.

Computational Infrastructure

Processing transcriptomic data requires robust computational infrastructure. Alignment, quantification, and downstream analysis for a typical scRNA-seq experiment with 10,000 cells require 16-64 GB of RAM and multiple CPU cores. Cloud computing resources and workflow management systems like Snakemake and Nextflow provide scalable and reproducible processing [13]. For a guide to workflow environments, see Cloud Computing in Modern Bioinformatics.

Limitations and Future Directions

Current transcriptomic methods face several limitations. Bulk RNA-seq sacrifices cellular resolution. scRNA-seq captures only a small fraction of total mRNA per cell (typically 5-20%) and requires cell dissociation, which can alter gene expression. Spatial methods have trade-offs between resolution, throughput, and gene detection efficiency. Future developments include the integration of transcriptomics with proteomics at single-cell resolution, improved spatial resolution through expansion microscopy, and the application of machine learning for gene regulatory network inference [14]. For deep learning approaches to network reconstruction, see Deep Learning for Gene Regulatory Network Reconstruction.

References

[1] Wagner GP, Kin K, Lynch VJ. Measurement of mRNA abundance using RNA-seq data: RPKM measure is inconsistent among samples. Theory Biosci. 2012;131(4):281-285.

[2] Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.

[3] Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139-140.

[4] Islam S, Zeisel A, Joost S, et al. Quantitative single-cell RNA-seq with unique molecular identifiers. Nat Methods. 2014;11(2):163-166.

[5] Lun ATL, Bach K, Marioni JC. Pooling across cells to normalize single-cell RNA sequencing data with many zero counts. Genome Biol. 2016;17:75.

[6] Korsunsky I, Millard N, Fan J, et al. Fast, sensitive and accurate integration of single-cell data with Harmony. Nat Methods. 2019;16(12):1289-1296.

[7] Becht E, McInnes L, Healy J, et al. Dimensionality reduction for visualizing single-cell data using UMAP. Nat Biotechnol. 2019;37(1):38-44.

[8] Traag VA, Waltman L, van Eck NJ. From Louvain to Leiden: guaranteeing well-connected communities. Sci Rep. 2019;9(1):5233.

[9] Trapnell C, Cacchiarelli D, Grimsby J, et al. The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nat Biotechnol. 2014;32(4):381-386.

[10] Ståhl PL, Salmén F, Vickovic S, et al. Visualization and analysis of gene expression in tissue sections by spatial transcriptomics. Science. 2016;353(6294):78-82.

[11] Chen KH, Boettiger AN, Moffitt JR, et al. Spatially resolved, highly multiplexed RNA profiling in single cells. Science. 2015;348(6233):aaa6090.

[12] Cao J, Cusanovich DA, Ramani V, et al. Joint profiling of chromatin accessibility and gene expression in thousands of single cells. Science. 2018;361(6409):1380-1385.

[13] Köster J, Rahmann S. Snakemake-a scalable bioinformatics workflow engine. Bioinformatics. 2012;28(19):2520-2522.

[14] Efremova M, Vento-Tormo M, Teichmann SA, et al. CellPhoneDB: inferring cell-cell communication from combined expression of multi-subunit ligand-receptor complexes. Nat Protoc. 2020;15(4):1484-1506. *** Disclaimer: This article is for educational and informational purposes only. It is not intended to substitute for professional veterinary advice, diagnosis, treatment, or regulatory guidance. Always consult a licensed veterinarian or qualified specialist regarding animal health, disease diagnosis, and therapeutic decisions.