Summary Statistics for Genomic Data

A genomic summary statistic is any descriptive quantity — mean, median, variance, standard deviation, coefficient of variation — computed per gene (or per genomic feature) across a set of samples, rather than computed once across an entire dataset. In a typical RNA-seq experiment, the raw data table has one row per gene and one column per sample, with entries representing the number of sequencing reads mapped to that gene (a read count), or a normalized version of that count such as TPM (transcripts per million) or CPM (counts per million), which correct for differences in sequencing depth and gene length between samples. Because there are thousands of genes, genomic analysis means computing thousands of parallel univariate summaries, not one. The statistical difficulty is that raw counts are almost never symmetric. A gene expressed at low, noisy levels might read 0, 1, 3, 0, 2 across five replicate samples; a highly expressed housekeeping gene might read 8,500, 9,200, 7,800, 11,000, 9,600. Both distributions are right-skewed to different degrees, and their variances scale with their means — a property called mean-variance dependence, formalized in RNA-seq statistics by the observation that count variance is approximately proportional to the mean plus a term proportional to the mean squared (the negative binomial variance structure used by tools such as edgeR and DESeq2). This is the opposite of the constant-variance assumption behind classical statistics such as the ordinary $t$-test, and it means a raw arithmetic mean and standard deviation, computed without correction, systematically misrepresent low-count genes as "more variable" purely because they are lowly expressed, not because they are biologically unstable. The standard remedy is the logarithmic transformation: replacing each expression value $x$ with $\log_2(x + c)$, where $c$ is a small constant (often 1) added to avoid taking the log of zero. On the log scale, a two-fold increase in expression — from 100 to 200, or from 10,000 to 20,000 — becomes the same additive step (an increase of exactly 1, since $\log_2 2 = 1$), regardless of the gene's baseline expression level. This is the key property that justifies the transform: biological regulation acts multiplicatively (a transcription factor doubling a promoter's activity doubles expression whether the gene starts at 10 copies or 10,000), and the log transform converts this multiplicative structure into an additive one that the mean, variance, and coefficient of variation were built to describe.

Step 1 — Raw counts and normalization. Sequencing produces integer read counts per gene per sample. These counts depend not only on true expression but on total sequencing depth (a sample sequenced twice as deeply gets roughly twice the counts for every gene) and on transcript length (longer transcripts capture more reads). Before any descriptive statistic is meaningful, counts must be normalized — commonly to CPM (dividing by total library size and multiplying by $10^6$) or TPM (which also corrects for gene length). Skipping this step confounds technical sequencing depth with biological expression differences, the single most common error in early genomic pipelines.

Step 2 — Per-gene distributions across samples. For each gene $g$, collect its normalized expression value across all $n$ samples: $x_{g1}, x_{g2}, \dots, x_{gn}$. This is the raw material for a per-gene mean $\bar{x}_g$, variance $s_g^2$, and coefficient of variation $CV_g = s_g / \bar{x}_g$. Because a typical dataset has 20,000+ genes, this produces 20,000+ parallel distributions rather than one.

Step 3 — Log transformation. Apply $y_{gi} = \log_2(x_{gi} + 1)$ to every value. The pseudocount of 1 is a pragmatic compromise: it makes $\log_2(0+1) = 0$ so zero counts do not produce undefined values, but it also compresses the extreme low end nonlinearly, a distortion analysts must keep in mind for very lowly expressed genes.

Step 4 — Summary statistics on the log scale. Compute $\bar{y}_g$ and $s_{y,g}$ on the transformed data. Because addition on the log scale corresponds to multiplication on the raw scale, differences in $\bar{y}_g$ between two conditions are directly interpretable as log fold-changes:

A $\log_2 FC$ of $+1$ means a doubling; $-1$ means a halving; $+3.32$ means a ten-fold increase ($\log_2 10 \approx 3.32$).

Step 5 — The coefficient of variation as a filtering tool. Before testing thousands of genes for differential expression, analysts often rank genes by $CV_g$ (on the log scale, or equivalently by variance since means are more comparable after log transform) to identify genes with real, biologically meaningful variability across samples versus genes that are flat and uninformative — a step central to selecting "highly variable genes" in single-cell RNA-seq analysis, where thousands of genes are essentially noise around zero. A gene with $CV_g$ near the noise floor is unlikely to be a useful marker regardless of its mean.

