7  Expression analysis

NoteFrom the exam description

Expression analysis of genes and proteins: methods, algorithms, statistics.

7.1 Overview

Measuring which genes are expressed, and how much, links the genome to the phenotype. Expression analysis was the first area of genomics that produced large, quantitative data sets about many genes in many samples, and it forced biologists to confront statistical problems that are now everywhere: normalisation, multiple testing, batch effects and the exploration of high-dimensional data. This chapter follows the development from microarrays to single-cell and spatial methods:

• microarrays, their preprocessing and what they taught us (Section 7.2);
• RNA-seq: library preparation and experimental design, spliced alignment, transcript assembly, quantification and normalisation (Section 7.3);
• statistics: differential expression with count models, multiple testing, clustering and its validation, classification with many genes and few samples, and dimensionality reduction (Section 7.4);
• interpretation: gene set analysis, co-expression and eQTLs, and the relation between mRNA and protein levels (Section 7.5);
• single-cell genomics: technologies, analysis, multimodal data, spatial transcriptomics, foundation models and cell types (Section 7.6);
• single-cell sequencing of cancer and tumour evolution (Section 7.7).

The databases that archive expression data are described in Section 2.8, and the measurement of proteins in Chapter 8.

7.2 The first genome-wide methods: microarrays

Measuring which genes are expressed, and how much, links the genome to the phenotype. Before the 1990s, expression was measured one gene at a time, for example with Northern blots. Two approaches made it possible to measure all genes at once: counting short sequence tags, as in serial analysis of gene expression (SAGE) (Velculescu et al. 1995), and hybridisation to DNA microarrays. Microarrays dominated for about fifteen years, and they introduced most of the ideas that are still used to analyse expression data.

7.2.1 Principle

A DNA microarray is a solid surface with thousands to millions of spots (features), each containing many copies of a DNA probe for one transcript. RNA from a sample is converted to labelled cDNA or cRNA and hybridised to the array. Each labelled molecule binds to the probe with the complementary sequence, and the fluorescence intensity of each spot is a measure of how abundant that transcript was in the sample. The relation is monotonic but not simply proportional: background dominates at low abundance, the signal saturates at high abundance, and each probe has its own binding affinity. Intensities are therefore relative measures, best compared for the same probe between samples. Three main technologies were used:

• Spotted cDNA arrays, developed in Patrick Brown’s laboratory at Stanford, were made by robotically printing PCR products onto glass slides (Schena et al. 1995). They were usually used in two-colour experiments: two samples, labelled with different dyes (Cy3 and Cy5), were hybridised to the same array, and the result was the ratio between them. Because they could be built in academic laboratories, they spread rapidly.
• Affymetrix GeneChips were made by synthesising short (25-nucleotide) oligonucleotides directly on the chip with photolithography, the technique used to make computer chips (Fodor et al. 1991; Lockhart et al. 1996). Each gene was measured with a probe set of about 11–20 different probes, each with a mismatch control probe that differed in the central base. These were one-colour arrays, each measuring one sample, and the high reproducibility of the industrial process made them the most widely used platform.
• Illumina BeadArrays use a different principle. Probes are attached to microscopic beads, which are randomly distributed in wells etched in a fibre-optic bundle or a glass slide. Since the position of each bead is random, the array is decoded after manufacture by a series of hybridisations that identify which probe each bead carries (Gunderson et al. 2004). Each probe is present on many beads, which gives built-in replication. BeadArrays were used for expression profiling (Kuhn et al. 2004), and, above all, for SNP genotyping: the genotyping chips used in most genome-wide association studies, measuring up to millions of SNPs per sample, are Illumina bead arrays.

7.2.2 Preprocessing and normalisation

The raw data from an array are images, and many steps are needed before intensities can be compared between genes and samples (Quackenbush 2001):

1. Image analysis: locating the spots and estimating their intensity and the local background.
2. Background correction: removing the signal from non-specific hybridisation and the optical background.
3. Transformation: working on the logarithmic scale, usually \(\log_2\), because intensities span several orders of magnitude, errors are roughly multiplicative, and fold changes become differences.
4. Normalisation: removing systematic differences between arrays or dyes, caused by differences in the amount of RNA, labelling efficiency or scanner settings. For two-colour arrays, the log-ratio \(M = \log_2(R/G)\) is plotted against the average log intensity \(A\) (an MA plot), and a curve fitted with local regression (loess) is subtracted, so that most genes have \(M\) close to zero. For one-colour arrays, quantile normalisation forces all arrays to have the same distribution of intensities (Bolstad et al. 2003).
5. Summarisation: for Affymetrix arrays, combining the probes of a probe set into one expression value. The robust multi-array average (RMA) fits a model across arrays that accounts for probe-specific affinities, and ignores the mismatch probes, which turned out to measure specific signal as well as background (Irizarry 2003).

These steps are not just technicalities. Normalisation methods assume that most genes do not change between samples, an assumption that fails in some experiments, and the choice of method can change the results.

7.2.3 What microarrays showed, and their limitations

Microarrays made it possible to ask questions about whole transcriptomes. Clustering genes by their expression profiles showed that genes with related functions are co-expressed (Eisen et al. 1998). Expression profiles could classify cancers: acute leukaemias could be told apart by expression alone (Golub et al. 1999), breast tumours fell into molecular subtypes with different prognoses (Perou et al. 2000), and a 70-gene signature predicted the outcome of breast cancer (Veer et al. 2002) and became a clinical test. Atlases of expression across tissues showed which genes are ubiquitous and which are tissue specific (Su et al. 2004). The analysis of these experiments raised the statistical problems that are discussed in Section 7.4: how to test thousands of genes with few replicates, for which empirical Bayes methods such as limma were developed (Smyth 2004), and how to control false discoveries.

The early years also showed how easily such experiments go wrong. Results were hard to reproduce between laboratories and platforms until standards for reporting experiments, MIAME (Brazma et al. 2001), and large comparison studies (MAQC Consortium 2006) were established. Microarrays have several built-in limitations:

• they can only measure transcripts whose sequences are known and represented on the array;
• signals are affected by cross-hybridisation of similar sequences and by probe-specific affinities, so intensities of different genes cannot be compared directly;
• the dynamic range is limited by background at the low end and saturation at the high end;
• they cannot easily distinguish isoforms, detect new transcripts or measure allele-specific expression.

RNA-seq (Section 7.3) removed most of these limitations, and replaced microarrays for expression profiling from about 2010. Arrays remain in use for genotyping and DNA methylation profiling, where they are cheap and robust.

See also: further explanation in Section D.7.1; related Wikipedia articles in Section E.7.1.

NoteSources

Parts of this section are adapted from the Wikipedia articles DNA microarray, Gene expression profiling, Serial analysis of gene expression, Quantile normalization and MA plot (CC BY-SA 4.0, retrieved 4 October 2026).

7.3 RNA sequencing

RNA-seq measures expression by sequencing cDNA made from the RNA in a sample and counting the reads that come from each gene or transcript. It was introduced in 2008 (Mortazavi et al. 2008; Nagalakshmi et al. 2008), shown to agree well with microarrays while being more reproducible and sensitive (Marioni et al. 2008), and quickly became the standard (Wang et al. 2009; Stark et al. 2019). Because it does not depend on predefined probes, RNA-seq can discover new genes and transcripts, measure isoforms, detect allele-specific expression, fusion genes and RNA editing, and be used in any organism, with or without a reference genome. Its dynamic range is limited only by the sequencing depth.

7.3.1 Library preparation and experimental design

Most of the RNA in a cell is ribosomal. To sequence mRNA, either polyadenylated RNA is selected with oligo(dT), or rRNA is depleted, which also keeps non-polyadenylated RNAs. The RNA is fragmented, reverse transcribed into cDNA, ligated to adapters and amplified. Strand-specific protocols record which strand the transcript came from, which matters where genes overlap on opposite strands. Typical experiments use 20–50 million short reads per sample for gene-level analysis. Long-read sequencing (Section 5.2) reads entire transcripts and is increasingly used to define isoforms (see below).

The most important design decision is the number of biological replicates. The key concept is the experimental unit: the smallest unit that is independently assigned to (or sampled from) a condition, such as an animal, a patient or an independently grown culture. Biological replicates are different experimental units in the same condition, and they capture the natural variation between individuals. Technical replicates are repeated measurements of the same unit, for example the same RNA sequenced in two libraries or two lanes; they measure only the noise of the protocol. Technical replicates show little variation, while biological replicates vary much more, and differential expression can only be assessed relative to biological variation. Treating technical replicates, or several samples from the same individual, as independent biological replicates is called pseudoreplication and gives far too many significant results. Studies have shown that with three biological replicates many true changes are missed, and that six or more are needed for good sensitivity (Schurch et al. 2016).

Good design also controls confounding. Technical batch effects (different days, operators, reagent lots or sequencing runs) can be larger than the biological differences of interest (Leek et al. 2010). If all control samples are processed in one batch and all treated samples in another, batch and condition are confounded, and no statistical method can separate them. The remedies are classical: randomise the assignment of samples to conditions and to processing order; block, by processing equal numbers of samples from each condition in every batch, so that batch can be included as a covariate in the model (Section 7.4); and use paired designs where possible, for example tumour and normal tissue from the same patient, or before and after treatment in the same individual, so that the comparison is made within each individual and differences between individuals cancel out. Such design decisions must be made before the experiment; they cannot be repaired in the analysis.

7.3.2 The analysis workflow

A standard analysis has the following steps (Conesa et al. 2016):

1. Quality control of the reads: base qualities, adapter content, duplication and rRNA contamination.
2. Alignment or quantification, described below.
3. Counting reads per gene, or estimating abundances per transcript.
4. Normalisation for sequencing depth and composition.
5. Exploratory analysis, such as PCA and clustering, to detect outliers and batch effects (Section 7.4).
6. Statistical testing for differential expression, and interpretation with gene sets and pathways (Section 7.5).

Each step rests on assumptions, and an observed difference can have several explanations besides regulation within cells: a change in cell composition, poor RNA quality, protocol bias, errors in the annotation or a batch effect (Figure 7.1). A longer gene list does not repair a flawed design.

Figure 7.1: An expression result as a chain of evidence, from the biological question to validation. Each step has its own question that must be answered before the next can be trusted. Figure: Matteo Tadiello, Modern Bioinformatics: From Foundations to Frontiers (CC BY 4.0).

Quality control

Quality control is a diagnosis of how the data were generated, not an automatic filter. Before alignment, base qualities, adapter content, over-represented sequences and per-cycle base composition are checked (Section 5.4). After alignment, the useful statistics are specific to RNA: the fraction of uniquely mapped reads, the fractions of ribosomal and mitochondrial reads, strand specificity, splice-junction support and the coverage along gene bodies. Degraded RNA gives a characteristic 3′ bias, since fragments near the poly(A) tail survive best. At the level of samples, PCA, sample correlations and clustering show whether samples group by condition or by processing batch, and genotypes called from the reads can confirm that each sample comes from the expected individual. An outlying sample is a reason to investigate, not to exclude: exclusion rules should be decided before the results are known, and conclusions checked with and without doubtful samples. The SEQC/MAQC-III consortium sequenced reference RNA mixtures at many sites and showed that platform, site, expression level and analysis pipeline all affect the results, while relative expression measures were reproducible across sites (“A Comprehensive Assessment of RNA-Seq Accuracy, Reproducibility and Information Content by the Sequencing Quality Control Consortium” 2014).

Spliced alignment

RNA-seq reads come from spliced transcripts, so a read that spans an exon–exon junction aligns to the genome in two parts separated by an intron, which may be tens of kilobases long. Ordinary DNA read mappers (Section 5.4) cannot handle this. TopHat first mapped reads that align without gaps, predicted junctions from them and then aligned the remaining reads across those junctions (Trapnell et al. 2009). STAR searches for the longest prefix of the read that matches the genome exactly (the maximal mappable prefix) using an uncompressed suffix array, then repeats the search for the rest of the read, and stitches the pieces together (Dobin et al. 2012). It is very fast but needs about 30 GB of memory for the human genome. HISAT2 uses a hierarchical graph FM-index with much smaller memory requirements (Kim et al. 2019).

Transcript assembly

Alignments can be used to reconstruct the transcripts themselves. Reference-guided assembly, with Cufflinks (Trapnell et al. 2010) or StringTie (Pertea et al. 2015), builds a graph of overlapping alignments for each gene and finds a small set of transcripts that explain the reads, which can reveal new isoforms. Cufflinks used a minimum path cover of the overlap graph, the most parsimonious set of transcripts. In organisms without a reference genome, transcripts are assembled de novo with de Bruijn graph methods such as Trinity (Grabherr et al. 2011) (Section 5.3).

Quantification and pseudo-alignment

The simplest quantification counts the reads that overlap each gene. Estimating the abundance of individual transcripts is harder, because isoforms of the same gene share most of their exons, so many reads are compatible with several transcripts. RSEM solves this with the expectation–maximisation (EM) algorithm: it alternates between assigning each read fractionally to transcripts in proportion to their current abundance estimates, and re-estimating the abundances from these assignments (Li and Dewey 2011).

Aligning every read base by base is not necessary for quantification. It is enough to know which transcripts a read is compatible with. kallisto introduced pseudo-alignment: it builds a de Bruijn graph of the transcriptome, looks up the k-mers of each read, and intersects the sets of transcripts that contain them (Bray et al. 2016). Reads with the same set of compatible transcripts form an equivalence class, and the EM algorithm then runs on the counts of these classes. Salmon uses a similar selective alignment and models sequence-specific, fragment-GC and positional biases (Patro et al. 2017). These methods are one to two orders of magnitude faster than alignment and similarly accurate for known transcripts, but they cannot discover new ones.

Long-read RNA-seq

Short reads cover only a small part of a transcript, so the full isoform must be inferred by assembly or by EM, and different combinations of exons far apart in the same gene are hard to resolve. Long-read platforms (PacBio and Oxford Nanopore; Section 5.2) can sequence complete cDNA molecules, and nanopore sequencing can also read native RNA directly. Each read then shows one full-length isoform, with its exact combination of exons, its start and end sites and its poly(A) tail, so no assembly is needed. Reads are aligned to the genome with spliced aligners such as minimap2 (Li 2018), and grouped into known and new transcripts. This has made long reads the method of choice for transcript discovery and for building isoform catalogues, and has shown that many genes have more isoforms than annotated (Pardo-Palacios et al. 2024). The drawbacks are fewer reads per sample, higher error rates on some platforms (notably direct RNA sequencing), which make exact splice sites harder to call, truncated reads from degraded RNA, and a bias towards shorter molecules. Short reads are therefore still often used for quantification, sometimes in combination with long reads.

7.3.3 Normalisation

Raw counts are not comparable between genes or samples. Longer transcripts produce more fragments, and samples sequenced more deeply produce more reads. A gene with twice as many reads as another is therefore not necessarily twice as highly expressed. Several units correct for these effects:

