Here I will introduce a popular statistical framework for working with the results of many parallel statistical tests, a setting that occurs frequently in genomic analyses. I won’t go into depth here, but will discuss more in lecture. For background you can also read the lecture notes from the HarvardX Biological Data Science series (also the Wikipedia article on multiple comparisons is a good introduction).
While up to this point, we have minimized the use of simulation data to show concepts in Computational Biology, multiple testing is often introduced using simulation data, and there is a practical reason why. We can use real biological datasets to generate test statistics from the null hypothesis, by taking a large group of samples and randomly dividing into two groups with no known biological difference. Unless we happen to choose a partition which corresponds to biological or technical factors, this is pretty good for observing null test statistics. And we can limit this probability by choosing a large dataset, or by balancing our random partition with respect to technical factors, or by removing technical factors using factor analysis methods such as SVA, RUV, or PEER.
However, it’s not easy to generate many test statistics from real data which combine known null and alternative hypotheses, because we often do not know which hypotheses are truly null, e.g. \(\Delta = 0\), and which are alternative, \(\Delta \ne 0\). So we often resort to simulation, or a combination of real null data with simulation used to create alternative hypotheses, to demonstrate the concepts of multiple test correction. Here we will use the latter, starting with a real dataset where we can create a null distribution of test statistics.
About the dataset
We use a publicly available gene expression dataset from the HipSci project, accessible via the recount3 Bioconductor package. The study (accession ERP020977) profiled gene expression in macrophages derived from induced pluripotent stem cells (iPSCs) of 86 unrelated healthy donors under four conditions: naive, interferon-gamma (IFN-γ), Salmonella infection, and IFN-γ followed by Salmonella. The publication is:
Alasoo et al. (2018) “Shared genetic effects on chromatin and gene expression indicate a role for enhancer priming in immune response.” Nature Genetics.
For our null distribution, we focus on the naive (untreated) samples, where no biological difference between random groups of donors is expected.
We parse condition and donor from the sample attributes and title, grab the fragment count for later use with voom, then drop the bulk SRA and project metadata columns from colData:
library(stringr)library(tidyr)attrs <-colData(rse)$sra.sample_attributescolData(rse)$frag_count <-colData(rse)$`recount_qc.bc_frag.count`colData(rse)$donor <-str_extract(colData(rse)$sra.sample_title, "HPSI[^ ]+")# Extract the treatment field; combined IFNg+Salmonella samples only have# "interferon gamma" there — their Sample Name sub-attribute is the only# indicator, so we detect them in a second passcolData(rse)$treatment <-str_extract(attrs, "(?<=treatment;;)[^|]+") |>replace_na("naive")colData(rse)$treatment[str_detect(attrs, "salmonella and IFNg")] <-"IFNg_Salmonella"# Drop the heavy SRA and project metadata now that we have what we needcd <-colData(rse)colData(rse) <- cd[, !grepl("^sra\\.|^recount_project\\.", colnames(cd))]table(rse$treatment)
We normalize the raw coverage counts using transform_counts, which scales each sample to 40 million mapped reads (analogous to RPM normalization). We keep a raw count matrix for use with limma::voom, then filter to genes expressed in more than 70% of samples.
I take a quick look at the histogram of values for the first four genes:
par(mfrow =c(2, 2))for (i in1:4) hist(mat[i, ], col ="grey", main =rownames(mat)[i])
Also, it’s always a good idea to make a PCA plot of the samples, to look for any large-scale structure in the data:
par(mfrow =c(1, 1))pc <-prcomp(t(mat))plot(pc$x[, 1:2], main ="PCA (unlabeled)")
We can see that the first PC captures substantial variation. Coloring by donor reveals that this structure is driven by individual-to-individual differences — each donor was profiled at a single time point, so there is no within-donor replication across conditions in the naive subset.
plot(pc$x[, 1:2], col =as.integer(rse$donor),pch =20, main ="PCA colored by donor")
If you square the standard deviation values and sum them, you get the total variance in the data. The top PC explains a notable fraction of variance, consistent with the donor structure we see in the plot.
plot(pc$sdev[1:10]^2/sum(pc$sdev^2), type ="b", ylab ="% Var")
Multiple testing
To investigate multiple testing, we first generate a mock comparison: I choose two random groups of the data, indicated by a factor variable called fake. Rather than a simple row-wise t-test (which would be inflated by the donor structure visible in the PCA), we use limma::voom to fit a linear model that includes donor as a covariate. This removes the donor-to-donor variance from the residuals, giving a much more uniform null distribution of p-values.
library(limma)set.seed(5)n <-ncol(mat_raw)fake <-factor(sample(c(rep(1, n %/%2), rep(2, n - n %/%2))))design <-model.matrix(~ rse$donor + fake)v <-voom(mat_raw, design)fit <-lmFit(v, design)fit <-eBayes(fit)res <-topTable(fit, coef ="fake2", number =Inf, sort.by ="none")nbins <-20brks <-0:nbins / nbinshist(res$P.Value, col ="grey", breaks = brks,main ="Null p-values (donor-corrected voom)")
Here I make this plot again, and show a line for what we would expect on average for each bin for null hypotheses. We can see we are pretty close to the expected values for this set of bins.
p <- res$P.Valuehist(p, col ="grey", breaks = brks)abline(h =length(p) / nbins, col ="dodgerblue", lwd =3)
1/ nbins
[1] 0.05
sum(p <1/ nbins)
[1] 874
length(p) / nbins
[1] 897.5
Now, because we want to explore how multiple testing frameworks deal with combined null and alternative hypotheses, I will “spike in” 1,000 p-values corresponding to alternative hypotheses. I don’t generate the data, just p-values which, after log10 transformation, are uniformly distributed on [-6,-1].
p2 <- pp2[1:1000] <-10^runif(2000, -6, -1)
Warning in p2[1:1000] <- 10^runif(2000, -6, -1): number of items to replace is not a multiple of
replacement length
Correction with p.adjust
I demonstrate two types of correction for multiple testing, the Bonferroni method and the Benjamini-Hochberg method, which are both implemented in the p.adjust function in R. The p.adjust function provides adjusted p-values, such that thresholding on the adjusted p-values delivers a set obeying certain statistics properties. As we went over in lecture, the Bonferroni method bounds the family-wise error rate (FWER) while the Benjamini-Hochberg method bounds the false discovery rate (FDR) in expectation. A reminder: these are critically different bounds: FWER is the number of false positives over all null hypotheses, while FDR is the number of false positives over the set which is called positive. And the Benjamini-Hochberg method provides the FDR bound in expectation.
Below we show the histogram of adjusted p-values, and the minimal adjusted p-value. The adjustment is simply the original p-value multiplied by the number of tests. So we could provide FWER control for this particular simulated data for a few hundred tests, but only at a high rate.
padj.bonf <-p.adjust(p2, method ="bonferroni")hist(padj.bonf, col ="grey", breaks = brks)
min(padj.bonf)
[1] 0.01799849
min(p2) *length(p2) ==min(padj.bonf)
[1] TRUE
sum(padj.bonf < .5)
[1] 272
The false discovery rate is a much more practical and desirable bound for genomic data analysis, compared to the family-wise error rate. Investigators are often interested in the rate of false positives among the set of most promising features (genes, or in this case, genomic regions) that are identified by a statistical test.
Here, we pick up on many tests where we can control the FDR at 5%, for example.
padj.bh <-p.adjust(p2, method ="BH")hist(padj.bh, col ="grey", breaks = brks)
min(padj.bh)
[1] 0.0005184323
sum(padj.bh < .05)
[1] 670
One way to think about what our adjusted p-values are delivering, is to examine a histogram of the original p-values for various bin sizes along [0,1]. We see an enrichment of small p-values in the smallest bin. If we draw the line of expected counts if all the hypotheses were null, we see that the first bin could be expected to contain more than half null hypotheses. This corresponds to the largest adjusted p-value in this bin, which makes sense: by thresholding on an adjusted p-value, we should obtain a set of hypotheses bounded by a given false discovery rate in expectation.
nbins <-20brks <-0:nbins / nbinshist(p2, col ="grey", breaks = brks)abline(h =length(p) / nbins, col ="dodgerblue", lwd =3)
max(padj.bh[p2 <1/ nbins])
[1] 0.5136901
If we make the bins smaller, the ratio of expected nulls decreases relative to the height of the bar for the first bin, and again the maximum adjusted p-value matches the expected proportion in that bin.
nbins <-100brks <-0:nbins / nbinshist(p2, col ="grey", breaks = brks)abline(h =length(p) / nbins, col ="dodgerblue", lwd =3)
max(padj.bh[p2 <1/ nbins])
[1] 0.1863315
Runs of identical adjusted p-values
Finally, I want to show a property of adjusted p-values using the Benjamini-Hochberg (BH) method, which sometimes surprises users who aren’t familiar with the procedure. See the runs of identical adjusted p-values for the first 500 sorted p-values:
The publication of the BH method describes a procedure, whereby one finds the largest i such that the i-th smallest p-value is less than \(\frac{i}{m} q\), where q is the desired FDR bound. The adjusted p-values are the smallest value of q for each test, so a reverse of the procedure where you start with q. Geometrically, the BH procedure can be thought of putting the tests in order of p-value along the x-axis, with the value of p on the y-axis. For the first set of tests, this looks like:
We then can find the smallest q to define a set, by drawing the line \(y(i) = \frac{q}{m} i\). So we draw a line that goes through the origin with increasing slope. This first touches our sorted values of p at an asymptotic point, and so all the p-values less than this point are assigned the same adjusted p-value. They could only have a smaller adjusted p-value if the line had touches the sorted p-values earlier.
q <-min(padj.bh)i <-sum(padj.bh == q)m <-length(p2)plot(1:n_p, p.sort[1:n_p], ylim =c(0, p.sort[n_p]), type ="l")abline(0, q / m, col ="red")
q-values
Finally, we mention that there is a method which improves upon the BH method by attempting to estimate the proportion of nulls in the set of all hypotheses, called \(\pi_0 = \frac{m_0}{m}\). The BH method controls the FDR in expectation for any \(m_0 \le m\), but if we can come up with a good estimate of \(\pi_0\) when it is less than 1, we can reject more hypotheses and still control the FDR in expectation. Without getting into details, the qvalue function in the qvalue package does this estimation, and provides q-values which operate analogously to the adjusted p-values we defined above: using a threshold on q-values produces a set with FDR which should be bounded by that value in expectation. We can see that for this simulated data, the estimate of \(\pi_0\) is close to the truth:
library(qvalue)qres <-qvalue(p2)qres$pi0
[1] 0.9717311
1-1000/length(p) # we spiked in 1k alternatives
[1] 0.9442897
The q-values are multiplicatively scaled relative to the BH adjusted p-values, and are strictly smaller, although here only slightly so. For smaller \(\pi_0\) the gain in sensitivity would be more substantial.
