Reading a Differential Expression Table: Why the Missing Rows Carry as Much Information as the Ones You Got

Reading a Differential Expression Table: Why the Missing Rows Carry as Much Information as the Ones You Got
Zetobit Reading the Report
NA NA NA NA READING THE REPORT Reading a Differential Expression Table Why the missing rows carry as much information as the ones you got Kanna Nandakumar, PhD Zetobit
Reading the Report

Reading a Differential Expression Table: Why the Missing Rows Carry as Much Information as the Ones You Got

The rows you received are a ranked list of evidence, not a ranked list of biology — and the rows that are absent are absent for four different reasons the file does not distinguish.

A spreadsheet arrives with a few hundred rows in it. Each row is a gene, and beside it are two numbers: how much it moved, and how sure the analysis is that it moved. The file looks like a finished answer — the list of things that changed. Most of the conversations that follow treat it as one.

It is something narrower. The file ranks the genes for which the evidence of a difference was strongest, among the genes that were tested, under the comparison that was specified. Each of those three qualifiers does real work, and none of them appears in the columns you were sent.

Nothing that follows implies the analysis was done badly. Every behavior described here is documented, deliberate, and usually the statistically correct choice. The difficulty is only that the choices are invisible in the delivered file — and the delivered file is the thing that circulates, gets pasted into slides, and sets the direction of the next six months of work.

The ranking ranks confidence

Sort a results table by adjusted p-value, as nearly everyone does, and you have ordered the genes by how well the difference was measured. Not by how large it is, and not by how much it matters. The American Statistical Association's statement on p-values makes the point in a single line: statistical significance does not measure the size of an effect or the importance of a result.1

In count data, confidence tracks abundance. A gene sitting at 50,000 counts can produce a persuasive p-value from a 15% shift, because a 15% shift is precisely measurable at that depth. A transcription factor sitting at 30 counts can genuinely double and still land in the middle of the table, because at 30 counts a doubling is barely separable from sampling noise. The top of a significance-ranked table is therefore enriched for the abundant and the well-behaved — which in most tissues means metabolic and structural genes rather than the regulators that were probably the point of the experiment.

This is not a technicality. It changes which lists replicate. In the FDA-led MicroArray Quality Control work, gene lists built by ranking on fold change with a non-stringent p-value cutoff were more reproducible across laboratories and platforms than lists ranked on the test statistic.2 Ranking by significance is ranking partly by measurement quality, and measurement quality is a property of the site as well as the sample.

A blank is not a result

Now the harder half of the file: the genes that are not in it. Absence has at least four distinct causes, and the delivered table almost never tells them apart.

  1. The gene was never in the annotation used for quantification, so it was never counted at all.
  2. The gene was in the count matrix but had essentially no reads in any sample. No test was possible.
  3. The gene was tested and the evidence did not clear the threshold.
  4. The gene was tested, produced a raw p-value, and was then removed from the multiple-testing correction before an adjusted value was ever calculated.

The fourth surprises people. DESeq2's own documentation states it plainly: a row filtered out for having a low mean normalized count has only its adjusted p-value set to NA, while a row containing an extreme count outlier has both the p-value and the adjusted value set to NA.3 Three different upstream events all print as an empty cell — and if whoever prepared the file dropped the NA rows before sending it, they print as nothing at all.

RESULTS TABLE — EXCERPT AS DELIVERED Illustrative. Gene identifiers replaced. gene baseMean log2FC lfcSE pvalue padj GENE A 48210.4 0.42 0.09 1.1e-06 8.4e-04 Top of the list Ranked first because it is measured well, not because it moved much: 34%. GENE B 1205.7 2.87 0.61 3.0e-05 9.7e-03 GENE C 318.2 -1.94 0.88 2.4e-02 1.9e-01 Tested, not significant The only absence that means what people assume it means. GENE D 6.1 3.10 2.44 8.1e-02 NA Filtered before correction Low mean count: a raw p-value exists, an adjusted one never did. GENE E 0.0 NA NA NA NA Never testable — no reads at all. GENE F 940.3 -0.77 0.20 NA NA Set aside as an outlier One extreme sample. Well expressed, plausible, no result. ... WHAT THE FILE IS A SUBSET OF 21,806 rows in the count matrix 13,241 rows that received a padj 214 rows you were sent Both narrowings happened before you opened the file.
An excerpt of a differential expression results table, with the four ways a gene can fail to appear marked in the rows where they occur. The strip beneath shows the two narrowings that precede delivery: the statistical filter, and the export.

The set of genes tested was chosen using the data