• CPM (counts per million) divides by the library size and corrects for depth only.
• RPKM/FPKM (reads or fragments per kilobase of transcript per million mapped reads) also divides by the transcript length (Mortazavi et al. 2008).
• TPM (transcripts per million) first divides each count by the transcript length, and then scales so that the values in each sample sum to one million: \[ \text{TPM}_i = \frac{c_i / \ell_i}{\sum_j c_j / \ell_j} \times 10^6, \] where \(c_i\) is the count and \(\ell_i\) the effective length of transcript \(i\). TPM estimates the relative molar concentration of each transcript and, unlike RPKM, sums to the same total in every sample (Wagner et al. 2012).

These units do not handle composition bias. RNA-seq data are compositional: a library contains a fixed number of reads sampled from the RNA in it, so the counts measure each gene’s share of the library, not its absolute amount. If a few genes are very highly expressed in one sample, they take up a large share of the reads, and all other genes appear lower even if their true expression is unchanged. Methods for differential expression therefore estimate normalisation factors under the assumption that most genes are not differentially expressed:

• TMM (trimmed mean of M-values), used by edgeR, computes the weighted mean of log-ratios between a sample and a reference sample after trimming the most extreme genes (Robinson and Oshlack 2010).
• The median-of-ratios method of DESeq computes, for each gene, the ratio of its count to the geometric mean of its counts across all samples. The median of these ratios in a sample is its size factor (Anders and Huber 2010; Love et al. 2014).

This assumption can fail. If the total amount of mRNA per cell changes, for example when a transcription factor such as MYC raises the expression of most genes, every gene’s share of the library stays about the same, and all of these normalisations wrongly conclude that little has changed (Lovén et al. 2012). Detecting such global shifts requires an external reference, such as synthetic spike-in RNAs added in proportion to the number of cells.

TPM values are useful for comparing the relative abundance of transcripts within a sample. Because TPMs sum to the same total in every sample, they are still proportions, and they do not remove composition effects between samples; they are suitable for descriptive comparisons but are not a general normalisation. For differential expression, tools such as DESeq2 and edgeR take the raw (unnormalised) integer counts as input, not TPM or CPM values, and include the normalisation factors as offsets in the statistical model (Section 7.4).

See also: further explanation in Section D.7.2; related Wikipedia articles in Section E.7.2.

NoteSources

Parts of this section are adapted from the Wikipedia articles RNA-Seq, Transcriptomics technologies and Expectation–maximization algorithm (CC BY-SA 4.0, retrieved 4 October 2026). Parts of this section are adapted from Matteo Tadiello, Modern Bioinformatics: From Foundations to Frontiers, chapter “Gene-Expression Analysis” (CC BY 4.0), used with permission.

7.4 Testing and statistics

Expression studies measure tens of thousands of genes in a handful of samples. This raises three statistical problems that appear throughout genomics: modelling the variability of the measurements with few replicates, testing thousands of hypotheses at once, and exploring high-dimensional data. The section starts with the general statistical background that is expected of every student, whatever their topic.

7.4.1 Statistical background

Most bioinformatics results end with a statistical statement: an E-value, a p-value for differential expression, a confidence in a phylogenetic split, an accuracy on a test set. You should be comfortable with the basic concepts and know their common pitfalls. The Nature Methods Points of Significance series (Statistics for Biologists) covers all of this in two-page articles and is the recommended starting point (Krzywinski and Altman 2013).

Probability and distributions

You should be able to explain conditional probability and Bayes’ theorem, \[P(H \mid D) = \frac{P(D \mid H)\, P(H)}{P(D)},\] and use it on a simple example. A classic one: a test that detects a rare condition (1 in 1,000) with 99% sensitivity and 99% specificity still gives mostly false positives, because the few true cases are outnumbered by the 1% false positives among the many healthy people. The same base-rate argument explains why E-values must account for database size (Section 1.4) and why most predicted interactions or binding sites in a genome-wide screen are false.

Know the common distributions and when they arise: the binomial (number of successes in \(n\) trials), the Poisson (rare events, read counts at low coverage), the normal (sums of many small effects), the negative binomial (overdispersed counts, as in RNA-seq; Chapter 7) and the extreme value distribution (maxima, as in alignment scores; Section 1.4). Understand likelihood, the probability of the data given a model, and maximum-likelihood estimation. It underlies phylogenetics (Chapter 4), HMM training (Section 1.6) and much of machine learning.

Hypothesis tests and p-values

A statistical test compares the data with a null hypothesis. The p-value is the probability, assuming the null hypothesis is true, of observing data at least as extreme as those seen. It is not the probability that the null hypothesis is true, and a large p-value does not show that the null is true. A small p-value says nothing about how large or important an effect is: with enough data, a trivial difference becomes “significant”. Report effect sizes with confidence intervals, not only p-values (Krzywinski and Altman 2013; Nuzzo 2014; Wasserstein and Lazar 2016).

You should know the most common tests and their assumptions:

• the t-test compares means. The one-sample test compares a mean with a fixed value, the two-sample test compares two groups, and the paired test compares two measurements on the same units. It assumes approximately normal data (or large samples) and, in its classic form, equal variances; Welch’s t-test drops the equal-variance assumption and is usually the better default. The test statistic is the difference in means divided by its standard error;
• non-parametric alternatives, such as the Mann–Whitney (Wilcoxon rank-sum) and Wilcoxon signed-rank tests, use ranks and make fewer assumptions, at some cost in power;
• χ² and Fisher’s exact tests for counts in contingency tables, for example enrichment of a GO term in a gene list (Chapter 10);
• correlation (Pearson for linear relationships, Spearman for monotonic ones), and why correlation does not imply causation;
• permutation tests, which build the null distribution by shuffling labels, and are useful when no standard test fits.

Know also the difference between type I errors (false positives, controlled by the significance level \(\alpha\)) and type II errors (false negatives), and what determines statistical power: effect size, variability and sample size. Too few replicates is the most common reason for unreliable conclusions in genomics experiments.

7.4.2 Differential expression

The basic question is which genes change in expression between conditions. A t-test for each gene works poorly. With three replicates per group, the variance of each gene is estimated very imprecisely, and genes that happen to have a small sample variance get spuriously large t-statistics.

Moderated statistics. The solution, developed for microarrays, is to borrow information across genes. In the limma framework (Smyth 2004), an empirical Bayes model treats the true variances of all genes as drawn from a common prior distribution, estimated from all genes together. Each gene’s variance is then shrunk towards the common value. The resulting moderated t-statistic is much more stable, and limma remains in use, also for RNA-seq data after a transformation that estimates the mean–variance relationship (voom) (Law et al. 2014).

Count models. RNA-seq produces counts. If the same library were sequenced repeatedly, the count of a gene would follow a Poisson distribution, with variance equal to the mean (Marioni et al. 2008). Between biological replicates, the variance is larger (overdispersion). The standard model is the negative binomial distribution, with variance

\[ \operatorname{Var}(K) = \mu + \alpha \mu^2, \]

where \(\mu\) is the mean and \(\alpha\) the gene-specific dispersion. For highly expressed genes, the second term dominates, and \(\sqrt{\alpha}\) is the biological coefficient of variation. edgeR (Robinson et al. 2009) and DESeq2 (Love et al. 2014) fit a generalised linear model for each gene, with a log link,

\[ \log \mu_{ij} = \log s_j + \sum_r x_{jr} \beta_{ir}, \]

where \(s_j\) is the normalisation factor of sample \(j\) (Section 7.3), \(x_{jr}\) the design matrix describing the samples (condition, batch, paired individuals) and \(\beta_{ir}\) the coefficients, such as the log fold change between conditions. As in limma, the dispersion of each gene is estimated with shrinkage towards a trend fitted across all genes. DESeq2 also shrinks the fold changes, so that large fold changes of low-count, noisy genes are not over-interpreted. Significance is assessed with Wald or likelihood-ratio tests. Including known batches and other covariates in the design matrix is the main way to deal with batch effects (Leek et al. 2010).

Figure 7.2: The negative binomial model of RNA-seq counts. The expected count of each gene in each sample combines a sample-specific offset for depth and composition with the design matrix and coefficients; the observed count varies around it with a gene-specific dispersion; a contrast of the coefficients answers the biological question. Figure: Matteo Tadiello, Modern Bioinformatics: From Foundations to Frontiers (CC BY 4.0).

The tests answer slightly different questions. A Wald test divides an estimated coefficient, or a combination of coefficients, by its standard error. A likelihood-ratio test compares a full model with a reduced model without the terms of interest, and is the natural choice when the question involves several coefficients at once, such as any difference between several time points. The quasi-likelihood F-test of edgeR also takes into account the uncertainty in the estimated dispersion, and gives better control of false positives when there are few replicates (Lund et al. 2012). The test should follow from the question, not be chosen because it gives the most significant genes.

7.4.3 Design matrices, contrasts and interactions

The design matrix has one row per sample and one column per coefficient, and it states which comparisons the model can make. In a paired design, such as treated and untreated samples from the same donors, the model ~ donor + treatment includes one coefficient per donor, so that the treatment effect is estimated within donors, and stable differences between donors are removed. A known batch is added in the same way, as long as it is not confounded with the condition. If two columns of the design matrix are identical, for example when all treated samples were processed in one batch, the matrix is not of full rank and the two effects cannot be separated. Software that refuses to fit such a model is doing the right thing.

A contrast is a linear combination of coefficients, \(c^T\beta\), that answers one question. In an experiment with two genotypes, A and B, each with and without treatment, the model ~ genotype + treatment + genotype:treatment includes an interaction term. The interaction is a difference of differences:

\[ (\text{B}_\text{treated} - \text{B}_\text{control}) - (\text{A}_\text{treated} - \text{A}_\text{control}), \]

and asks whether the response to treatment depends on the genotype. It is not the treatment effect in either genotype. A common error is to test the treatment in each genotype separately and conclude that the genotypes respond differently because the effect is significant in one and not in the other. The difference between “significant” and “not significant” is not itself significant (Gelman and Stern 2006); a non-significant result may simply reflect low power. Such errors were found in about half of the neuroscience papers that compared two effects in this way (Nieuwenhuis et al. 2011). The correct analysis tests the interaction directly. Contrasts should be written down, in words and as coefficients, before the analysis.

7.4.4 Differential transcript usage

Differential expression of a gene is not the only kind of change. A gene can switch between isoforms without changing its total expression: differential transcript usage (or exon usage) asks whether the proportions of the isoforms differ between conditions. This requires transcript-level quantification (Section 7.3), whose estimates are uncertain when isoforms share most of their sequence. Quantifiers such as Salmon and kallisto estimate this uncertainty by bootstrapping or Gibbs sampling, and transcript-level estimates can be summarised to genes for gene-level tests, which also improves those tests when isoform usage changes (Soneson et al. 2016). A claim about isoforms should only be made when the reads, for example unique exons or junctions, can actually distinguish them.

7.4.5 Multiple testing

If 20,000 genes are tested at a significance level of 0.05, about 1,000 will be called significant even if nothing changes. Two kinds of error rates are controlled:

• The family-wise error rate (FWER) is the probability of at least one false positive. The Bonferroni correction controls it by testing each of \(m\) hypotheses at level \(\alpha/m\). Holm’s step-down method (Holm 1979) is uniformly more powerful and equally valid. FWER control is appropriate when any false positive is costly, but it is very conservative for genome-scale data (Noble 2009).
• The false discovery rate (FDR) is the expected proportion of false positives among the genes called significant (Benjamini and Hochberg 1995). It is the natural quantity when the result is a list of candidates for follow-up. The Benjamini–Hochberg procedure sorts the \(m\) p-values, \(p_{(1)} \le \dots \le p_{(m)}\), finds the largest \(k\) for which \(p_{(k)} \le \frac{k}{m}\alpha\), and rejects the hypotheses \(1, \dots, k\). This controls the FDR at level \(\alpha\) when the tests are independent or positively correlated.

The FDR is a property of the whole list, and an expected value over repeated experiments, not a statement about individual genes. A list called at FDR 5% is expected to contain about 5% false discoveries, but it does not mean that each gene in it has a 5% chance of being a false positive: genes near the cutoff are much more likely to be false than those with very small p-values.

Storey’s q-value is the smallest FDR at which a gene would be called significant (Storey and Tibshirani 2003). It also estimates the proportion of true null hypotheses, \(\pi_0\), from the flat part of the p-value histogram, which gives more power when many genes change. Efron’s empirical Bayes approach estimates the local false discovery rate, the probability that a particular gene is null given its test statistic (Efron and Tibshirani 2002). The histogram of p-values is a useful diagnostic: under the null hypothesis p-values are uniformly distributed, so a peak near zero indicates true signal, while other shapes indicate problems with the model.

Statistical significance is not biological importance. The p-value depends on the effect size, the variability and the number of replicates, so a large fold change in a noisy, lowly expressed gene may not be significant, while with enough replicates very small fold changes become significant. A volcano plot, the \(-\log_{10}\) p-value against the log fold change, shows both at once. Results are often filtered by both FDR and fold change, but there is no universal cutoff (such as a two-fold change): what matters depends on the gene and the biological question. Filtering on the estimated fold change after testing also undermines the FDR guarantee; it is better to test directly against a minimal relevant fold change, which both limma and DESeq2 allow (McCarthy and Smyth 2009; Love et al. 2014).

Removing genes before testing can increase power, because fewer tests are corrected for. Independent filtering removes genes with too few reads to give a significant result, using a statistic, such as the mean count over all samples, that is independent of the test statistic when the null hypothesis is true (Bourgon et al. 2010). DESeq2 does this automatically. A filter that uses the condition labels, such as keeping only genes that look different between the groups, is not independent and invalidates the FDR.

E-values in sequence searching (Section 1.4) are a built-in multiple-testing correction: they count the expected number of chance hits in the whole database. Be aware also of the hidden forms of multiple testing: trying several analyses and reporting the one that “worked”, testing many subgroups, or choosing thresholds after seeing the data. These practices are a major reason why many published findings fail to replicate (Ioannidis 2005).

7.4.6 Clustering and dimensionality reduction

Expression data are high-dimensional: each sample is a point in a space with one dimension per gene, and each gene a point in a space with one dimension per sample. Unsupervised methods find structure in these data without using the sample labels.

Clustering groups genes or samples with similar profiles. It requires a distance measure, such as Euclidean distance or one minus the correlation, which groups genes that go up and down together regardless of their absolute level.

• Hierarchical clustering repeatedly joins the two closest clusters, as in UPGMA (Section 4.3), and produces a tree (dendrogram). The distance between clusters is defined by the linkage: single (closest members), complete (most distant members) or average. Eisen’s heat maps of hierarchically clustered genes and samples became the iconic display of expression data (Eisen et al. 1998).
• k-means partitions the data into \(k\) clusters by alternating between assigning each point to the nearest cluster centre and recomputing the centres. The number of clusters must be chosen, and the result depends on the starting point.
• Networks of co-expressed genes, built by thresholding correlations, are an alternative to clusters, and modules of co-expressed genes can be related to traits (Langfelder and Horvath 2008).

Principal component analysis (PCA) finds the orthogonal directions of largest variance in the data, computed from the singular value decomposition of the centred data matrix (Lever et al. 2017). Plotting samples on the first two or three components shows the main sources of variation. If samples group by processing date rather than by condition, there is a batch effect. PCA is linear, so the first components keep the large-scale structure of the data, and the fraction of variance explained by each component tells how much of the data the plot shows.

