DiPhyxDiPhyx
Workspace
The shape of the problem

Twenty thousand tests,and six samples

Genomic data is wide and shallow. Almost every statistical method in the field exists to survive that shape.

Many features, few samples

Tens of thousands of genes against a handful of replicates. Classical statistics assumes the opposite.

Counts, not measurements

Discrete, non-negative, and with variance that grows with the mean. Normal-theory tests do not fit.

Structure you did not intend

Batch, sex, ancestry, run date. All of it correlates with something, and none of it is what you asked about.

Multiple testing

At p < 0.05,a thousand genes lie

Run twenty thousand tests on pure noise and a thousand come back significant. The first thing to correct for, and the easiest to botch.

BonferroniDivide the threshold by the number of tests. Correct, extremely conservative, and rarely right for a screen.
Benjamini-HochbergControls the false discovery rate — the fraction of your hits that are wrong. The standard for genomics.
Filter before you correctDropping genes with almost no counts raises power, because you are no longer paying for tests that could not resolve.

Bonferroni controls the chance of any false positive at all, which is the right goal for a handful of confirmatory tests and far too strict for a screen. Benjamini-Hochberg controls the proportion of your hits that are false, which is what you actually want when the output is a candidate list. Report the adjusted value, and say which method produced it.

Counts are not normal

Why RNA-seqneeds its own tests

Read counts are discrete, non-negative, and more variable than a Poisson allows. A t-test on them is the wrong tool.

Negative binomial

Poisson plus an over-dispersion term, which is what real biological replicates actually look like.

Shared dispersion

Estimates are borrowed across genes, which is the trick that makes three replicates usable.

Give it raw counts

DESeq2 and edgeR normalise internally. Feeding them TPM or FPKM invalidates the model.

DESeq2 and edgeR both model counts as negative binomial and both borrow information across genes to estimate dispersion, which is what makes three replicates workable at all. Give them raw counts. They normalise internally, and handing them TPM or FPKM breaks the assumption the whole model rests on.

Before any of that

Batch effectsoutrank your biology

The largest source of variation in a genomics dataset is frequently the day it was processed, and no test will tell you that.

Plot PCAFirstbefore any hypothesis test
By batchIf it clustersthat is your dominant signal
In the designFixconfounded batches cannot be modelled out

Plot the first two principal components before testing anything. If samples cluster by processing date or sequencing run rather than by condition, that is your dominant signal. It can be modelled — include batch as a covariate — but only if the design allows it. A batch perfectly confounded with your condition cannot be rescued by statistics, only by a better experiment.

Effect size, not just pWith enough samples everything is significant. Report the fold change beside the adjusted p-value.
Decide the test firstChoosing the method after seeing the result is how a screen produces a finding that never replicates.
Enrichment is not evidenceA pathway term coming up significant is a hypothesis. Treat it as the start of the next experiment.
On a machine

RStudio Server isalready built

Nothing here needs compiling. Pick a machine, and these two commands put it on there.

hub://rstudio

$ dxflow workflow create --identity rstudio hub://rstudio

$ dxflow workflow start rstudio

Start a machineWhat RStudio Server needs
Read it, then run itPrepaid and by the hour, on a machine that is yours about a minute after you ask.
Start a machineAll notes