2 Analysis with 8 replicate samples from an experiment
Load required packages.
library("Biobase")
library("limma")
library("ggplot2")
library("beyonce")
library("magrittr")
library("reshape2")Load an example microarray gene expression dataset, the ApoAI data from the limma user’s guide.
if (!exists("MA")) {
zipfile <- "apoai.zip"
download.file("https://www.bioconductor.org/help/course-materials/2005/BioC2005/labs/lab01/Data/apoai.zip",
destfile = zipfile)
unzip(zipfile)
load("ApoAI.RData")
MA <- normalizeWithinArrays(RG)
}The dataset MA comprises 16 samples, which come in two groups:
table(MA$targets$Cy5)##
## ApoAI-/- C57BL/6
## 8 8Let’s focus on the wildtype samples (controls).
MAwt <- MA[, MA$targets$Cy5=="C57BL/6"]We use the function combn to generate all possible subdivisions of these 8 samples into two groups of equal size.
splits <- combn(ncol(MAwt), ncol(MAwt) %/% 2)
dim(splits)## [1] 4 70splits[, 1:12]## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12]
## [1,] 1 1 1 1 1 1 1 1 1 1 1 1
## [2,] 2 2 2 2 2 2 2 2 2 2 2 2
## [3,] 3 3 3 3 3 4 4 4 4 5 5 5
## [4,] 4 5 6 7 8 5 6 7 8 6 7 8The following code loops over the columns of splits and calls limma’s lmFit/eBayes differential expression analysis on each of the 70 subdivisions. Since the samples are all replicates, all null hypotheses are true, i.e. there is no true true differential expression. The resulting p-values are obtained in the matrix pvalues.
pvalues <- apply(splits, 2, function(i) {
mock <- seq_len(ncol(MAwt)) %in% i
fit <- lmFit(MAwt$M, design = model.matrix(~ mock)) %>% eBayes
fit$p.value[, "mockTRUE"]
})
str(pvalues)## num [1:6384, 1:70] 0.4365 0.2706 0.0474 0.1847 0.4188 ...Let’s look at the histogram of the p-values. For ggplot2 it is convenient to stick the data into a data.frame df. We also define 4 equal-sized groups of genes, according to their overall variance, in the column strat_var, and similarly, 4 equal-sized groups according to quartiles of the overall mean.
df <- melt(pvalues, varnames = c("probe", "replicate"), value.name = "p")
means <- base::rowMeans(MAwt$M) %>% rank %>% cut(4, labels = paste("mean quartile", 1:4))
vars <- genefilter::rowVars(MAwt$M) %>% rank %>% cut(4, labels = paste("variance quartile", 1:4))
df$strat_mean <- means[df$probe]
df$strat_var <- vars[df$probe]h <- ggplot(df , aes(x = p, ..density..)) + geom_histogram(binwidth = 0.025, fill = beyonce_palette(9)[3]) +
geom_hline(yintercept = 1)
h
If all null hypotheses are true, the distribution of each p-value should be uniform. Indeed, we see that the mixture distribution of the p-values for all genes is pretty much uniform. Now let’s stratify by overall variance, i.e. the factor strat_var that we defined above.
h + facet_wrap( ~ strat_var, ncol = 2)
The histograms -which correspond to the conditional distributions of the p-values, given the stratification variable- are no longer uniform. We see that limma tends to be conservative for genes with small overall variance, and anticonservative for genes with large overall variance. This is, of course, as intended by design of the method, due to the shrinkage estimation of the variance: genes that have a small overall (empirical) variance also tend to have a small empirical within-groups variance, which then gets shrunk upward for the purpose of the moderated t-test. And vice versa for genes with a large within-groups variance.
In contrast, stratifying (or conditioning) by mean leads to only weak (and probably negligible) deviations from uniformity.
h + facet_wrap( ~ strat_mean, ncol = 2)