t-SNE (Maaten and Hinton 2008) and UMAP (McInnes et al. 2018; Becht et al. 2018) are non-linear methods that place similar points close together in two dimensions. They preserve local neighbourhoods and reveal clusters, which has made them the standard way to visualise single-cell data (Section 7.6). However, distances between clusters, cluster sizes and densities in such plots are not meaningful, and the layout depends on parameters and random initialisation. Two clusters that lie close together, or are joined by a bridge of cells, are not necessarily related by lineage or development, and well-separated islands can be artefacts of the method. They are tools for visualisation, not for quantitative analysis; clustering and trajectory inference should be done on the underlying data, such as the PCA space, not on the two-dimensional embedding (Chari and Pachter 2023).

7.4.7 Choosing and validating a clustering

Clustering always finds clusters. Hierarchical clustering returns a tree and k-means returns \(k\) groups even for random noise. The result depends on choices made before any algorithm runs, so a clustering is a hypothesis to be checked, not a finding in itself (D’haeseleer 2005).

Distance measure. The distance decides what “similar” means. Euclidean distance between two gene profiles is dominated by the absolute level: two highly expressed genes are close even if one rises when the other falls. Correlation distance, \(d = 1 - r\), compares the shape of the profiles: genes that go up and down together are close, whatever their level. With \(1 - |r|\), anticorrelated genes are also grouped together, which may be wanted for a repressor and its targets but mixes opposite responses. Spearman correlation is less sensitive to outliers. The two families are related: for profiles over \(n\) samples that are centred and scaled to unit variance, the squared Euclidean distance is \(2(n-1)(1-r)\).

Scaling. Clustering is done on log-transformed or variance-stabilised values, so that a few very highly expressed genes do not dominate. When genes are clustered by the shape of their response, each gene is standardised (centred and divided by its standard deviation). When samples are clustered, the genes with the largest variance dominate the distance, which is why the analysis is usually restricted to the most variable genes; standardising every gene instead gives the noise of flat, unexpressed genes the same weight as real signal. Strong batch effects become clusters of their own.

Linkage. In hierarchical clustering, single linkage tends to chain points into long, straggling clusters and is sensitive to noise; complete linkage gives compact clusters of similar diameter; average linkage (UPGMA) lies between them and is common for expression data; Ward’s method joins the pair of clusters that least increases the within-cluster variance and, like k-means, favours compact, spherical clusters. A dendrogram becomes clusters only when it is cut at some height, and the left–right order of its leaves is partly arbitrary, because each node can be flipped. Neighbouring rows in a clustered heat map are therefore not necessarily similar.

How many clusters? There is rarely one correct number. Internal criteria compare the tightness of clusters with their separation. The elbow in the within-cluster sum of squares as a function of \(k\) is often unclear. The gap statistic compares this sum with that of reference data without cluster structure (Tibshirani et al. 2001). The silhouette of a point is \(s = (b - a)/\max(a, b)\), where \(a\) is its mean distance to the other members of its own cluster and \(b\) its mean distance to the members of the nearest other cluster; \(s\) lies between \(-1\) and \(1\), and an average close to zero means that the data have little cluster structure (Rousseeuw 1987). Model-based clustering with Gaussian mixtures can choose \(k\) with a criterion such as BIC.

Validation. Stability is tested by resampling: the samples or genes are subsampled or bootstrapped, the data re-clustered, and pairs of items that are grouped together consistently indicate robust clusters. Consensus clustering summarises this as a matrix of co-clustering frequencies (Monti et al. 2003). Stable clusters can still be artefacts, for example of batch. External validation asks whether the clusters mean something biologically: whether gene clusters are enriched in functional annotation such as GO terms (Gibbons and Roth 2002) (but see the caveats in Section 7.5), or whether sample clusters agree with known labels, measured for example with the adjusted Rand index. The strongest test is that the same clusters appear in an independent data set.

Biclustering. Ordinary clustering groups genes by their behaviour across all samples, but many genes are co-regulated only under some conditions or in some tumours. Biclustering searches for submatrices: subsets of genes with a coherent pattern over subsets of samples. It was introduced to expression data by Cheng and Church in 2000, and many algorithms followed (Madeira and Oliveira 2004). Matrix factorisations, such as non-negative matrix factorisation, address the same problem by letting a gene contribute to several overlapping programs.

7.4.8 Classification with many genes and few samples

Clustering is class discovery; predicting a known class, such as tumour type, response to a drug or survival, is class prediction, a supervised learning problem (Section 12.2). The early leukaemia study did both (Golub et al. 1999). Expression data are extreme for supervised learning: the number of features \(p\) (about 20,000 genes) is far larger than the number of samples \(n\) (often tens to a few hundred), the \(p \gg n\) problem. This has several consequences:

• With more genes than samples, a linear classifier can separate any labelling of the training samples perfectly, including random labels. Accuracy on the training data means nothing.
• Among thousands of genes, many correlate with the outcome by chance.
• Many different gene sets predict equally well, so a signature is not unique. Most random sets of genes are significantly associated with breast cancer outcome, because so many genes correlate with proliferation (Venet et al. 2011).

Models for this situation are simple and regularised: nearest shrunken centroids shrinks each class mean towards the overall mean, and genes for which every class mean is shrunk all the way to it drop out of the classifier (Tibshirani et al. 2002); penalised regression (lasso, elastic net) and linear discriminant analysis with a diagonal covariance matrix are also common. How the penalties work is described in Section 12.2.

The most common error is selection bias. If the genes that best separate the classes are chosen using all samples, and a classifier on those genes is then cross-validated, the test samples have already influenced the model, and the estimated error can be close to zero even for random labels (Ambroise and McLachlan 2002; Simon et al. 2003). Gene selection, and the tuning of any penalty, are part of the model and must be repeated inside each cross-validation fold (nested cross-validation). Even a correctly cross-validated estimate comes from the same cohort, laboratory and platform. A clinical signature needs validation in an independent cohort, and must improve on the established clinical predictors, such as tumour grade and size, not only on chance.

See also: further explanation in Section D.7.3; related Wikipedia articles in Section E.7.3.

NoteSources

Parts of this section are adapted from the Wikipedia articles Negative binomial distribution, Multiple comparisons problem, False discovery rate, Hierarchical clustering, K-means clustering, Principal component analysis and T-distributed stochastic neighbor embedding (CC BY-SA 4.0, retrieved 4 October 2026). The statistical background is adapted from Bayes’ theorem, Student’s t-test and p-value. Parts of this section are adapted from Matteo Tadiello, Modern Bioinformatics: From Foundations to Frontiers, chapter “Gene-Expression Analysis” (CC BY 4.0), used with permission.

7.5 From gene lists to biology

A differential expression analysis typically ends with a list of hundreds or thousands of genes. Interpreting such a list gene by gene is impossible, and individual genes are often not significant even when a whole pathway changes modestly. This section describes how gene lists are interpreted with prior knowledge, and how mRNA levels relate to protein levels.

7.5.1 What a change in expression means

A difference in the measured amount of a transcript does not by itself show that its transcription has changed. The amount of an RNA, \(M\), is a balance between synthesis and degradation,

\[ \frac{dM}{dt} = k_\text{syn} - k_\text{deg} M, \]

so that at steady state \(M = k_\text{syn}/k_\text{deg}\). A higher level can result from faster transcription, from changes in processing or export, or from slower degradation, and an ordinary RNA-seq experiment cannot tell these apart. Methods that label newly made RNA, or that compare unspliced and spliced reads (as RNA velocity does; Section 7.6), are needed for that.

A bulk sample is also a mixture of cell types. The measured expression of gene \(g\) is approximately

\[ A_g = \sum_c \pi_c\, a_{gc}, \]

where \(\pi_c\) is the contribution of cell type \(c\) and \(a_{gc}\) the expression of the gene in that cell type. \(A_g\) can change because the expression within a cell type changes, or because the proportions of the cell types change, for example when a disease leads to an influx of immune cells into a tissue. These are different biological conclusions. Bulk data can be deconvolved into estimated proportions of cell types, using reference expression profiles of purified cells or single-cell data; CIBERSORT, for example, estimates the proportions of immune cell types from tissue expression profiles with support vector regression (Newman et al. 2015). Deconvolution depends on the reference being appropriate for the tissue. Single-cell data separate the two kinds of change directly (Section 7.6).

Finally, an expression study is usually exploratory, and its key results should be validated with independent samples or with a different method, such as quantitative PCR, in situ hybridisation, protein measurements or perturbation experiments. Gene set analysis summarises the same data in another way and is not an independent confirmation. Agreement between RNA and protein strengthens an association, but neither shows that the gene causes the phenotype.

7.5.2 Gene set analysis

Gene sets are groups of genes that share a function, a location or a regulatory mechanism: Gene Ontology terms (Section 2.6), pathways from KEGG or Reactome, targets of a transcription factor or microRNA, or signatures from earlier experiments. The Molecular Signatures Database (MSigDB) collects many such sets, including a small curated “hallmark” collection (Liberzon et al. 2015).

Over-representation analysis (ORA) asks whether a gene set is over-represented among the significant genes. With \(N\) genes in total, \(K\) of which belong to the set, and \(n\) significant genes, of which \(k\) are in the set, the probability of seeing \(k\) or more by chance is given by the hypergeometric distribution (equivalent to a one-sided Fisher’s exact test):

\[ P(X \ge k) = \sum_{i \ge k} \frac{\binom{K}{i}\binom{N-K}{n-i}}{\binom{N}{n}}. \]

Tools such as DAVID made this routine for many annotation types (Huang et al. 2008). ORA is simple, but it ignores the size of the changes and assumes that genes are independent. Two choices are part of the null hypothesis itself, not technical details. The first is the background (the universe of \(N\) genes): it should be the genes that were actually measured and could have been detected, not the whole genome. With the whole genome as background, any set of genes expressed in the tissue studied looks enriched, simply because unexpressed genes could never have been selected. The second is the threshold that defines the gene list: changing the FDR or fold-change cutoff changes \(n\) and \(k\), and can change which sets are significant. Results that hold only for one background or one threshold should be treated with caution. Since thousands of gene sets are tested, multiple-testing correction is needed (Section 7.4), and because gene sets overlap heavily, the significant sets are often redundant.

Gene set enrichment analysis (GSEA) avoids the cutoff (Subramanian et al. 2005). All genes are ranked by a statistic of differential expression. For each gene set, a running-sum statistic walks down the ranked list, increasing when it meets a gene in the set and decreasing otherwise. The maximum deviation from zero, the enrichment score, is large when the members of the set are concentrated at the top or bottom of the list, similar to a Kolmogorov–Smirnov statistic. Significance is assessed by permuting the sample labels, which keeps the correlation between genes, or the gene labels, which is faster but ignores it. GSEA can detect pathways in which many genes change modestly, the typical situation in biology.

Methods differ in the null hypothesis they test. Competitive tests ask whether the genes in a set are more differentially expressed than other genes. Self-contained tests ask whether any of the genes in the set are differentially expressed at all. Pathway methods that also use the topology of the network, and network-based approaches, are discussed in Chapter 11.

7.5.3 Co-expression and regulatory variation

Expression data also reveal relationships between genes. Genes with correlated expression across many conditions often share functions, which is used to predict function (Chapter 10) and to define modules (Langfelder and Horvath 2008). Combined with genotypes, expression becomes a quantitative trait: expression quantitative trait loci (eQTLs) are genetic variants associated with the expression of a gene. The GTEx project mapped eQTLs in about 50 tissues from nearly a thousand donors and found that almost all genes have genetic variants that affect their expression (The GTEx Consortium 2020). eQTLs are an important tool for interpreting the non-coding variants found in genome-wide association studies.

7.5.4 Protein-level expression

What matters for most cellular functions is the amount of protein, not mRNA. How well do mRNA levels predict protein levels? Early studies found correlations that were surprisingly low. A large study in mouse cells measured mRNA and protein levels and turnover for thousands of genes, and concluded that protein levels are determined mainly by translation rates (Schwanhäusser et al. 2011). A later reanalysis argued that much of the apparent discrepancy is due to measurement noise, and that across genes in steady state, mRNA levels explain most of the variation in protein levels, while during dynamic changes, and for comparisons of the same gene between conditions, post-transcriptional regulation is more important (Liu et al. 2016).

The Human Protein Atlas (HPA), a Swedish project started in 2003, combines RNA-seq of human tissues and cell lines with antibody-based imaging of the corresponding proteins, and classifies each gene by where it is expressed, from ubiquitous housekeeping genes to genes enriched in a single tissue (Uhlén et al. 2015). Mass-spectrometry and antibody-based measurement of proteins are described in Chapter 8, and the HPA in Section 2.8.

See also: further explanation in Section D.7.4; related Wikipedia articles in Section E.7.4.

NoteSources

Parts of this section are adapted from the Wikipedia articles Gene set enrichment analysis, Hypergeometric distribution, Expression quantitative trait loci and Human Protein Atlas (CC BY-SA 4.0, retrieved 4 October 2026). Parts of this section are adapted from Matteo Tadiello, Modern Bioinformatics: From Foundations to Frontiers, chapter “Gene-Expression Analysis” (CC BY 4.0), used with permission.

7.6 Single-cell genomics

Bulk RNA-seq measures the average expression of millions of cells. A tissue is a mixture of cell types, and even cells of the same type differ, so the average can hide the most interesting biology. Single-cell RNA-seq (scRNA-seq) measures the transcriptome of each cell separately. It has become one of the fastest-developing areas of genomics, with its own experimental methods, statistical models and analysis workflows (Luecken and Theis 2019; Heumos et al. 2023).

7.6.1 Technologies

The first single-cell transcriptome was published in 2009 (Tang et al. 2009). The number of cells that can be profiled in one experiment has since grown exponentially, from a handful to millions (Svensson et al. 2018):

• Plate-based methods, such as Smart-seq2 (Picelli et al. 2013), sort single cells into wells and produce full-length cDNA. They give high sensitivity and full transcript coverage for hundreds to thousands of cells.
• Droplet-based methods, such as Drop-seq (Macosko et al. 2015), inDrop (Klein et al. 2015) and the commercial 10x Genomics Chromium platform (Zheng et al. 2017), encapsulate single cells in nanolitre droplets together with beads carrying barcoded primers. Every transcript from the same cell gets the same cell barcode, so tens of thousands of cells can be sequenced together and separated computationally.
• Combinatorial indexing methods label cells in several rounds of split-pool barcoding without physically isolating them, and reach millions of cells.

Most high-throughput methods sequence only the 3’ or 5’ end of each transcript, and tag each original molecule with a unique molecular identifier (UMI), a random sequence added before amplification (Islam et al. 2013). Reads with the same UMI come from the same molecule, so counting UMIs instead of reads removes amplification bias.

7.6.2 Properties of the data

