Checking q values against a known truth

Why we need a ground truth¶

In the previous chapter we estimated false discovery rates. The B6 against B6 comparison let us check that the procedure reports nothing when every gene follows the null hypothesis, but such a null experiment only checks the null side: it contains no true findings, and in the comparison of the two strains we do not know which genes are truly differentially expressed. A procedure that reports $q=0.01$ makes a promise about the long run frequency of errors among its findings, and that promise should be testable.

A spike-in experiment, which contains known true positives as well as true nulls, lets us test it. One prepares a complex but constant biological background, and adds a small number of foreign proteins at concentrations chosen, and therefore known, by the experimenter. Only the spiked-in proteins should differ between the samples. We can then run a complete differential abundance analysis and afterwards count how many of our findings were actually changed. The fraction of errors among the findings is the false discovery proportion (FDP),

which is what $\hat{\rm FDR}(t)$, and hence $\hat{q}$, is trying to predict. The FDR is an expectation over hypothetical repetitions of the experiment; the FDP is a number we can read off a single data set, if we know the truth.

The iPRG 2015 spike-in study¶

We will use the benchmark set of the ABRF Proteome Informatics Research Group, described by Choi et al., J. Proteome Res. 2017, and deposited at MassIVE under accession MSV000079843.

A tryptic digest of a Saccharomyces cerevisiae lysate was used as background, identical in all samples. Into it six proteins were spiked at four different combinations of concentrations, giving four sample groups, each analysed in triplicate on a Q-Exactive instrument: twelve LC-MS/MS runs in total. Table 1 of Choi et al. lists the amounts, in fmol per injection:

None of the six spiked proteins is a yeast protein. In the data set they carry disguised, yeast-like names and accession numbers, so that the participants in the study could not spot them by looking at the protein list. Every other protein in the experiment is a yeast protein and should, by construction, be unchanged.

We start from the protein-level table that the MSstats pipeline produces after summarising the Skyline peptide intensities into one log$_2$ intensity per protein and run.

The twelve columns are named after the sample group and the technical replicate, which makes the design easy to recover.

The design is visible in the data¶

Before doing any statistics we look at the quantities themselves. For every protein we subtract its own mean log$_2$ intensity, putting all proteins on a comparable scale, and plot the profile over the four conditions.

The grey band of background proteins is essentially flat — 95% of them vary by less than a factor of two between the extreme conditions — while the six coloured profiles swing by factors of forty to five hundred, following the concentration table above.

Differential abundance between two conditions¶

We now perform the analysis we would have performed without knowing the truth: a $t$ test of every protein between condition 1 and condition 2, as in the differential expression chapter. One detail differs. With only three replicates per group we use the pooled-variance $t$ test rather than Welch’s, since Welch’s correction buys robustness against unequal variances at a price in degrees of freedom that we cannot afford at $n=3$.

The $p$ value histogram of the background proteins is the diagnostic that the null distribution behaves. If the yeast proteins really are unchanged, and the $t$ test is correctly calibrated, their $p$ values should be uniformly distributed on $[0,1]$.

The histogram is close to flat, if somewhat depleted near zero, and the fraction of background proteins below $p=0.05$ is a little under the 5% that a uniform null would give. The background is a proper — indeed mildly conservative — null. Five of the six spiked proteins sit at the extreme left edge. The sixth is VAC2, which was spiked at 65 and 55 fmol in the two conditions — a genuine but only 1.2-fold difference that three replicates cannot resolve. It will end up as a false negative, but not one that any method could reasonably be expected to avoid.

$q$ values¶

We estimate $\pi_0$ and the $q$ values with the same bootstrap procedure as in the previous chapter, here imported from the book’s common module rather than retyped.

The estimate $\hat{\pi}_0\approx 1$ says that essentially none of the proteins are differentially abundant, which is the correct answer: 6 out of 3025 is $0.2\%$. Contrast this with the comparison of the two mouse strains in the previous chapter, where $\hat{\pi}_0$ came out at about 0.5.

The payoff: estimated FDR against the true FDP¶

Because we know which proteins were spiked, we can walk down the list in order of increasing $p$ value and, at every position, compute both what the procedure claims the error rate is and what it actually is.

The four highest ranked proteins are all spiked, and there the $q$ value and the FDP agree at essentially zero. The first background protein enters at rank five, where the observed FDP jumps to $1/5=0.2$ against an estimate of 0.04. From then on the two curves track each other loosely, with the estimate rising the faster of the two: by the tenth acceptance it has already saturated at one, while the true FDP is only one half. Note also how granular the FDP is: among $k$ accepted proteins it can take only the values $0, 1/k, 2/k, \ldots$, and with only six true positives $k$ stays small in the interesting part of the list, so one unlucky background protein at rank five moves it by twenty percentage points. That granularity is a property of the benchmark, not of the method.

To get a less noisy picture we repeat the analysis for all six pairs of conditions and record, at a set of $q$ value thresholds, how many proteins were accepted and what fraction of them were background.

Points on the dashed line are perfectly calibrated, points below it conservative and points above it anti-conservative. At the strictest settings the procedure makes almost no mistakes: five of the six comparisons have an FDP of exactly zero at $q\le0.02$, and four of them still do at $q\le0.1$. Further out the points scatter around and mostly above the line, with FDPs between 0.14 and 0.33 at $q\le0.2$. Part of that drift is unavoidable arithmetic rather than a real bias — one false discovery among five accepted proteins is already an FDP of 0.2, so a benchmark with six true positives simply cannot resolve an FDR of 0.05.

Condition 1 versus 3 is the real outlier. At $q\le0.1$ roughly two thirds of the accepted proteins are background, and its $\hat{\pi}_0$ comes out around 0.6 — the estimator concluded that more than a third of the yeast proteins were differentially abundant. Printing $\pi_0$ and the behaviour of the background for every comparison shows where the problem lies.

This is not a failure of the $q$ value machinery but of our assumption about the ground truth. For condition 1 versus 3 more than one in ten of the yeast proteins falls below $p=0.05$, about twice what a uniform null would give, and the estimator is right to conclude that something is changing. The yeast background is simply not perfectly identical between those two sample groups. All four samples were made from the same yeast digest and the twelve runs were acquired in random order, but each sample was spiked separately, and with only three runs per sample, randomisation cannot fully average out drift of the instrument. Small differences in the amount of digest that ended up in each sample, or in the conditions under which its runs were acquired, leave systematic traces in the quantities. No FDR estimator can protect against a null that is not a null, and the proteins it flags are not errors of statistics but of our labelling of the truth. A benchmark tells us about the whole pipeline — sample preparation, quantification, normalisation and statistics together — not only about the last step.