Step 6 — Outlier detection at scale. With thousands of parallel distributions, a small fraction of genes will show extreme values by pure sampling variation (multiple-comparison logic exactly analogous to running thousands of simultaneous hypothesis tests). Genomic pipelines flag both sample-level outliers (an entire sample deviates from all others, often signaling technical failure — degraded RNA, contamination, sample swap) and gene-level outliers (a single measurement within a gene's distribution is implausibly extreme, often a PCR duplication artifact or a sequencing alignment error), and treat them differently: sample-level outliers are usually excluded entirely; gene-level single-point outliers are handled with robust statistics (median, trimmed mean) or explicit outlier-aware count models.

The case for logging expression data, and for treating gene-by-gene descriptive statistics carefully, rests on both a theoretical result and repeated empirical demonstrations from the early microarray era through modern RNA-seq.

The theoretical foundation: variance-stabilizing transformations. The statistical logic was formalized well before genomics existed — the delta method shows that if a variable's variance is proportional to the square of its mean (as is approximately true for expression counts, a consequence of multiplicative noise sources like PCR amplification efficiency and cell-to-cell variability), then the log transform makes the variance approximately constant across the range of the mean. This is not a genomics-specific trick; it is the same reasoning behind logging income data in economics or logging population sizes in ecology, wherever a quantity grows multiplicatively.

Microarray-era validation (Affymetrix and cDNA array studies, late 1990s–2000s). Early microarray analyses of thousands of genes simultaneously showed directly what happens without a log transform: housekeeping and highly expressed genes dominated variance rankings purely due to scale, while biologically dynamic but modestly expressed genes were statistically invisible. Studies comparing raw versus log-transformed intensities (for example, in analyses underlying the MIAME standards and early limma/Bioconductor methodology developed by Gordon Smyth and colleagues in the early 2000s) demonstrated that the log-scale residual variance was far more homogeneous across the expression range, and that statistical tests performed on log ratios gave far fewer spurious "differentially expressed" calls driven purely by scale-dependent noise.

RNA-seq count modeling (Anders & Huber, 2010; Love, Huber & Anders, 2014). The development of DESeq and later DESeq2 directly measured the mean–variance relationship across thousands of real RNA-seq genes from replicate samples and confirmed the negative binomial pattern: variance rises with the mean faster than a Poisson (purely counting) process would predict, due to extra biological variability between replicates. This empirical mean-variance curve is precisely why these tools use a variance-stabilizing transformation (closely related to $\log_2$) before clustering samples or computing per-gene summary statistics such as means used in quality-control plots, and precisely why raw-count PCA plots look qualitatively different — and less biologically interpretable — from log-transformed PCA plots on the same data.

The MAQC (MicroArray Quality Control) Consortium, 2006. This multi-lab study, published in Nature Biotechnology, compared expression measurements for the same RNA samples across multiple microarray platforms and RT-PCR benchmarks. A central finding was that agreement between platforms, and between technical replicates, was dramatically higher when comparisons were made on log-ratio (fold-change) scales rather than raw intensity scales — direct empirical evidence that the log transform captures the reproducible biological signal while raw-scale summary statistics amplify platform-specific technical noise. The MAQC results became a founding argument for standardizing log fold-change as the default unit of differential expression reporting across the field.

What these lines of evidence jointly established: (1) expression variance genuinely scales with the mean, confirmed empirically across microarray and RNA-seq platforms; (2) log transformation measurably stabilizes this variance and improves reproducibility across labs and platforms; (3) descriptive statistics computed on raw counts systematically misrepresent both variability and outlier status for lowly expressed genes. None of this is an assumption imported wholesale from other fields — it is a repeatedly re-verified empirical property of gene expression data itself.

Example 1 — Comparing two genes' variability using the CV, on the correct scale. Suppose two genes are measured across five tumor samples (normalized counts, CPM):

Raw means: $\bar{x}_A = 26$, $\bar{x}_B = 1100$. Raw standard deviations: $s_A \approx 15.7$, $s_B \approx 660$. Raw CVs:

The two genes have identical raw coefficients of variation — a coincidence chosen here to make a point: their raw scales differ 40-fold, but their relative variability, as measured by CV, is the same, exactly the situation where the CV is doing its job correctly by making magnitude-independent variability comparable. Now log-transform ($y=\log_2(x+1)$):

Log-scale means: $\bar{y}_A \approx 4.54$, $\bar{y}_B \approx 9.91$. Log-scale SDs: $s_{yA} \approx 0.82$, $s_{yB} \approx 0.75$. The log-scale spreads (0.82 vs 0.75) are close but not identical to the raw CVs (0.60 and 0.60) — this is expected, because for CV $< \sim$0.3–0.5 the two measures are approximately equal, but they diverge as the CV grows, since $s_y \approx CV$ only in the small-CV limit (a first-order Taylor approximation). Here both genes show fold-change-type spread of roughly $2^{0.8} \approx 1.7$-fold around their means — a genuinely comparable, scale-free statement that the raw means and SDs alone could not make transparent without first computing the CV.

Example 2 — A single outlier sample distorting a fold-change calculation. A researcher compares expression of a stress-response gene between four control and four treated samples (normalized counts): Control: 102, 98, 110, 95 Treated: 205, 215, 198, 9800 Naively computed raw means:

This looks like a dramatic, headline-worthy induction. But inspecting the treated group (a step this article insists must accompany any summary statistic, following the logic of Displaying Data Honestly) shows three values tightly clustered near 200 and one wildly discrepant value of 9800 — almost certainly a technical artifact (e.g., a PCR jackpot event, a contaminating sample, or an alignment error placing reads from a different gene into this one). Recomputing with the median, or with the outlier flagged and excluded per the diagnostic principles in Outliers in Biological Data (e.g., studentized residual far beyond expected range, or a value more than 5 median absolute deviations from the other three):

The true effect is a modest, biologically plausible two-fold induction, not a 25-fold outlier-driven artifact. On the log scale the distortion is naturally suppressed even before exclusion: $\log_2(9800+1) \approx 13.3$ versus $\log_2(205) \approx 7.68$, so the log-scale mean of the treated group is pulled toward the outlier far less severely than the raw mean is, though careful practice still flags and investigates the point rather than relying on the transform alone to fix it.

The modern practice of log-transforming and rigorously summarizing genomic data does not trace to a single discoverer in the way Mendelian ratios trace to Mendel; it emerged from statisticians and computational biologists who adapted classical descriptive statistics to the specific pathologies of high-throughput count data. The figure most closely associated with formalizing variance-stabilization and moderated statistics for genomic data is Gordon K. Smyth, whose limma framework (Bioconductor, first released 2004, with the foundational empirical Bayes moderated statistics paper published in 2004–2005) built the standard statistical machinery — including log-ratio summaries and per-gene variance moderation — used across essentially all microarray and much RNA-seq differential expression analysis since. Smyth, a statistician at the Walter and Eliza Hall Institute of Medical Research in Melbourne, Australia, worked at the interface of classical linear-model theory and the specific small-sample, many-variable structure of genomic experiments — often only 3–6 replicate arrays per condition, but tens of thousands of genes. His key methodological insight, building on empirical Bayes shrinkage ideas, was that per-gene variance estimates from so few replicates are themselves unreliable, and can be improved by "borrowing strength" across all genes simultaneously — a genomics-scale extension of the same logic that motivates using pooled, robust, and shrinkage-based summary statistics rather than naive per-gene numbers taken at face value.

Cancer genomics and precision oncology. Tumor expression profiling (e.g., in breast cancer subtyping via the PAM50 classifier) depends on correctly summarized, log-scale expression signatures to assign patients to molecular subtypes that guide chemotherapy versus targeted therapy decisions; a raw-scale summary would misclassify tumors due to scale artifacts unrelated to biology.

Single-cell RNA-seq and cell-type discovery. Modern single-cell atlases (e.g., the Human Cell Atlas) rely on per-gene coefficient of variation and log-normalized variance to select "highly variable genes" before clustering cells into types — a direct, large-scale descendant of the CV-based filtering described here, now applied across tens of thousands of individual cells rather than bulk tissue samples.

Diagnostic gene expression panels. Clinical assays such as Oncotype DX (breast cancer recurrence risk) and similar multi-gene panels report results as normalized, log-scale expression scores precisely because raw fluorescence or count units are not portable or reproducible across laboratories and instruments — the same reproducibility argument the MAQC study established at genome scale.

Agricultural and crop genomics. Expression profiling of stress-response genes in crop breeding programs (drought or pathogen resistance in maize, rice, wheat) uses the same log-fold-change framework to compare resistant versus susceptible cultivars, guiding marker-assisted selection.

Conservation genomics. Transcriptomic monitoring of wildlife populations (e.g., coral bleaching stress-response studies, salmon thermal-tolerance studies) applies identical per-gene summary statistics to identify biomarkers of environmental stress, where sample sizes are often even smaller than in clinical studies, making robust, correctly-scaled summary statistics especially critical.

Public health and infectious disease. Host transcriptomic signatures distinguishing bacterial from viral infection (used in some point-of-care diagnostic development) are built from log-scale differential expression summaries across curated gene panels, directly applying the descriptive statistics machinery described above to a clinical decision problem.

Ask the tutor: How would you decide whether to filter a low-count gene out of an analysis versus keep it and just interpret its statistics cautiously?