A typical cell contains a few hundred thousand mRNA molecules, of which only a few thousand to tens of thousands are captured. The data are therefore a very sparse matrix of counts, with most entries zero. Many zeros were first thought to be technical dropouts requiring special zero-inflated models, but for UMI data the number of zeros is largely what is expected from sampling a small number of molecules at the given sequencing depth, and negative binomial or Poisson models fit well (Svensson 2020). Whether extra zeros need to be modelled therefore depends on the protocol and the model; it may matter more for read counts from full-length, non-UMI protocols. Many observed zeros also reflect real biological differences between cell types. Other technical problems are specific to single cells: empty droplets, doublets (two cells with the same barcode), ambient RNA from lysed cells, and dying cells with high fractions of mitochondrial reads.

7.6.3 The analysis workflow

The standard workflow, implemented in the R package Seurat (Hao et al. 2021) and the Python package Scanpy (Wolf et al. 2018), consists of (Luecken and Theis 2019; Heumos et al. 2023):

1. Quality control: removing barcodes with too few genes, too many counts (possible doublets) or a high mitochondrial fraction.
2. Normalisation for the number of molecules captured per cell, for example by scaling and log-transforming, by pooling cells to estimate size factors (Lun et al. 2016), or with regularised negative binomial regression (Hafemeister and Satija 2019).
3. Feature selection: keeping the most variable genes.
4. Dimensionality reduction with PCA (Section 7.4).
5. Building a k-nearest-neighbour graph of cells in PCA space and clustering it with community-detection algorithms such as Louvain or Leiden (Traag et al. 2019).
6. Visualisation with UMAP or t-SNE.
7. Finding marker genes for each cluster and annotating the clusters as cell types, by comparison with known markers or by transferring labels from annotated reference atlases.

A cluster is a computational result, not automatically a cell type. The number of clusters depends on the resolution parameter, on the choice of genes and on batch correction; a continuous process such as differentiation can be cut into several arbitrary clusters; and residual batch effects, doublets or stressed cells can form clusters of their own. Annotation is therefore an inference that should be supported by several known markers, by biological context and, where possible, by independent validation, for example with imaging or protein markers.

Further analyses depend on the question:

• Data integration. Samples processed separately differ by batch effects, which must be removed before cells from different samples, individuals or technologies can be compared, without removing real biological differences. Harmony iteratively corrects the PCA embedding (Korsunsky et al. 2019). scVI uses a variational autoencoder, a deep generative model of the counts with batch as a covariate (Lopez et al. 2018). Integration is a trade-off between mixing the batches and keeping the biological variation, and a large benchmark found that the balance depends on the method and the preprocessing, in particular when biology and batch overlap (Luecken et al. 2021). If a disease state is present in only one batch, forcing the batches to mix can remove the disease signal. An integrated embedding is useful for clustering, visualisation and transfer of labels, but tests between conditions should go back to the uncorrected counts and a design that respects donors (Figure 7.3).
• Trajectory inference. Cells in development or differentiation form continua rather than discrete clusters. Pseudotime methods order cells along a path through expression space (Trapnell et al. 2014). RNA velocity uses the ratio of unspliced to spliced reads to estimate whether each gene is being switched on or off, and so predicts the future state of each cell (La Manno et al. 2018).
• Differential expression between conditions should be tested between samples, not between individual cells. Cells from the same individual share its genotype, environment and sample processing, so they are not independent biological replicates; the experimental unit is the individual (Section 7.3). Treating thousands of cells as replicates is pseudoreplication and gives very many false positives (Squair et al. 2021; Zimmerman et al. 2021). In pseudobulk analysis, the counts of all cells of one cell type are summed for each sample, and the resulting sample-level counts are analysed with bulk methods such as edgeR, DESeq2 or limma. Mixed models with the individual as a random effect are an alternative; a comparison of the two strategies found that pseudobulk methods performed well (Crowell et al. 2020). More cells improve the description of each individual, but only more individuals support conclusions about a population.
• Differential abundance asks a different question: whether the proportions of cell types or states differ between conditions, for example an expansion of a T-cell population in disease, rather than whether genes change within a cell type. Proportions are compositional, since an increase in one cell type must lower the share of the others, and they must also be compared between samples, not between cells. Methods test either the cluster proportions or local neighbourhoods in the kNN graph, which avoids depending on a fixed clustering (Dann et al. 2021).

Figure 7.3: One single-cell count matrix supports two different analyses. Exploratory analyses of an integrated representation are used for clustering, visualisation and label transfer; tests between conditions aggregate the uncorrected counts by donor and cell type, or use mixed models with the donor as a random effect. Figure: Matteo Tadiello, Modern Bioinformatics: From Foundations to Frontiers (CC BY 4.0).

7.6.4 Multimodal single-cell data

RNA is only one layer of a cell’s state. Multimodal (multiome) methods measure more than one layer in the same cell, using the shared cell barcode to link them. CITE-seq adds oligonucleotide-tagged antibodies, so that surface proteins are counted together with the transcriptome (Stoeckius et al. 2017). Joint RNA and ATAC methods, such as the 10x Multiome kit, measure gene expression and chromatin accessibility (Section 6.9) in the same nucleus, which links candidate regulatory regions to the genes whose expression they correlate with. Because each modality has its own noise and sparsity, they are first analysed separately and then combined, for example with Seurat’s weighted nearest neighbour method, which learns for each cell how much to trust each modality (Hao et al. 2021). When modalities are measured in different cells, they can still be integrated by mapping both into a common space, but the link between them is then inferred rather than observed.

7.6.5 Spatial transcriptomics

Dissociating a tissue into single cells loses the information about where each cell was. Spatial transcriptomics measures expression while keeping its position, and was named Method of the Year 2020 by Nature Methods (Marx 2021). There are two families of methods:

• Sequencing-based methods capture RNA from a tissue section on an array of spatially barcoded spots. The method was developed in Stockholm by Patrik Ståhl, Joakim Lundeberg and colleagues (Ståhl et al. 2016) and commercialised as 10x Visium. Early spots were 100 µm across and contained several cells. Later methods, such as Slide-seq, reach near-cellular resolution (Rodriques et al. 2019).
• Imaging-based methods, such as MERFISH, detect individual RNA molecules in place with multiple rounds of fluorescent hybridisation, and measure hundreds to thousands of chosen genes at subcellular resolution (Chen et al. 2015).

Analysis tasks include identifying spatially variable genes and tissue domains, deconvolution of spots into mixtures of cell types using single-cell references, and the study of cell–cell communication between neighbouring cells. Spatial data bring their own statistical problems. Neighbouring spots are correlated (spatial autocorrelation), so methods that treat locations as independent overstate the evidence; tissue edges, histology and local cell density can mimic spatial patterns; and a single section is one sample, so conclusions about a tissue need replicate sections, ideally from several individuals. Platforms also differ in sensitivity as well as resolution, and the nominal spot size alone is not a sufficient basis for choosing between them (You et al. 2024). Spatial and multimodal methods are converging: spatial measurements of proteins and chromatin accessibility are being added to spatial RNA, and dissociated single-cell multiome data are mapped onto tissue sections, so that each cell’s type, regulatory state and position can be studied together.

7.6.6 Foundation models for single-cell data

Following the success of language models for proteins (Section 1.10), foundation models have been trained on tens of millions of single-cell transcriptomes. Geneformer represents each cell as a list of its genes ranked by expression, and learns with masked-gene prediction (Theodoris et al. 2023). scGPT uses a generative transformer (Cui et al. 2024). The models produce embeddings of genes and cells that can be fine-tuned for cell-type annotation, batch integration or the prediction of the effects of perturbations. Their value is still debated: independent evaluations found that, without fine-tuning, their embeddings often do not outperform simple baselines such as highly variable genes and PCA (Kedzierska et al. 2025), and that for predicting the effects of genetic perturbations they do not yet beat simple linear models (Ahlmann-Eltze et al. 2025). This is an area where students should expect rapid change, and where critical evaluation against simple baselines is essential (Chapter 12).

7.6.7 Cell types and cell atlases

Single-cell data have revived an old question: what is a cell type? Classically, cell types were defined by morphology, location, function or a few marker proteins. Single-cell transcriptomics proposes them as clusters in expression space, which raises new questions: at which resolution should clusters be cut, how do discrete types relate to continuous states, and how should types defined by expression be matched to those defined by other properties (Zeng 2022)? The Human Cell Atlas aims to map all cell types in the human body (Regev et al. 2017), and atlases of whole organisms, such as Tabula Muris (The Tabula Muris Consortium 2018) and Tabula Sapiens (The Tabula Sapiens Consortium 2022), provide references across organs. Such references, together with ontologies of cell types, make it possible to annotate new data sets automatically. The atlases, and the efforts to build models of cells from them, are discussed further in Section 13.7.

See also: further explanation in Section D.7.5; related Wikipedia articles in Section E.7.5.

NoteSources

Parts of this section are adapted from the Wikipedia articles Single-cell transcriptomics, Single-cell sequencing, Unique molecular identifier, Spatial transcriptomics, RNA velocity, Foundation model and Human Cell Atlas (CC BY-SA 4.0, retrieved 4 October 2026). Parts of this section are adapted from Matteo Tadiello, Modern Bioinformatics: From Foundations to Frontiers, chapter “Gene-Expression Analysis” (CC BY 4.0), used with permission.

7.7 Single-cell sequencing and cancer evolution

Cancer is an evolutionary process within the body. Single-cell methods have made it possible to follow this process cell by cell, and to apply to tumours the phylogenetic methods developed for species (Chapter 4).

7.7.1 Cancer as clonal evolution

In 1976 Peter Nowell proposed that tumours evolve like populations of organisms: they start from a single cell, and successive mutations followed by selection of the fittest subclones lead to increasingly malignant cell populations (Nowell 1976). The concepts of Chapter 3 apply directly: mutation, selection, drift and population size all shape a tumour, and the result is intratumour heterogeneity, genetically different subclones within the same tumour. This heterogeneity matters clinically, because a treatment that kills the dominant clone selects for resistant subclones, which then cause relapse.

Bulk sequencing of tumours measures the frequency of each mutation in a mixture of normal and tumour cells. Clonal mutations, present in all tumour cells, occurred early, while subclonal mutations, present in only some, occurred later. Sequencing several regions of the same tumour showed that a biopsy from one region misses many of the mutations, and that tumours evolve in a branched rather than a linear way (Gerlinger et al. 2012; McGranahan and Swanton 2017). Reconstructing the subclones and their phylogeny from bulk data is a deconvolution problem with many possible solutions (Schwartz and Schäffer 2017).

7.7.2 Single-cell DNA sequencing

Single-cell sequencing removes the need for deconvolution. A cell contains only two copies of its genome, about 6 pg of DNA, which must be amplified before sequencing, and whole-genome amplification introduces errors, uneven coverage and the loss of one allele (allelic dropout) (Gawad et al. 2016). Copy-number changes, which affect large regions, can nevertheless be measured robustly even with very shallow sequencing of each cell. Navin and colleagues used single-nucleus sequencing of breast tumours to show that copy-number profiles grouped cells into a few major clones, and that a tumour can arise through a short burst of large genomic changes followed by stable clonal expansion (punctuated evolution) (Navin et al. 2011). High-throughput single-cell DNA sequencing of thousands of cells later showed that breast tumours maintain a reservoir of subclonal diversity during their growth (Minussi et al. 2021).

Trees of single cells are built with methods adapted from phylogenetics (Section 4.3), with models that account for the specific errors of single-cell data, such as allelic dropout and false-positive mutations, and for the fact that copy-number changes are not independent events at single sites (Schwartz and Schäffer 2017). Many mutation-tree methods rest on the infinite-sites assumption: each site mutates at most once in the history of the tumour, and a mutation, once gained, is never lost. Cancer genomes violate this. A deletion or a loss of heterozygosity can remove a mutation from a subclone, and copy-number changes alter the number of copies of a mutated allele, which also distorts the variant allele frequencies used in bulk deconvolution. In addition, a single biopsy or a small sample of cells may miss subclones altogether. Trees inferred under the infinite-sites model should therefore be read as hypotheses, and newer methods relax the assumption, for example by allowing mutation losses or by modelling single-nucleotide and copy-number changes jointly.

7.7.3 Single-cell transcriptomics of tumours

Single-cell RNA-seq of tumours measures both the cancer cells and the tumour microenvironment: immune cells, fibroblasts and blood vessels (Lawson et al. 2018). In glioblastoma, cells from the same tumour expressed programmes characteristic of different molecular subtypes, which had been defined from bulk data as properties of whole tumours (Patel et al. 2014). In metastatic melanoma, scRNA-seq separated malignant from non-malignant cells, identified cancer cell states associated with drug resistance and characterised the state of the infiltrating T cells (Tirosh et al. 2016).

Cancer cells are usually identified in scRNA-seq data by inferring large copy-number changes from the expression of neighbouring genes along each chromosome, since an amplified region raises the average expression of all its genes (Patel et al. 2014; Gao et al. 2021). This also makes it possible to group cancer cells into genetic subclones and link genotype to transcriptional state.

Single-cell studies have shown that not all heterogeneity is genetic. Cancer cells with the same genome can occupy different transcriptional or epigenetic states, and switch between them. Some of these states, such as slowly dividing persister cells, survive treatment without any resistance mutation. Integrating genetic and non-genetic heterogeneity, using multi-omic methods that measure DNA, RNA and chromatin in the same cell, and lineage tracing, is a major current goal (Nam et al. 2020).

See also: further explanation in Section D.7.6; related Wikipedia articles in Section E.7.6.

NoteSources

Parts of this section are adapted from the Wikipedia articles Somatic evolution in cancer, Tumour heterogeneity and Single-cell sequencing (CC BY-SA 4.0, retrieved 4 October 2026).

7.8 Key concepts

A checklist for revision.