The filter in the fourth case is worth understanding, because it behaves in a way that catches people out. The expression threshold is not fixed in advance. It is selected to maximize the number of genes that pass at a stated significance level — a level whose default value in DESeq2 is 0.10, not the 0.05 most reports quote.3

The boundary of the table is a function of the whole dataset. Change the samples in the comparison and the boundary moves, without any individual gene's numbers changing much.

That has a practical consequence you have probably already experienced. A gene appears with an adjusted p-value in one analysis and comes back blank in the next, after the analyst dropped two samples, added a covariate, or split the cohort. Nothing went wrong and nothing about that gene changed. The population of genes eligible for correction changed, and it took the gene with it. That disappearance is worth recognizing on sight, because it is otherwise very easy to interpret as biology.

What the design could see at all

Before any of this, the table is bounded by replication. The largest experiment of its kind — 48 biological replicates per condition in yeast — allows a direct measurement of what smaller designs miss. Analyses using three replicates recovered only 20–40% of the genes identified with the full set of 42 clean replicates. For genes changing by more than fourfold, three replicates recovered over 85%. Recovering more than 85% of differentially expressed genes regardless of fold change required more than 20 replicates.4

Use that as a translation table for your own file. A result set from a three-versus-three design is, approximately, the large-effect subset of the real answer. It is not wrong; it is partial in a specific and predictable direction, and what it omits is disproportionately the modest, regulated change that most biology actually consists of. The corollary matters most in the sentence people say without thinking: we saw no effect on pathway X. At three replicates per group, that sentence is largely a statement about the number three. The Nature comment signed by more than 800 researchers in 2019 made the general form of the argument — that reading non-significance as evidence of no difference routinely dismisses effects that may be crucial.5

The fold change has a version too

One last column. The log fold change you were given may or may not be the raw ratio between group means. Genes with low counts produce wild ratio estimates — a gene going from 1 count to 8 is an eightfold change and also nothing at all — so the standard practice is to shrink those estimates toward zero using information shared across genes, an approach built into DESeq2 since 2014.6 Both the raw and the shrunken estimate are legitimate outputs of the same run on the same model, and for a low-expressed gene they can differ severalfold. If you are applying a fold-change cutoff, comparing against a qPCR result, or putting a number into a figure legend, which one you are holding is not a detail.

What survives all of this is still useful, and it is worth being precise about what it is. A differential expression table is a well-constructed device for ordering hypotheses: it tells you where, in a space of twenty thousand genes, the evidence for a difference was concentrated under one specific comparison. That is a genuinely valuable thing to be handed. It is simply not a census of what changed, and the three or four decisions separating those two readings all happened before the file reached you.

Three questions to ask

  1. Can I have the unfiltered table — every gene tested, with raw p-values and blanks preserved? The file you were given has had its boundary drawn twice: once by the statistical filter, once by whoever exported the significant rows. Only the full table tells you which absences mean what.
  2. What exactly was compared, and what else was in the model? The direction of every fold change depends on which group was set as the reference. A model that adjusts for batch, sex, or donor answers a different question than one that does not. Neither fact is visible in the columns.
  3. How many independent biological replicates, and is the fold-change column shrunken? The first tells you the size of effect the experiment was able to detect. The second tells you whether the magnitudes can be compared with anything else you are holding.

References

  1. Wasserstein RL, Lazar NA. The ASA Statement on p-Values: Context, Process, and Purpose. The American Statistician. 2016;70(2):129–133. tandfonline.com
  2. Guo L, Lobenhofer EK, Wang C, et al. Rat toxicogenomic study reveals analytical consistency across microarray platforms. Nature Biotechnology. 2006;24:1162–1169. nature.com
  3. Love MI, Anders S, Huber W. Analyzing RNA-seq data with DESeq2 — package vignette and results() documentation. Bioconductor. bioconductor.org
  4. Schurch NJ, Schofield P, Gierliński M, et al. How many biological replicates are needed in an RNA-seq experiment and which differential expression tool should you use? RNA. 2016;22(6):839–851. rnajournal.cshlp.org
  5. Amrhein V, Greenland S, McShane B. Scientists rise up against statistical significance. Nature. 2019;567:305–307. nature.com
  6. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology. 2014;15:550. genomebiology.biomedcentral.com
Reading the Report is a Zetobit series for people who receive analytical output rather than produce it. Each piece takes one document as it arrives and asks what it establishes, what it leaves open, and what to ask next. Figures are illustrative and do not reproduce any client result.
Previous
Previous

“We Ran the Standard Pipeline”: Four Questions That Make That Sentence Mean Something

Next
Next

Reading a Pathway Enrichment Figure: Why Twenty Bars Are Not Twenty Findings