• Before microarrays: Northern blots, ESTs and SAGE
• Microarrays: spotted cDNA (two-colour) vs. Affymetrix oligonucleotide (one-colour) arrays; Illumina BeadArrays and SNP genotyping arrays; background correction, log transformation, MA plots and loess, quantile normalisation, RMA; MIAME and reproducibility; limitations (known sequences, cross-hybridisation, non-linear intensity response and dynamic range)
• Landmark microarray studies: co-expression clusters; cancer classification and breast cancer subtypes; expression signatures in the clinic; tissue atlases
• RNA-seq experiments: poly(A) selection vs. rRNA depletion; strand-specific libraries; sequencing depth; the experimental unit; biological vs. technical replicates, pseudoreplication and how many replicates are needed; batch effects and confounding; randomisation, blocking and paired designs
• RNA-seq quality control: RNA-specific statistics (unique mapping, rRNA and mitochondrial fractions, strandedness, gene-body coverage and 3′ bias); sample-level PCA and identity checks; outliers as a reason to investigate; the SEQC/MAQC-III reproducibility study
• RNA-seq analysis: spliced alignment (TopHat, STAR, HISAT2); reference-guided and de novo transcript assembly (Cufflinks, StringTie, Trinity); quantification with EM (RSEM); pseudo-alignment and equivalence classes (kallisto, Salmon); long-read RNA-seq for full-length isoforms and transcript discovery
• Normalisation: CPM, RPKM/FPKM, TPM and what TPM does not correct; compositional data and composition bias; TMM and median-of-ratios size factors; global shifts and spike-ins; raw counts as input to DE tools
• Differential expression: why a plain t-test fails; empirical Bayes and moderated statistics (limma, voom); Poisson vs. negative binomial models; dispersion and its shrinkage; GLMs and design matrices; fold-change shrinkage (DESeq2, edgeR); Wald, likelihood-ratio and quasi-likelihood F-tests
• Statistical background: Bayes’ theorem and base rates; binomial, Poisson, normal, negative binomial and extreme value distributions; likelihood; p-values and what they are not; effect sizes and confidence intervals; t-tests, non-parametric, χ²/Fisher and permutation tests; correlation; type I/II errors and power
• Design and contrasts: design matrices and full rank; paired designs (~ donor + treatment); contrasts; interactions as a difference of differences; why “significant in one group but not the other” is not a test of a difference
• Differential transcript usage: isoform proportions vs. total gene expression; uncertainty of transcript estimates; summarising transcripts to genes
• Multiple testing: FWER and Bonferroni; FDR and the Benjamini–Hochberg procedure; q-values and \(\pi_0\); local FDR; why FDR is not a per-gene probability; p-value histograms; significance vs. effect size; volcano plots and fold-change thresholds; independent filtering
• Clustering and dimensionality reduction: distance measures; hierarchical clustering and linkage; k-means; co-expression networks; PCA; t-SNE and UMAP and their pitfalls (distances and proximity are not lineage)
• Choosing and validating clusters: clustering always finds clusters; Euclidean vs. correlation distance and what each groups together; log transformation, standardising genes and selecting variable genes; single, complete, average and Ward linkage; cutting a dendrogram; choosing \(k\) (elbow, gap statistic, silhouette, BIC); stability by resampling and consensus clustering; external validation with annotation and independent data; biclustering
• Classification from expression data: class discovery vs. class prediction; the \(p \gg n\) problem and perfect separation of random labels; unstable signatures and random signatures; nearest shrunken centroids and penalised regression; selection bias when genes are chosen outside cross-validation; independent validation cohorts and clinical baselines
• What an expression change means: synthesis and degradation at steady state; bulk samples as mixtures of cell types; changes in composition vs. changes within cells; deconvolution (CIBERSORT); validation with independent samples and orthogonal methods
• Interpretation: gene sets and MSigDB; over-representation analysis and the hypergeometric test; choice of background and threshold; GSEA; competitive vs. self-contained tests; eQTLs and GTEx
• mRNA and protein: how well mRNA predicts protein levels; the Human Protein Atlas
• Single-cell RNA-seq: plate-based vs. droplet methods; cell barcodes and UMIs; sparsity and whether zeros are inflated; doublets, ambient RNA and QC; normalisation, HVGs, PCA, kNN graphs, Leiden clustering, UMAP; marker genes and annotation, and why clusters are not automatically cell types; integration (Harmony, scVI) and the trade-off between batch mixing and conserved biology; pseudotime and RNA velocity; pseudobulk differential expression and why cells are not replicates; differential abundance of cell types; multimodal data (CITE-seq, joint RNA+ATAC, weighted nearest neighbours)
• Spatial transcriptomics: sequencing-based (spatial barcodes, Visium, Slide-seq) vs. imaging-based (MERFISH); deconvolution; spatial autocorrelation and replicate sections
• Foundation models: Geneformer and scGPT; evaluations against simple baselines
• Cell types: what is a cell type; the Human Cell Atlas, Tabula Muris and Tabula Sapiens
• Cancer evolution: Nowell’s clonal evolution; intratumour heterogeneity; clonal vs. subclonal mutations; branched and punctuated evolution; single-cell DNA sequencing and its errors; the infinite-sites assumption and how copy-number changes and loss of heterozygosity violate it; inferring copy number from scRNA-seq; tumour microenvironment; non-genetic heterogeneity and persister cells

7.9 Recommended reading

These papers are recommended for reading in full:

• Wang et al. (2009) — RNA-seq, a revolutionary tool for transcriptomics; why sequencing replaced arrays.
• Conesa et al. (2016) — a survey of best practices for RNA-seq data analysis, from design to interpretation.
• Love et al. (2014) — DESeq2; negative binomial models, dispersion and fold-change shrinkage.
• Benjamini and Hochberg (1995) — the false discovery rate; the most important idea in genome-scale statistics.
• Luecken and Theis (2019) — best practices in single-cell RNA-seq analysis, step by step.
• Nam et al. (2020) — cancer evolution studied with single-cell multi-omics.

Titles, abstracts, links and the full texts are in Section A.7; a longer list of further papers is in Section B.7.

7.10 Background reading (open access)

• DNA microarray (Wikipedia) — A DNA microarray is a collection of microscopic DNA spots attached to a solid surface. Scientists use DNA microarrays to measure the expression levels of large numbers of genes simultaneously or to genotype multiple regions of a genome. Each DNA spot contains picomoles of a specific DNA sequence, known as probes.
  https://en.wikipedia.org/wiki/DNA_microarray
• RNA-Seq (Wikipedia) — RNA-Seq is a next-generation sequencing (NGS) technique used to quantify and identify RNA molecules in a biological sample, providing a snapshot of the transcriptome at a specific time. It enables transcriptome-wide analysis by sequencing cDNA derived from RNA.
  https://en.wikipedia.org/wiki/RNA-Seq
• Transcriptomics technologies (Wikipedia) — Transcriptomics technologies are the techniques used to study an organism’s transcriptome, the sum of all of its RNA transcripts. The information content of an organism is recorded in the DNA of its genome and expressed through transcription.
  https://en.wikipedia.org/wiki/Transcriptomics_technologies
• Multiple comparisons problem (Wikipedia) — Multiple comparisons, multiplicity or multiple testing problem occurs when many statistical tests are performed on the same dataset. Each test has its own chance of a Type I error, so the overall probability of making at least one false positive increases as the number of tests grows.
  https://en.wikipedia.org/wiki/Multiple_comparisons_problem
• Bonferroni correction (Wikipedia) — Bonferroni correction is a method to counteract the multiple comparisons problem in statistics. It is named after the mathematician Carlo Emilio Bonferroni.
  https://en.wikipedia.org/wiki/Bonferroni_correction
• False discovery rate (Wikipedia) — In statistics, the false discovery rate (FDR) is a method of conceptualizing the rate of type I errors in null hypothesis testing when conducting multiple comparisons. FDR-controlling procedures are designed to control the FDR, which is the expected proportion of “discoveries” that are false.
  https://en.wikipedia.org/wiki/False_discovery_rate
• Bayes’ theorem (Wikipedia) — In probability theory, Bayes’ theorem, named after Thomas Bayes, gives a mathematical rule for inverting conditional probabilities, allowing the probability of a cause to be found given its effect.
  https://en.wikipedia.org/wiki/Bayes%27_theorem
• Student’s t-test (Wikipedia) — Student’s t-test is a statistical test used to test whether the difference between the response of two groups is statistically significant or not. It is any statistical hypothesis test in which the test statistic follows a Student’s t-distribution under the null hypothesis.
  https://en.wikipedia.org/wiki/Student%27s_t-test
• p-value (Wikipedia) — In null-hypothesis significance testing, the p-value is the probability of obtaining test results at least as extreme as the result actually observed, under the assumption that the null hypothesis is correct.
  https://en.wikipedia.org/wiki/P-value
• Statistics for Biologists — Nature’s collection of short articles for biologists, including the Points of Significance column from Nature Methods, explaining p-values, power, multiple testing, regression, PCA and other statistical concepts.
  https://www.nature.com/collections/qghhqm
• Think Stats — Allen Downey’s free book introducing probability and statistics through Python programming and exploratory data analysis.
  https://greenteapress.com/wp/think-stats-3e/
• Principal component analysis (Wikipedia) — Principal component analysis (PCA) is a linear dimensionality reduction technique with applications in exploratory data analysis, visualization and data preprocessing.
  https://en.wikipedia.org/wiki/Principal_component_analysis
• Gene set enrichment analysis (Wikipedia) — Gene set enrichment analysis (GSEA) (also called functional enrichment analysis or pathway enrichment analysis) is a method to identify classes of genes or proteins that are over-represented in a large set of genes or proteins, and may have an association with different phenotypes (e.g. different organism growth patterns or diseases).
  https://en.wikipedia.org/wiki/Gene_set_enrichment_analysis
• Single-cell transcriptomics (Wikipedia) — Single-cell transcriptomics refers to the quantification and analysis of the transcriptomes of individual cells. Single-cell transcriptomics makes it possible to unravel heterogeneous cell populations, reconstruct cellular developmental pathways, and model transcriptional dynamics.
  https://en.wikipedia.org/wiki/Single-cell_transcriptomics
• Spatial transcriptomics (Wikipedia) — Spatial transcriptomics, or spatially resolved transcriptomics, is a method that captures positional context of transcriptional activity within intact tissue.
  https://en.wikipedia.org/wiki/Spatial_transcriptomics
• Somatic evolution in cancer (Wikipedia) — Somatic evolution is the accumulation of mutations and epimutations in somatic cells during a lifetime, and the effects of those mutations and epimutations on the fitness of those cells. This evolutionary process has first been shown by the studies of Bert Vogelstein in colon cancer.
  https://en.wikipedia.org/wiki/Somatic_evolution_in_cancer

Descriptions of Wikipedia articles are taken from the introductions of the articles (CC BY-SA 4.0).

7.11 Student recommendations

Recommended by students preparing for the exam: papers, reviews, online lectures, tutorials and short explanations. To suggest something, use the suggestion form (SU login).

No recommendations yet.

7.12 Review problems

These problems are adapted from Matteo Tadiello, Modern Bioinformatics: From Foundations to Frontiers (CC BY 4.0), used with permission. For a longer treatment of this topic, see the chapters “Gene-Expression Analysis” and “Statistics and Machine Learning for Bioinformatics” of that book. In the HTML version, click a box to show the extensions, solution outline and common errors.

7.12.1 Designing a bulk RNA-seq experiment

A study aims to identify genes affected by a drug in cultures derived from patients. Formulate an experimental and statistical design that supports a conclusion about the patient population.

TipExtensions, solution outline and common errors

Extensions:

• What is a biological replicate, and is a technical replicate useful?
• Would you prefer more samples or more sequencing depth?
• How do you handle paired control and treated cultures from each patient?
• What if all controls were prepared in one batch and all treated samples in another?

Solution outline:

1. Define what is to be estimated: the average change associated with the drug within patients, in the target population.
2. Treat patients, not wells or reads, as the independent experimental units; base the sample size on the expected effect, the variability and the number of tests.
3. Pair treatment within patients, randomise and balance the processing, and decide the exclusion rules and the primary contrast in advance.
4. Use a design such as ~ patient + treatment; add nuisance covariates that can be estimated, and an interaction only if the question requires it.
5. Once the coverage is adequate, prefer more biological replicates to deeper sequencing; technical replicates describe the measurement, not the population.
6. State that perfect confounding of treatment and batch cannot be repaired computationally.

Common errors:

• Calling several wells from one culture independent biological replicates.
• Choosing the number of samples from “typical RNA-seq practice” without considering the effect size or the dispersion.
• Promising that batch correction will recover a perfectly confounded treatment effect.
• Letting a PCA plot decide the design after the results are visible.

Connections: experimental design and replication (Section 7.3); power and statistical background (Section 7.4.1); sequencing technologies (Section 5.2).

7.12.2 Alignment versus pseudo-alignment

Decide whether RNA-seq reads should be aligned to a genome or assigned to a transcriptome by pseudo-alignment or lightweight mapping, for a specified biological analysis.

TipExtensions, solution outline and common errors

Extensions:

• What does pseudo-alignment keep, and what does it discard?
• How are reads compatible with several isoforms assigned?
• When is genome alignment necessary?
• What happens if the annotation is incomplete?

Solution outline:

1. Start from the feature of interest: gene abundance, transcript abundance, new splice junctions, fusion genes, allele-specific expression or evidence to inspect in a genome browser.
2. Explain spliced genome alignment as placing reads on genomic coordinates across introns, followed by counting or transcript assembly.
3. Explain equivalence classes of compatible transcripts and the probabilistic assignment of reads with a mixture likelihood and the EM algorithm.
4. Prefer fast transcriptome quantification for the abundance of annotated transcripts, and genome alignment for new features or questions that depend on coordinates.
5. Record the versions of genome, transcriptome and annotation, the strandedness and the treatment of multi-mapping reads.
6. Validate isoforms that are hard to distinguish with junction-specific or targeted evidence, and take the uncertainty of the estimates into account.

Common errors:

• Saying that pseudo-alignment means “no mapping”.
• Assuming that paired-end reads uniquely identify every isoform.
• Ignoring mismatches between the reference genome and the annotation.
• Treating transcript estimates from shared reads as direct observations.

Connections: spliced alignment and quantification (Section 7.3); differential transcript usage (Section 7.4); exact matching and indexes (Section 1.4); read mapping (Section 5.4); gene annotation (Section 5.6).

7.12.3 Normalisation and composition effects

Explain why raw counts cannot be compared directly between libraries, and why testing TPM values does not in general give valid differential expression results that take replicates into account.

TipExtensions, solution outline and common errors

Extensions:

• What do TMM and the median-of-ratios normalisation of DESeq2 estimate?
• What is the difference between correcting for library depth and for transcript length?
• What if nearly all genes truly increase in every cell?
• Are spike-ins a complete solution?

Solution outline:

1. Counts depend on sequencing depth, library composition, feature length and protocol bias, as well as on abundance.
2. Describe TMM and median-of-ratios scaling as robust estimates of the relative depth, under the assumption that most genes do not change.
3. Include the scaling factors in the count model as offsets, and use the raw counts for inference.
4. Describe TPM as a length-adjusted proportion within a sample, useful for describing abundance but not a model of the variance between replicates.
5. Explain that a global change per cell cannot be detected from standard relative libraries without an external reference.
6. Discuss spike-ins and cell counts as additional assumptions rather than complete solutions.

Common errors:

• Equating normalised counts with absolute numbers of molecules per cell.
• Saying that normalisation makes all samples biologically identical.
• Applying batch correction or quantile normalisation before asking what signal it removes.
• Comparing RPKM or TPM values with a t-test, ignoring the dependence of the variance on the mean.

Connections: normalisation (Section 7.3); microarray normalisation (Section 7.2); relative and absolute quantification in proteomics (Section 8.4).

7.12.4 Negative binomial differential expression

Derive and interpret the negative binomial generalised linear model used by DESeq2 or edgeR for differential expression with replicates.

TipExtensions, solution outline and common errors

Extensions:

• Why not use a Poisson model?
• What is dispersion, and why is it shrunk?
• How do Wald, likelihood-ratio and quasi-likelihood tests differ?
• How should fold changes and the FDR be reported?

Solution outline:

1. Write \(Y_{gi} \sim \text{NB}(\mu_{gi}, \alpha_g)\) and \(\operatorname{Var}(Y) = \mu + \alpha\mu^2\).
2. Write \(\log\mu_{gi} = \log s_i + x_i^T\beta_g\), and define the offset, the row of the design matrix, the coefficients and the contrast.
3. Explain that the Poisson variance equals the mean and misses the overdispersion between biological replicates.
4. Explain how noisy dispersion estimates of individual genes are shrunk towards a trend fitted across all genes (empirical Bayes).
5. Use a Wald test for a single coefficient or contrast, a likelihood-ratio test for a question involving several coefficients, or a quasi-likelihood test to account for uncertainty in the dispersion.
6. Report the effect size, its uncertainty, the adjusted p-value and the abundance; shrinkage stabilises noisy estimates but does not create evidence.

Common errors:

• Calling dispersion only technical noise.
• Interpreting a shrunken log fold change as the observed change.
• Defining the FDR as the probability that a particular gene is a false positive.
• Choosing the test that gives the longest gene list.

Connections: differential expression and multiple testing (Section 7.4); linear models (Section 12.2); peptide- and protein-level FDR (Section 8.4).

7.12.5 Contrasts, interactions and batch

A drug effect is statistically significant in wild type but not in a mutant. Does this show that the effect of the drug depends on the genotype? Formulate the appropriate test.

TipExtensions, solution outline and common errors

Extensions:

• Which contrast tests whether the effect differs?
• What if genotype and batch are correlated?
• What makes a design matrix rank-deficient?
• How would you analyse an experiment with several time points?

Solution outline:

1. Reject the comparison of two significance labels; test a difference of differences.
2. Specify \((\text{mut}_\text{drug} - \text{mut}_\text{control}) - (\text{WT}_\text{drug} - \text{WT}_\text{control})\) in a model with main effects and an interaction.
3. Include batch when it can be estimated; inspect the design matrix, and explain that perfect confounding makes separate coefficients impossible.
4. For time courses, test the condition-by-time terms jointly, for example with a likelihood-ratio test.
5. Keep the independent experimental unit and any repeated measurements explicit.
6. Interpret an interaction as a context-dependent association, not automatically as a mechanism.

Common errors:

• “Significant here, not significant there” as a test of interaction.
• Removing a batch effect until the groups look mixed in a PCA plot.
• Including more coefficients than the design supports.
• Treating repeated time points from the same individual as independent.

Connections: design matrices, contrasts and interactions (Section 7.4); batch effects and confounding (Section 7.3); linear models (Section 12.2).

7.12.6 Single-cell differential expression

A data set contains 10,000 cells from each of three treated and three control patients. What is the sample size for a conclusion about the effect of treatment in the patient population? Formulate a valid analysis of differential expression within a cell type.

TipExtensions, solution outline and common errors

Extensions:

• Why are cells not independent replicates?
• What is pseudobulk analysis?
• Are zeros “dropouts” that require a zero-inflated model?
• How do you separate differential abundance from differential expression within a cell type?

Solution outline:

1. The sample size is six patients; the cells are nested measurements within patients.
2. Sum the UMI counts per patient and cell type (pseudobulk) and use a bulk count model, or use a mixed model with the patient as a random effect.
3. Treat UMI zeros as non-detections; do not assume extra zero inflation without evidence.
4. Distinguish a change of expression within a cell population from a change in the frequency of that population.
5. Check that each patient contributes enough cells, and show the results per patient.
6. Distinguish exploratory marker genes from tests between conditions.

Common errors:

• Reporting a tiny p-value driven by thousands of correlated cells.
• Downsampling cells and claiming that this creates biological independence.
• Confusing more immune cells with activation within immune cells.

Connections: single-cell analysis, pseudobulk and differential abundance (Section 7.6); experimental units (Section 7.3); pitfalls in evaluation (Section 12.2).

7.12.7 Clustering, annotation and integration

After integration of several single-cell RNA-seq batches, the cells appear well mixed in a low-dimensional embedding. Is this enough to show that the integration was successful?

TipExtensions, solution outline and common errors

Extensions:

• What should be preserved, as well as mixed?
• Is a cluster a cell type?
• Can the corrected matrix be used for differential expression?
• What if a disease state occurs in only one batch?

Solution outline:

1. Define the purpose of the integration: visualisation, clustering, mapping to a reference or label transfer.
2. Assess batch mixing and the conservation of known cell types, continuous states, rare populations and marker genes separately.
3. Explain that biology confounded with batch cannot be recovered; forced mixing may erase a disease signal.
4. Annotate clusters with several positive and negative markers, reference data and consistency across samples.
5. Keep the uncorrected counts and use models that take patients into account for tests between conditions.
6. Compare several methods and parameters; an attractive UMAP is not validation.

Common errors:

• Treating maximal batch mixing as the only goal.
• Naming a cluster after one marker gene.
• Reading distances between clusters or cluster sizes literally from a UMAP.

Connections: single-cell analysis and integration (Section 7.6); clustering and dimensionality reduction (Section 7.4); unsupervised learning (Section 12.2).

7.12.8 Spatial transcriptomics and validation

Design an analysis to determine whether an immune-response programme is localised at the boundary of a tumour, and distinguish spatial association from evidence of signalling between cells.

TipExtensions, solution outline and common errors

Extensions:

• Would you choose sequencing-based spatial transcriptomics or an imaging-based method?
• What does deconvolution of spots infer?
• Why is ordinary differential expression between spots risky?
• How would you validate a proposed ligand–receptor interaction?

Solution outline:

1. Define the boundary, the target genes, the spatial scale, and whether the aim is transcriptome-wide discovery or targeted measurement.
2. Choose sequencing-based capture for broad discovery, or imaging for targeted measurement of single molecules, balancing field of view, sensitivity and resolution.
3. Align expression to histology, account for spots that contain several cells, segmentation errors and spatial autocorrelation.
4. Use replicate sections from several patients or animals; neighbouring spots are not independent individuals.
5. Describe deconvolution as inference of mixtures that depends on a reference, not as direct observation of single cells.
6. Validate expression with an independent imaging method, and signalling with protein localisation, perturbation or a functional response; co-expression alone is a hypothesis.

Common errors:

• Equating the nominal spot size with single-cell resolution.
• Treating an inferred cell type as observed ground truth.
• Inferring physical signalling from ligand and receptor transcripts alone.
• Generalising from one striking tissue section.

Connections: spatial transcriptomics (Section 7.6); what an expression change means (Section 7.5); protein interactions and networks (Section 11.3).

7.12.9 Statistical evidence and biological importance

A transcript has an estimated treatment log fold change of 0.35 and a nominal \(p\)-value of 0.003. Interpret these quantities without overstating the evidence.

TipExtensions, solution outline and common errors

Extensions:

Define the null model and tail event. Contrast a frequentist confidence interval with a Bayesian credible interval. Explain how low power, pseudoreplication, and selective reporting change the interpretation. State what additional evidence is required before describing the transcript as biologically important.

Solution outline:

The \(p\)-value is the probability, under the null model and its assumptions, of a test statistic at least as extreme as the observed statistic. It is neither the probability that the null is true nor the probability that the finding arose by chance. The coefficient estimates a conditional contrast on the model scale; its biological magnitude must be evaluated using an uncertainty interval, assay characteristics, covariate structure, and an appropriate multiplicity procedure. A 95% frequentist confidence procedure has 95% repeated-sampling coverage under the model, whereas a Bayesian credible interval has a posterior-probability interpretation conditional on the likelihood and prior. Independent replication and functional evidence are needed for a mechanistic claim.

Common errors:

Interpreting \(p=0.003\) as a 0.3% probability that the null is true; equating statistical significance with biological relevance; treating non-significance as proof of no effect; counting cells or technical replicates as independent organisms.

Connections: Section 7.4.1 (hypothesis tests and p-values); Section 7.4 (differential expression); Section 8.4 (error control in proteomics); Section 12.8 (prediction and causation).

7.12.10 Multiple testing and false-discovery control

A study tests 20,000 genes. Derive the Benjamini–Hochberg procedure and explain precisely what its false-discovery-rate guarantee does and does not mean.

TipExtensions, solution outline and common errors

Extensions:

Contrast FDR with family-wise error rate. Discuss independence, positive dependence, and arbitrary dependence. Explain adjusted \(p\)-values, \(q\)-values, independent filtering, and the need to define a hypothesis family before inspecting results.

Solution outline:

Let \(V\) be false rejections and \(R\) total rejections. FDR is \(E[V/\max(R,1)]\), whereas FWER is \(P(V\ge 1)\). Sort the \(m\) values so that \(p_{(1)}\le\cdots\le p_{(m)}\), find the largest \(k\) satisfying \(p_{(k)}\le kq/m\), and reject through \(k\). The original BH result assumes independent tests, with extensions to specified positive-dependence structures; Benjamini–Yekutieli provides a more conservative option under arbitrary dependence. The guarantee concerns a repeated-study expectation for a pre-specified family, not the probability that a particular reported gene is false.

Common errors:

Declaring that exactly 5% of reported genes are false; applying BH only after favourable tests have been selected; assuming correlated molecular features either invalidate every FDR method or remove the need for correction.

Connections: Section 7.4 (multiple testing); Section 8.4 (target–decoy FDR); Section 10.4 (function annotation).

7.12.11 Tumour subtypes and a prognostic signature

A group profiles 80 breast tumours by RNA-seq. Hierarchical clustering of the tumours gives three clear groups in a heat map. The group then selects the 50 genes most associated with survival, trains a classifier on them and reports a cross-validated accuracy of 95%. Evaluate both claims.

TipExtensions, solution outline and common errors

Extensions:

• Which distance measure, scaling and linkage would you choose, and how do they change the result?
• How would you decide whether there are three groups rather than two or five?
• How should the cross-validation have been done?
• What would convince you that the signature is clinically useful?

Solution outline:

1. Clustering always returns groups; the heat map alone is not evidence of subtypes.
2. State the choices: variance-stabilised data, the most variable genes, correlation or Euclidean distance on standardised genes, and the linkage; check that the groups survive other reasonable choices.
3. Assess the number of clusters with silhouettes or the gap statistic, and stability by resampling tumours and genes (consensus clustering).
4. Check that the groups do not follow batch, processing date or tumour purity, and validate them externally: known subtypes, enrichment of annotation, and replication in an independent cohort.
5. With \(p \gg n\), selecting the 50 genes on all 80 tumours before cross-validation is selection bias; the selection must be repeated inside each fold, with nested cross-validation for any tuning.
6. Validate the final model on an independent cohort and show that it adds to clinical predictors such as grade, size and nodal status; compare with random gene sets of the same size.

Common errors:

• Treating a dendrogram as proof of discrete subtypes.
• Choosing \(k\) because the heat map “looks like” three groups.
• Reporting training or leaky cross-validated accuracy as performance on new patients.
• Interpreting the 50 selected genes as the causal drivers of outcome.

Connections: clustering and its validation, and classification with many genes (Section 7.4); landmark microarray studies (Section 7.2); regularisation, cross-validation and leakage (Section 12.2).

References

“A Comprehensive Assessment of RNA-Seq Accuracy, Reproducibility and Information Content by the Sequencing Quality Control Consortium.” 2014. Nature Biotechnology 32 (9): 903–14. https://doi.org/10.1038/nbt.2957.

Ahlmann-Eltze, Constantin, Wolfgang Huber, and Simon Anders. 2025. “Deep-Learning-Based Gene Perturbation Effect Prediction Does Not yet Outperform Simple Linear Baselines.” Nature Methods 22 (8): 1657–61. https://doi.org/10.1038/s41592-025-02772-6.

Ambroise, Christophe, and Geoffrey J. McLachlan. 2002. “Selection Bias in Gene Extraction on the Basis of Microarray Gene-Expression Data.” Proceedings of the National Academy of Sciences 99 (10): 6562–66. https://doi.org/10.1073/pnas.102102699.

Anders, Simon, and Wolfgang Huber. 2010. “Differential Expression Analysis for Sequence Count Data.” Genome Biology 11 (10). https://doi.org/10.1186/gb-2010-11-10-r106.

Becht, Etienne, Leland McInnes, John Healy, et al. 2018. “Dimensionality Reduction for Visualizing Single-Cell Data Using UMAP.” Nature Biotechnology 37 (1): 38–44. https://doi.org/10.1038/nbt.4314.

Benjamini, Yoav, and Yosef Hochberg. 1995. “Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing.” Journal of the Royal Statistical Society Series B: Statistical Methodology 57 (1): 289–300. https://doi.org/10.1111/j.2517-6161.1995.tb02031.x.

Bolstad, B. M., R. A Irizarry, M. Åstrand, and T. P. Speed. 2003. “A Comparison of Normalization Methods for High Density Oligonucleotide Array Data Based on Variance and Bias.” Bioinformatics 19 (2): 185–93. https://doi.org/10.1093/bioinformatics/19.2.185.

Bourgon, Richard, Robert Gentleman, and Wolfgang Huber. 2010. “Independent Filtering Increases Detection Power for High-Throughput Experiments.” Proceedings of the National Academy of Sciences 107 (21): 9546–51. https://doi.org/10.1073/pnas.0914005107.

Bray, Nicolas L, Harold Pimentel, Páll Melsted, and Lior Pachter. 2016. “Near-Optimal Probabilistic RNA-Seq Quantification.” Nature Biotechnology 34 (5): 525–27. https://doi.org/10.1038/nbt.3519.

Brazma, Alvis, Pascal Hingamp, John Quackenbush, et al. 2001. “Minimum Information about a Microarray Experiment (MIAME)—Toward Standards for Microarray Data.” Nature Genetics 29 (4): 365–71. https://doi.org/10.1038/ng1201-365.

Chari, Tara, and Lior Pachter. 2023. “The Specious Art of Single-Cell Genomics.” PLOS Computational Biology 19 (8): e1011288. https://doi.org/10.1371/journal.pcbi.1011288.

Chen, Kok Hao, Alistair N. Boettiger, Jeffrey R. Moffitt, Siyuan Wang, and Xiaowei Zhuang. 2015. “Spatially Resolved, Highly Multiplexed RNA Profiling in Single Cells.” Science 348 (6233). https://doi.org/10.1126/science.aaa6090.

Conesa, Ana, Pedro Madrigal, Sonia Tarazona, et al. 2016. “A Survey of Best Practices for RNA-Seq Data Analysis.” Genome Biology 17 (1). https://doi.org/10.1186/s13059-016-0881-8.

Crowell, Helena L., Charlotte Soneson, Pierre-Luc Germain, et al. 2020. “Muscat Detects Subpopulation-Specific State Transitions from Multi-Sample Multi-Condition Single-Cell Transcriptomics Data.” Nature Communications 11 (1). https://doi.org/10.1038/s41467-020-19894-4.

Cui, Haotian, Chloe Wang, Hassaan Maan, et al. 2024. “scGPT: Toward Building a Foundation Model for Single-Cell Multi-Omics Using Generative AI.” Nature Methods 21 (8): 1470–80. https://doi.org/10.1038/s41592-024-02201-0.

D’haeseleer, Patrik. 2005. “How Does Gene Expression Clustering Work?” Nature Biotechnology 23 (12): 1499–501. https://doi.org/10.1038/nbt1205-1499.

Dann, Emma, Neil C. Henderson, Sarah A. Teichmann, Michael D. Morgan, and John C. Marioni. 2021. “Differential Abundance Testing on Single-Cell Data Using k-Nearest Neighbor Graphs.” Nature Biotechnology 40 (2): 245–53. https://doi.org/10.1038/s41587-021-01033-z.

Dobin, Alexander, Carrie A. Davis, Felix Schlesinger, et al. 2012. “STAR: Ultrafast Universal RNA-Seq Aligner.” Bioinformatics 29 (1): 15–21. https://doi.org/10.1093/bioinformatics/bts635.

Efron, Bradley, and Robert Tibshirani. 2002. “Empirical Bayes Methods and False Discovery Rates for Microarrays.” Genetic Epidemiology 23 (1): 70–86. https://doi.org/10.1002/gepi.1124.

Eisen, Michael B., Paul T. Spellman, Patrick O. Brown, and David Botstein. 1998. “Cluster Analysis and Display of Genome-Wide Expression Patterns.” Proceedings of the National Academy of Sciences 95 (25): 14863–68. https://doi.org/10.1073/pnas.95.25.14863.

Fodor, Stephen P. A., J. Leighton Read, Michael C. Pirrung, Lubert Stryer, Amy Tsai Lu, and Dennis Solas. 1991. “Light-Directed, Spatially Addressable Parallel Chemical Synthesis.” Science 251 (4995): 767–73. https://doi.org/10.1126/science.1990438.

Gao, Ruli, Shanshan Bai, Ying C. Henderson, et al. 2021. “Delineating Copy Number and Clonal Substructure in Human Tumors from Single-Cell Transcriptomes.” Nature Biotechnology 39 (5): 599–608. https://doi.org/10.1038/s41587-020-00795-2.

Gawad, Charles, Winston Koh, and Stephen R. Quake. 2016. “Single-Cell Genome Sequencing: Current State of the Science.” Nature Reviews Genetics 17 (3): 175–88. https://doi.org/10.1038/nrg.2015.16.

Gelman, Andrew, and Hal Stern. 2006. “The Difference Between ‘Significant’ and ‘Not Significant’ Is Not Itself Statistically Significant.” The American Statistician 60 (4): 328–31. https://doi.org/10.1198/000313006x152649.

Gerlinger, Marco, Andrew J. Rowan, Stuart Horswell, et al. 2012. “Intratumor Heterogeneity and Branched Evolution Revealed by Multiregion Sequencing.” New England Journal of Medicine 366 (10): 883–92. https://doi.org/10.1056/nejmoa1113205.

Gibbons, Francis D., and Frederick P. Roth. 2002. “Judging the Quality of Gene Expression-Based Clustering Methods Using Gene Annotation.” Genome Research 12 (10): 1574–81. https://doi.org/10.1101/gr.397002.

Golub, T. R., D. K. Slonim, P. Tamayo, et al. 1999. “Molecular Classification of Cancer: Class Discovery and Class Prediction by Gene Expression Monitoring.” Science 286 (5439): 531–37. https://doi.org/10.1126/science.286.5439.531.

Grabherr, Manfred G, Brian J Haas, Moran Yassour, et al. 2011. “Full-Length Transcriptome Assembly from RNA-Seq Data Without a Reference Genome.” Nature Biotechnology 29 (7): 644–52. https://doi.org/10.1038/nbt.1883.

Gunderson, Kevin L., Semyon Kruglyak, Michael S. Graige, et al. 2004. “Decoding Randomly Ordered DNA Arrays.” Genome Research 14 (5): 870–77. https://doi.org/10.1101/gr.2255804.

Hafemeister, Christoph, and Rahul Satija. 2019. “Normalization and Variance Stabilization of Single-Cell RNA-Seq Data Using Regularized Negative Binomial Regression.” Genome Biology 20 (1). https://doi.org/10.1186/s13059-019-1874-1.

Hao, Yuhan, Stephanie Hao, Erica Andersen-Nissen, et al. 2021. “Integrated Analysis of Multimodal Single-Cell Data.” Cell 184 (13): 3573–3587.e29. https://doi.org/10.1016/j.cell.2021.04.048.

Heumos, Lukas, Anna C. Schaar, Christopher Lance, et al. 2023. “Best Practices for Single-Cell Analysis Across Modalities.” Nature Reviews Genetics 24 (8): 550–72. https://doi.org/10.1038/s41576-023-00586-w.

Holm, Sture. 1979. “A Simple Sequentially Rejective Multiple Test Procedure.” Scandinavian Journal of Statistics 6 (2): 65–70.

Huang, Da Wei, Brad T Sherman, and Richard A Lempicki. 2008. “Systematic and Integrative Analysis of Large Gene Lists Using DAVID Bioinformatics Resources.” Nature Protocols 4 (1): 44–57. https://doi.org/10.1038/nprot.2008.211.

Ioannidis, John P. A. 2005. “Why Most Published Research Findings Are False.” PLoS Medicine 2 (8): e124. https://doi.org/10.1371/journal.pmed.0020124.

Irizarry, R. A. 2003. “Exploration, Normalization, and Summaries of High Density Oligonucleotide Array Probe Level Data.” Biostatistics 4 (2): 249–64. https://doi.org/10.1093/biostatistics/4.2.249.

Islam, Saiful, Amit Zeisel, Simon Joost, et al. 2013. “Quantitative Single-Cell RNA-Seq with Unique Molecular Identifiers.” Nature Methods 11 (2): 163–66. https://doi.org/10.1038/nmeth.2772.

Kedzierska, Kasia Z., Lorin Crawford, Ava P. Amini, and Alex X. Lu. 2025. “Zero-Shot Evaluation Reveals Limitations of Single-Cell Foundation Models.” Genome Biology 26 (1). https://doi.org/10.1186/s13059-025-03574-x.

Kim, Daehwan, Joseph M. Paggi, Chanhee Park, Christopher Bennett, and Steven L. Salzberg. 2019. “Graph-Based Genome Alignment and Genotyping with HISAT2 and HISAT-Genotype.” Nature Biotechnology 37 (8): 907–15. https://doi.org/10.1038/s41587-019-0201-4.

Klein, Allon M., Linas Mazutis, Ilke Akartuna, et al. 2015. “Droplet Barcoding for Single-Cell Transcriptomics Applied to Embryonic Stem Cells.” Cell 161 (5): 1187–201. https://doi.org/10.1016/j.cell.2015.04.044.

Korsunsky, Ilya, Nghia Millard, Jean Fan, et al. 2019. “Fast, Sensitive and Accurate Integration of Single-Cell Data with Harmony.” Nature Methods 16 (12): 1289–96. https://doi.org/10.1038/s41592-019-0619-0.

Krzywinski, Martin, and Naomi Altman. 2013. “Significance, p Values and t-Tests.” Nature Methods 10 (11): 1041–42. https://doi.org/10.1038/nmeth.2698.

Kuhn, Kenneth, Shawn C. Baker, Eugene Chudin, et al. 2004. “A Novel, High-Performance Random Array Platform for Quantitative Gene Expression Profiling.” Genome Research 14 (11): 2347–56. https://doi.org/10.1101/gr.2739104.

La Manno, Gioele, Ruslan Soldatov, Amit Zeisel, et al. 2018. “RNA Velocity of Single Cells.” Nature 560 (7719): 494–98. https://doi.org/10.1038/s41586-018-0414-6.

Langfelder, Peter, and Steve Horvath. 2008. “WGCNA: An r Package for Weighted Correlation Network Analysis.” BMC Bioinformatics 9 (1). https://doi.org/10.1186/1471-2105-9-559.

Law, Charity W, Yunshun Chen, Wei Shi, and Gordon K Smyth. 2014. “Voom: Precision Weights Unlock Linear Model Analysis Tools for RNA-Seq Read Counts.” Genome Biology 15 (2). https://doi.org/10.1186/gb-2014-15-2-r29.

Lawson, Devon A., Kai Kessenbrock, Ryan T. Davis, Nicholas Pervolarakis, and Zena Werb. 2018. “Tumour Heterogeneity and Metastasis at Single-Cell Resolution.” Nature Cell Biology 20 (12): 1349–60. https://doi.org/10.1038/s41556-018-0236-7.

Leek, Jeffrey T., Robert B. Scharpf, Héctor Corrada Bravo, et al. 2010. “Tackling the Widespread and Critical Impact of Batch Effects in High-Throughput Data.” Nature Reviews Genetics 11 (10): 733–39. https://doi.org/10.1038/nrg2825.

Lever, Jake, Martin Krzywinski, and Naomi Altman. 2017. “Principal Component Analysis.” Nature Methods 14 (7): 641–42. https://doi.org/10.1038/nmeth.4346.

Li, Bo, and Colin N Dewey. 2011. “RSEM: Accurate Transcript Quantification from RNA-Seq Data with or Without a Reference Genome.” BMC Bioinformatics 12 (1). https://doi.org/10.1186/1471-2105-12-323.

Li, Heng. 2018. “Minimap2: Pairwise Alignment for Nucleotide Sequences.” Bioinformatics 34 (18): 3094–100. https://doi.org/10.1093/bioinformatics/bty191.

Liberzon, Arthur, Chet Birger, Helga Thorvaldsdóttir, Mahmoud Ghandi, Jill P. Mesirov, and Pablo Tamayo. 2015. “The Molecular Signatures Database Hallmark Gene Set Collection.” Cell Systems 1 (6): 417–25. https://doi.org/10.1016/j.cels.2015.12.004.

Liu, Yansheng, Andreas Beyer, and Ruedi Aebersold. 2016. “On the Dependency of Cellular Protein Levels on mRNA Abundance.” Cell 165 (3): 535–50. https://doi.org/10.1016/j.cell.2016.03.014.

Lockhart, David J., Helin Dong, Michael C. Byrne, et al. 1996. “Expression Monitoring by Hybridization to High-Density Oligonucleotide Arrays.” Nature Biotechnology 14 (13): 1675–80. https://doi.org/10.1038/nbt1296-1675.

Lopez, Romain, Jeffrey Regier, Michael B. Cole, Michael I. Jordan, and Nir Yosef. 2018. “Deep Generative Modeling for Single-Cell Transcriptomics.” Nature Methods 15 (12): 1053–58. https://doi.org/10.1038/s41592-018-0229-2.

Love, Michael I, Wolfgang Huber, and Simon Anders. 2014. “Moderated Estimation of Fold Change and Dispersion for RNA-Seq Data with DESeq2.” Genome Biology 15 (12). https://doi.org/10.1186/s13059-014-0550-8.

Lovén, Jakob, David A. Orlando, Alla A. Sigova, et al. 2012. “Revisiting Global Gene Expression Analysis.” Cell 151 (3): 476–82. https://doi.org/10.1016/j.cell.2012.10.012.

Luecken, Malte D., M. Büttner, K. Chaichoompu, et al. 2021. “Benchmarking Atlas-Level Data Integration in Single-Cell Genomics.” Nature Methods 19 (1): 41–50. https://doi.org/10.1038/s41592-021-01336-8.

Luecken, Malte D, and Fabian J Theis. 2019. “Current Best Practices in Single-Cell RNA-Seq Analysis: A Tutorial.” Molecular Systems Biology 15 (6). https://doi.org/10.15252/msb.20188746.

Lun, Aaron T. L., Karsten Bach, and John C. Marioni. 2016. “Pooling Across Cells to Normalize Single-Cell RNA Sequencing Data with Many Zero Counts.” Genome Biology 17 (1). https://doi.org/10.1186/s13059-016-0947-7.

Lund, Steven P., Dan Nettleton, Davis J. McCarthy, and Gordon K. Smyth. 2012. “Detecting Differential Expression in RNA-Sequence Data Using Quasi-Likelihood with Shrunken Dispersion Estimates.” Statistical Applications in Genetics and Molecular Biology 11 (5). https://doi.org/10.1515/1544-6115.1826.

Maaten, Laurens van der, and Geoffrey Hinton. 2008. “Visualizing Data Using t-SNE.” Journal of Machine Learning Research 9: 2579–605. https://www.jmlr.org/papers/v9/vandermaaten08a.html.

Macosko, Evan Z., Anindita Basu, Rahul Satija, et al. 2015. “Highly Parallel Genome-Wide Expression Profiling of Individual Cells Using Nanoliter Droplets.” Cell 161 (5): 1202–14. https://doi.org/10.1016/j.cell.2015.05.002.

Madeira, S. C., and A. L. Oliveira. 2004. “Biclustering Algorithms for Biological Data Analysis: A Survey.” IEEE/ACM Transactions on Computational Biology and Bioinformatics 1 (1): 24–45. https://doi.org/10.1109/tcbb.2004.2.

MAQC Consortium. 2006. “The MicroArray Quality Control (MAQC) Project Shows Inter- and Intraplatform Reproducibility of Gene Expression Measurements.” Nature Biotechnology 24 (9): 1151–61. https://doi.org/10.1038/nbt1239.

Marioni, John C., Christopher E. Mason, Shrikant M. Mane, Matthew Stephens, and Yoav Gilad. 2008. “RNA-Seq: An Assessment of Technical Reproducibility and Comparison with Gene Expression Arrays.” Genome Research 18 (9): 1509–17. https://doi.org/10.1101/gr.079558.108.

Marx, Vivien. 2021. “Method of the Year: Spatially Resolved Transcriptomics.” Nature Methods 18 (1): 9–14. https://doi.org/10.1038/s41592-020-01033-y.

McCarthy, Davis J., and Gordon K. Smyth. 2009. “Testing Significance Relative to a Fold-Change Threshold Is a TREAT.” Bioinformatics 25 (6): 765–71. https://doi.org/10.1093/bioinformatics/btp053.

McGranahan, Nicholas, and Charles Swanton. 2017. “Clonal Heterogeneity and Tumor Evolution: Past, Present, and the Future.” Cell 168 (4): 613–28. https://doi.org/10.1016/j.cell.2017.01.018.

McInnes, Leland, John Healy, and James Melville. 2018. UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction. https://arxiv.org/abs/1802.03426.

Minussi, Darlan C., Michael D. Nicholson, Hanghui Ye, et al. 2021. “Breast Tumours Maintain a Reservoir of Subclonal Diversity During Expansion.” Nature 592 (7853): 302–8. https://doi.org/10.1038/s41586-021-03357-x.

Monti, Stefano, Pablo Tamayo, Jill Mesirov, and Todd Golub. 2003. “Consensus Clustering: A Resampling-Based Method for Class Discovery and Visualization of Gene Expression Microarray Data.” Machine Learning 52 (1-2): 91–118. https://doi.org/10.1023/a:1023949509487.

Mortazavi, Ali, Brian A Williams, Kenneth McCue, Lorian Schaeffer, and Barbara Wold. 2008. “Mapping and Quantifying Mammalian Transcriptomes by RNA-Seq.” Nature Methods 5 (7): 621–28. https://doi.org/10.1038/nmeth.1226.

Nagalakshmi, Ugrappa, Zhong Wang, Karl Waern, et al. 2008. “The Transcriptional Landscape of the Yeast Genome Defined by RNA Sequencing.” Science 320 (5881): 1344–49. https://doi.org/10.1126/science.1158441.

Nam, Anna S., Ronan Chaligne, and Dan A. Landau. 2020. “Integrating Genetic and Non-Genetic Determinants of Cancer Evolution by Single-Cell Multi-Omics.” Nature Reviews Genetics 22 (1): 3–18. https://doi.org/10.1038/s41576-020-0265-5.

Navin, Nicholas, Jude Kendall, Jennifer Troge, et al. 2011. “Tumour Evolution Inferred by Single-Cell Sequencing.” Nature 472 (7341): 90–94. https://doi.org/10.1038/nature09807.

Newman, Aaron M, Chih Long Liu, Michael R Green, et al. 2015. “Robust Enumeration of Cell Subsets from Tissue Expression Profiles.” Nature Methods 12 (5): 453–57. https://doi.org/10.1038/nmeth.3337.

Nieuwenhuis, Sander, Birte U Forstmann, and Eric-Jan Wagenmakers. 2011. “Erroneous Analyses of Interactions in Neuroscience: A Problem of Significance.” Nature Neuroscience 14 (9): 1105–7. https://doi.org/10.1038/nn.2886.

Noble, William S. 2009. “How Does Multiple Testing Correction Work?” Nature Biotechnology 27 (12): 1135–37. https://doi.org/10.1038/nbt1209-1135.

Nowell, Peter C. 1976. “The Clonal Evolution of Tumor Cell Populations: Acquired Genetic Lability Permits Stepwise Selection of Variant Sublines and Underlies Tumor Progression.” Science 194 (4260): 23–28. https://doi.org/10.1126/science.959840.

Nuzzo, Regina. 2014. “Scientific Method: Statistical Errors.” Nature 506 (7487): 150–52. https://doi.org/10.1038/506150a.

Pardo-Palacios, Francisco J., Dingjie Wang, Fairlie Reese, et al. 2024. “Systematic Assessment of Long-Read RNA-Seq Methods for Transcript Identification and Quantification.” Nature Methods 21 (7): 1349–63. https://doi.org/10.1038/s41592-024-02298-3.

Patel, Anoop P., Itay Tirosh, John J. Trombetta, et al. 2014. “Single-Cell RNA-Seq Highlights Intratumoral Heterogeneity in Primary Glioblastoma.” Science 344 (6190): 1396–401. https://doi.org/10.1126/science.1254257.

Patro, Rob, Geet Duggal, Michael I Love, Rafael A Irizarry, and Carl Kingsford. 2017. “Salmon Provides Fast and Bias-Aware Quantification of Transcript Expression.” Nature Methods 14 (4): 417–19. https://doi.org/10.1038/nmeth.4197.

Perou, Charles M., Therese Sørlie, Michael B. Eisen, et al. 2000. “Molecular Portraits of Human Breast Tumours.” Nature 406 (6797): 747–52. https://doi.org/10.1038/35021093.

Pertea, Mihaela, Geo M Pertea, Corina M Antonescu, Tsung-Cheng Chang, Joshua T Mendell, and Steven L Salzberg. 2015. “StringTie Enables Improved Reconstruction of a Transcriptome from RNA-Seq Reads.” Nature Biotechnology 33 (3): 290–95. https://doi.org/10.1038/nbt.3122.

Picelli, Simone, Åsa K Björklund, Omid R Faridani, Sven Sagasser, Gösta Winberg, and Rickard Sandberg. 2013. “Smart-Seq2 for Sensitive Full-Length Transcriptome Profiling in Single Cells.” Nature Methods 10 (11): 1096–98. https://doi.org/10.1038/nmeth.2639.

Quackenbush, John. 2001. “Computational Analysis of Microarray Data.” Nature Reviews Genetics 2 (6): 418–27. https://doi.org/10.1038/35076576.

Regev, Aviv, Sarah A Teichmann, Eric S Lander, et al. 2017. “The Human Cell Atlas.” eLife 6. https://doi.org/10.7554/elife.27041.

Robinson, Mark D., Davis J. McCarthy, and Gordon K. Smyth. 2009. “edgeR : A Bioconductor Package for Differential Expression Analysis of Digital Gene Expression Data.” Bioinformatics 26 (1): 139–40. https://doi.org/10.1093/bioinformatics/btp616.

Robinson, Mark D, and Alicia Oshlack. 2010. “A Scaling Normalization Method for Differential Expression Analysis of RNA-Seq Data.” Genome Biology 11 (3). https://doi.org/10.1186/gb-2010-11-3-r25.

Rodriques, Samuel G., Robert R. Stickels, Aleksandrina Goeva, et al. 2019. “Slide-Seq: A Scalable Technology for Measuring Genome-Wide Expression at High Spatial Resolution.” Science 363 (6434): 1463–67. https://doi.org/10.1126/science.aaw1219.

Rousseeuw, Peter J. 1987. “Silhouettes: A Graphical Aid to the Interpretation and Validation of Cluster Analysis.” Journal of Computational and Applied Mathematics 20: 53–65. https://doi.org/10.1016/0377-0427(87)90125-7.

Schena, Mark, Dari Shalon, Ronald W. Davis, and Patrick O. Brown. 1995. “Quantitative Monitoring of Gene Expression Patterns with a Complementary DNA Microarray.” Science 270 (5235): 467–70. https://doi.org/10.1126/science.270.5235.467.

Schurch, Nicholas J., Pietá Schofield, Marek Gierliński, et al. 2016. “How Many Biological Replicates Are Needed in an RNA-Seq Experiment and Which Differential Expression Tool Should You Use?” RNA 22 (6): 839–51. https://doi.org/10.1261/rna.053959.115.

Schwanhäusser, Björn, Dorothea Busse, Na Li, et al. 2011. “Global Quantification of Mammalian Gene Expression Control.” Nature 473 (7347): 337–42. https://doi.org/10.1038/nature10098.

Schwartz, Russell, and Alejandro A. Schäffer. 2017. “The Evolution of Tumour Phylogenetics: Principles and Practice.” Nature Reviews Genetics 18 (4): 213–29. https://doi.org/10.1038/nrg.2016.170.

Simon, R., M. D. Radmacher, K. Dobbin, and L. M. McShane. 2003. “Pitfalls in the Use of DNA Microarray Data for Diagnostic and Prognostic Classification.” JNCI Journal of the National Cancer Institute 95 (1): 14–18. https://doi.org/10.1093/jnci/95.1.14.

Smyth, Gordon K. 2004. “Linear Models and Empirical Bayes Methods for Assessing Differential Expression in Microarray Experiments.” Statistical Applications in Genetics and Molecular Biology 3 (1): 1–25. https://doi.org/10.2202/1544-6115.1027.

Soneson, Charlotte, Michael I. Love, and Mark D. Robinson. 2016. “Differential Analyses for RNA-Seq: Transcript-Level Estimates Improve Gene-Level Inferences.” F1000Research 4: 1521. https://doi.org/10.12688/f1000research.7563.2.

Squair, Jordan W., Matthieu Gautier, Claudia Kathe, et al. 2021. “Confronting False Discoveries in Single-Cell Differential Expression.” Nature Communications 12 (1). https://doi.org/10.1038/s41467-021-25960-2.

Ståhl, Patrik L., Fredrik Salmén, Sanja Vickovic, et al. 2016. “Visualization and Analysis of Gene Expression in Tissue Sections by Spatial Transcriptomics.” Science 353 (6294): 78–82. https://doi.org/10.1126/science.aaf2403.

Stark, Rory, Marta Grzelak, and James Hadfield. 2019. “RNA Sequencing: The Teenage Years.” Nature Reviews Genetics 20 (11): 631–56. https://doi.org/10.1038/s41576-019-0150-2.

Stoeckius, Marlon, Christoph Hafemeister, William Stephenson, et al. 2017. “Simultaneous Epitope and Transcriptome Measurement in Single Cells.” Nature Methods 14 (9): 865–68. https://doi.org/10.1038/nmeth.4380.

Storey, John D., and Robert Tibshirani. 2003. “Statistical Significance for Genomewide Studies.” Proceedings of the National Academy of Sciences 100 (16): 9440–45. https://doi.org/10.1073/pnas.1530509100.

Su, Andrew I., Tim Wiltshire, Serge Batalov, et al. 2004. “A Gene Atlas of the Mouse and Human Protein-Encoding Transcriptomes.” Proceedings of the National Academy of Sciences 101 (16): 6062–67. https://doi.org/10.1073/pnas.0400782101.

Subramanian, Aravind, Pablo Tamayo, Vamsi K. Mootha, et al. 2005. “Gene Set Enrichment Analysis: A Knowledge-Based Approach for Interpreting Genome-Wide Expression Profiles.” Proceedings of the National Academy of Sciences 102 (43): 15545–50. https://doi.org/10.1073/pnas.0506580102.

Svensson, Valentine. 2020. “Droplet scRNA-Seq Is Not Zero-Inflated.” Nature Biotechnology 38 (2): 147–50. https://doi.org/10.1038/s41587-019-0379-5.

Svensson, Valentine, Roser Vento-Tormo, and Sarah A Teichmann. 2018. “Exponential Scaling of Single-Cell RNA-Seq in the Past Decade.” Nature Protocols 13 (4): 599–604. https://doi.org/10.1038/nprot.2017.149.

Tang, Fuchou, Catalin Barbacioru, Yangzhou Wang, et al. 2009. “mRNA-Seq Whole-Transcriptome Analysis of a Single Cell.” Nature Methods 6 (5): 377–82. https://doi.org/10.1038/nmeth.1315.

The GTEx Consortium. 2020. “The GTEx Consortium Atlas of Genetic Regulatory Effects Across Human Tissues.” Science 369 (6509): 1318–30. https://doi.org/10.1126/science.aaz1776.

The Tabula Muris Consortium. 2018. “Single-Cell Transcriptomics of 20 Mouse Organs Creates a Tabula Muris.” Nature 562 (7727): 367–72. https://doi.org/10.1038/s41586-018-0590-4.

The Tabula Sapiens Consortium. 2022. “The Tabula Sapiens: A Multiple-Organ, Single-Cell Transcriptomic Atlas of Humans.” Science 376 (6594). https://doi.org/10.1126/science.abl4896.

Theodoris, Christina V., Ling Xiao, Anant Chopra, et al. 2023. “Transfer Learning Enables Predictions in Network Biology.” Nature 618 (7965): 616–24. https://doi.org/10.1038/s41586-023-06139-9.

Tibshirani, Robert, Trevor Hastie, Balasubramanian Narasimhan, and Gilbert Chu. 2002. “Diagnosis of Multiple Cancer Types by Shrunken Centroids of Gene Expression.” Proceedings of the National Academy of Sciences 99 (10): 6567–72. https://doi.org/10.1073/pnas.082099299.

Tibshirani, Robert, Guenther Walther, and Trevor Hastie. 2001. “Estimating the Number of Clusters in a Data Set via the Gap Statistic.” Journal of the Royal Statistical Society Series B: Statistical Methodology 63 (2): 411–23. https://doi.org/10.1111/1467-9868.00293.

Tirosh, Itay, Benjamin Izar, Sanjay M. Prakadan, et al. 2016. “Dissecting the Multicellular Ecosystem of Metastatic Melanoma by Single-Cell RNA-Seq.” Science 352 (6282): 189–96. https://doi.org/10.1126/science.aad0501.

Traag, V. A., L. Waltman, and N. J. van Eck. 2019. “From Louvain to Leiden: Guaranteeing Well-Connected Communities.” Scientific Reports 9 (1). https://doi.org/10.1038/s41598-019-41695-z.

Trapnell, Cole, Davide Cacchiarelli, Jonna Grimsby, et al. 2014. “The Dynamics and Regulators of Cell Fate Decisions Are Revealed by Pseudotemporal Ordering of Single Cells.” Nature Biotechnology 32 (4): 381–86. https://doi.org/10.1038/nbt.2859.

Trapnell, Cole, Lior Pachter, and Steven L. Salzberg. 2009. “TopHat: Discovering Splice Junctions with RNA-Seq.” Bioinformatics 25 (9): 1105–11. https://doi.org/10.1093/bioinformatics/btp120.

Trapnell, Cole, Brian A Williams, Geo Pertea, et al. 2010. “Transcript Assembly and Quantification by RNA-Seq Reveals Unannotated Transcripts and Isoform Switching During Cell Differentiation.” Nature Biotechnology 28 (5): 511–15. https://doi.org/10.1038/nbt.1621.

Uhlén, Mathias, Linn Fagerberg, Björn M. Hallström, et al. 2015. “Tissue-Based Map of the Human Proteome.” Science 347 (6220). https://doi.org/10.1126/science.1260419.

Veer, Laura J. van ’t, Hongyue Dai, Marc J. van de Vijver, et al. 2002. “Gene Expression Profiling Predicts Clinical Outcome of Breast Cancer.” Nature 415 (6871): 530–36. https://doi.org/10.1038/415530a.

Velculescu, Victor E., Lin Zhang, Bert Vogelstein, and Kenneth W. Kinzler. 1995. “Serial Analysis of Gene Expression.” Science 270 (5235): 484–87. https://doi.org/10.1126/science.270.5235.484.

Venet, David, Jacques E. Dumont, and Vincent Detours. 2011. “Most Random Gene Expression Signatures Are Significantly Associated with Breast Cancer Outcome.” PLoS Computational Biology 7 (10): e1002240. https://doi.org/10.1371/journal.pcbi.1002240.

Wagner, Günter P., Koryu Kin, and Vincent J. Lynch. 2012. “Measurement of mRNA Abundance Using RNA-Seq Data: RPKM Measure Is Inconsistent Among Samples.” Theory in Biosciences 131 (4): 281–85. https://doi.org/10.1007/s12064-012-0162-3.

Wang, Zhong, Mark Gerstein, and Michael Snyder. 2009. “RNA-Seq: A Revolutionary Tool for Transcriptomics.” Nature Reviews Genetics 10 (1): 57–63. https://doi.org/10.1038/nrg2484.

Wasserstein, Ronald L., and Nicole A. Lazar. 2016. “The ASA Statement on p -Values: Context, Process, and Purpose.” The American Statistician 70 (2): 129–33. https://doi.org/10.1080/00031305.2016.1154108.

Wolf, F. Alexander, Philipp Angerer, and Fabian J. Theis. 2018. “SCANPY: Large-Scale Single-Cell Gene Expression Data Analysis.” Genome Biology 19 (1). https://doi.org/10.1186/s13059-017-1382-0.

You, Yue, Yuting Fu, Lanxiang Li, et al. 2024. “Systematic Comparison of Sequencing-Based Spatial Transcriptomic Methods.” Nature Methods 21 (9): 1743–54. https://doi.org/10.1038/s41592-024-02325-3.

Zeng, Hongkui. 2022. “What Is a Cell Type and How to Define It?” Cell 185 (15): 2739–55. https://doi.org/10.1016/j.cell.2022.06.031.

Zheng, Grace X. Y., Jessica M. Terry, Phillip Belgrader, et al. 2017. “Massively Parallel Digital Transcriptional Profiling of Single Cells.” Nature Communications 8 (1). https://doi.org/10.1038/ncomms14049.

Zimmerman, Kip D., Mark A. Espeland, and Carl D. Langefeld. 2021. “A Practical Solution to Pseudoreplication Bias in Single-Cell Studies.” Nature Communications 12 (1). https://doi.org/10.1038/s41467-021-21038-1.