In stephaniehicks/benchmarkfdrData2019: Data and Benchmarking Results from Korthauer and Kimes et al. (2019)
knitr::opts_chunk$set(echo = TRUE)
Introduction
This is an analysis of the yeast data with 48 biological replicates in each of
two conditions (analyzed in this publication).
We chose this experiment because of the large number of biological replicates, which
will allow us to (1) implement null comparisons on random subsets of samples within
one condition, and (2) start with a null comparison and add in artificial differences
to a subset of genes to define 'true positives'.
We will simulate a strongly informative and weakly informative independent covariate.
We will also investigate the effects of an uninformative covariate by using a
randomly generated covariate.
Processed count table is made available by the authors in their paper codebase
GitHub repository.
In this Rmd we will carry out simulations for multiple
replicates, and plot results averaging over the replications. In contrast to
yeast-simulation.Rmd we will use an alternate simulation framework for generating
the effect sizes of the DE genes. Here we will use the fold changes of only genes
that were DE in a comparison of wild type versus knockout (with sample size 10).
Since this results in a distribution of effect sizes (log2 fold changes) under the
alternative that does not have a mode at zero, the assumptions for ashq are violated.
However, since this assumption is impossible to check, we still include this
method in the benchmark comparisons.
In this setting, we assume 500 DE genes and a bimodal alternative.
Set up workspace
# Load packages and source benchmark FDR
library(SummarizedBenchmark)
library(data.table)
library(tidyr)
library(dplyr)
library(readr)
library(ggplot2)
library(magrittr)
library(cowplot)
library(purrr)
library(DESeq2)
library(tibble)
library(ggthemes)
library(R.utils)
library(Hotelling)
## load helper functions
for (f in list.files("../R", "\\.(r|R)$", full.names = TRUE)) {
source(f)
}
# set up data / results directories
datdir <- "yeast-data"
resdir <- "../../results/RNAseq"
dir.create(datdir, showWarnings = FALSE)
dir.create(resdir, showWarnings = FALSE)
# results files that will be generated
# strong covariate
resfile_d5 <- file.path(resdir, "yeastIIH-results-de5.rds")
resfile_d10 <- file.path(resdir, "yeastIIH-results-de10.rds")
# weak covariate
resfile_d5_w <- file.path(resdir, "yeastIIHW-results-de5.rds")
resfile_d10_w <- file.path(resdir, "yeastIIHW-results-de10.rds")
# uninformative covariate
resfile_d5_uninfCov <- file.path(resdir, "yeastIIH-results-de5-uninfCov.rds")
resfile_d10_uninfCov <- file.path(resdir, "yeastIIH-results-de10-uninfCov.rds")
# set up parallel backend
library(parallel)
nCores <- 20
Data download
First, we download the processed count data from GitHub. There is one file for the
Snf2 condition and one file for the wild type condition. The Snf2 condition is a yeast
strain that has the transcriptional regulator gene Snf2 knocked out. Each of these
files is a compressed tar.gz archive that contains a
single bam file for each replicate in each condition, which we'll place in a subdirectory
called r datdir.
download.file(url = "https://github.com/bartongroup/profDGE48/raw/master/Preprocessed_data/Snf2_countdata.tar.gz",
destfile = file.path(datdir, "Snf2_countdata.tar.gz"))
download.file(url = "https://github.com/bartongroup/profDGE48/raw/master/Preprocessed_data/WT_countdata.tar.gz",
destfile = file.path(datdir, "WT_countdata.tar.gz"))
gunzip(file.path(datdir, "Snf2_countdata.tar.gz"))
gunzip(file.path(datdir, "WT_countdata.tar.gz"))
untar(file.path(datdir, "Snf2_countdata.tar"), exdir = datdir)
untar(file.path(datdir, "WT_countdata.tar"), exdir = datdir)
file.remove(file.path(datdir, "Snf2_countdata.tar"),
file.path(datdir, "WT_countdata.tar"))
Each of the data files contains two columns, one with a gene/feature name, and one
with the count value.
We'll also download the list of 'bad' replicates which were especially poorly
correlated to the others, as determined by the authors.
download.file(url = "https://github.com/bartongroup/profDGE48/raw/master/Bad_replicate_identification/exclude.lst", destfile = (file.path(datdir, "badreps.txt")))
Read data into R and create a count table
Here we make use of the map and map2 functions in the purrr package, to swiftly
apply the read_tsv function from readr to read in all of the 96 sample tables,
as well as add in the sample name (derived from the file name) to each subtable.
Finally, the reduce function is used to join all the replicates together.
We will remove samples that failed QC in the original study.
files <- dir(path = datdir, pattern = "*.bam.gbgout", full.names = TRUE)
sample_names <- sapply(strsplit(dir(path = datdir, pattern = "*.bam.gbgout"), "_MID"),
function(x) x[[1]])
badreps <- read_tsv(file.path(datdir, "badreps.txt"), col_names = FALSE)$X1
badreps <- unlist(lapply(strsplit(badreps, "_MID"), function(x) x[1]))
counts <- files %>%
purrr::map(read_tsv, col_names = FALSE) %>% # read in all the files individually
purrr::map2(sample_names, ~ dplyr:::rename(.x, !! .y := X2, feature = X1) ) %>% # add sample names
purrr::reduce(left_join, by = "feature") %>% # reduce with rbind into one dataframe
dplyr::select(-badreps ) # remove badreps
Empirical analysis of 10 vs 10 (non-null)
Here we'll carry out an analysis of a subset of the full dataset, comparing the controls to the
knockout samples. The fold changes observed in this comparison will be used when
generating the non-null simulated data in the following sections. Here we use a subset
of the samples (10 samples in each group) to mimic the sample sizes we use in
simulation.
We're ready to construct a DESeq2 object. First we pull out the feature names and
add them as rownames for the count table, and next we construct a column data object
that houses the sample names, replicate numbers, and condition factor (WT versus
Snf2 knockout).
Set up DESeq2 object
First we set up the DESeq2 object.
feats <- (counts %>% select(1))$feature
counts <- as.matrix(counts %>% select(-1))
rownames(counts) <- feats
coldat <- tibble(sample=colnames(counts)) %>%
separate(sample, sep="_", into=c("condition", "replicate"), remove=FALSE) %>%
mutate(condition = factor(condition))
# filter low count genes
counts <- counts[rowMeans(counts) > 1,]
dds_full <- DESeqDataSetFromMatrix(countData = counts,
colData = coldat,
design= ~ condition)
Run DESeq2
Next we run DESeq2, and subset results on genes with FDR < 0.05.
# results on full set
set.seed(9384)
subset.wt <- dds_full[,colData(dds_full)$condition == "WT"]
subset.sn <- dds_full[,colData(dds_full)$condition == "Snf2"]
subset <- cbind(subset.wt[, sample(1:ncol(subset.wt), 10)],
subset.sn[, sample(1:ncol(subset.sn), 10)])
dds <- DESeq(subset)
resultsNames(dds) # lists the coefficients
res_10 <- results(dds, name="condition_WT_vs_Snf2", independentFiltering = FALSE)
head(res_10)
sum(res_10$padj < 0.05, na.rm = TRUE)
Check densities of the test statistic and effect sizes for significant DE genes.
data.frame(res_10[res_10$padj < 0.05 & !is.na(res_10$padj),]) %>%
ggplot(aes(x = stat)) +
geom_density() +
xlab("Test Statistic")
data.frame(res_10[res_10$padj < 0.05 & !is.na(res_10$padj),]) %>%
ggplot(aes(x = log2FoldChange)) +
geom_density() +
xlab("Effect Size (log2 Fold Change)")
Check assumptions
There are several thousand genes detected using the full set. Next, we'll build input data frame for summarized benchmark.
geneExp <- tbl_df(data.frame(geneName=rownames(res_10),
pval=res_10$pvalue,
SE=res_10$lfcSE,
ind_covariate = res_10$baseMean,
effect_size=res_10$log2FoldChange,
test_statistic=res_10$stat,
pzero=rowSums(counts(subset)==0)/ncol(counts(subset))))
# filter NAs and those with less than 50% expressed
geneExp <- geneExp %>% na.omit() %>% dplyr::filter(pzero < 0.5)
We'll create a plot to examine the distribution of effect sizes, since
the ash method assumes that the distribution of true (unobserved) effect
sizes is unimodal.
ggplot(data=geneExp, aes(effect_size)) +
geom_histogram(bins=30)
We'll also explore how the standard error (used by ash)
correlates with the independent covariate (used by methods that incorporate
covariates), in order to get an idea of how these pieces of information relate
to one another.
ggplot(data=geneExp, aes(x = ind_covariate, y = SE)) +
geom_hex(bins = 100) +
scale_x_continuous(trans="log10") +
xlab("Covariate: Mean gene expression")
Look at covariate diagnostic plots.
strat_hist(geneExp, pvalue="pval", covariate="ind_covariate", maxy=30)
rank_scatter(geneExp, pvalue="pval", covariate="ind_covariate")
FDR benchmarking
Build common bench design. We also add in Scott's FDR Regression
since our test statistics are approximately t-distributed.
Note that the assumption of FDRreg-empirical that the distribution of non-null test
statistics does not have a significant mass at zero is violated here, but
we still include fdrreg-e in the evaluation since this assumption is impossible to check
outside of a simulation setting (since it is not possible to know which tests are
non-null). Thus, we'd like to evaluate how results change
when the assumption is violated, since it's plausible that it is violated in the
case studies.
bd <- initializeBenchDesign()
bd <- addBMethod(bd, "fdrreg-t",
FDRreg::FDRreg,
function(x) { x$FDR },
z = test_statistic,
features = model.matrix( ~ splines::bs(ind_covariate, df = 3) - 1),
nulltype = 'theoretical',
control = list(lambda = 0.01))
bd <- addBMethod(bd, "fdrreg-e",
FDRreg::FDRreg,
function(x) { x$FDR },
z = test_statistic,
features = model.matrix( ~ splines::bs(ind_covariate, df = 3) - 1),
nulltype = 'empirical',
control = list(lambda = 0.01))
Run benchmark methods.
sb <- bd %>% buildBench(data=geneExp, parallel = FALSE)
assayNames(sb) <- "qvalue"
sb <- addDefaultMetrics(sb)
Plot results.
rejections_scatter(sb, supplementary=FALSE)
plotFDRMethodsOverlap(sb, alpha=0.05, nsets=ncol(sb), order.by="freq", decreasing=TRUE, supplementary=FALSE)
Empirical analysis of 5 vs 5 (non-null)
Here we'll carry out an analysis of a subset of the full dataset, comparing the controls to the
knockout samples. The fold changes observed in this comparison will be used when
generating the non-null simulated data in the following sections. Here we use a subset
of the samples (5 samples in each group) to mimic the sample sizes we use in
simulation.
Run DESeq2
First we run DESeq2, and subset results on genes with FDR < 0.05.
# results on full set
set.seed(728)
subset.wt <- dds_full[,colData(dds_full)$condition == "WT"]
subset.sn <- dds_full[,colData(dds_full)$condition == "Snf2"]
subset <- cbind(subset.wt[, sample(1:ncol(subset.wt), 5)],
subset.sn[, sample(1:ncol(subset.sn), 5)])
dds <- DESeq(subset)
resultsNames(dds) # lists the coefficients
res_5 <- results(dds, name="condition_WT_vs_Snf2", independentFiltering = FALSE)
head(res_5)
sum(res_5$padj < 0.05, na.rm = TRUE)
Check densities of the test statistic and effect sizes for significant DE genes.
data.frame(res_5[res_5$padj < 0.05 & !is.na(res_5$padj),]) %>%
ggplot(aes(x = stat)) +
geom_density() +
xlab("Test Statistic")
data.frame(res_5[res_5$padj < 0.05 & !is.na(res_5$padj),]) %>%
ggplot(aes(x = log2FoldChange)) +
geom_density() +
xlab("Effect Size (log2 Fold Change)")
Check assumptions
There are several thousand genes detected using the full set. Next, we'll build input data frame for summarized benchmark.
geneExp <- tbl_df(data.frame(geneName=rownames(res_5),
pval=res_5$pvalue,
SE=res_5$lfcSE,
ind_covariate = res_5$baseMean,
effect_size=res_5$log2FoldChange,
test_statistic=res_5$stat,
pzero=rowSums(counts(subset)==0)/ncol(counts(subset))))
# filter NAs and those with less than 50% expressed
geneExp <- geneExp %>% na.omit() %>% dplyr::filter(pzero < 0.5)
We'll create a plot to examine the distribution of effect sizes, since
the ash method assumes that the distribution of true (unobserved) effect
sizes is unimodal.
ggplot(data=geneExp, aes(effect_size)) +
geom_histogram(bins=30)
We'll also explore how the standard error (used by ash)
correlates with the independent covariate (used by methods that incorporate
covariates), in order to get an idea of how these pieces of information relate
to one another.
ggplot(data=geneExp, aes(x = ind_covariate, y = SE)) +
geom_hex(bins = 100) +
scale_x_continuous(trans="log10") +
xlab("Covariate: Mean gene expression")
Look at covariate diagnostic plots.
strat_hist(geneExp, pvalue="pval", covariate="ind_covariate", maxy=23)
rank_scatter(geneExp, pvalue="pval", covariate="ind_covariate")
FDR benchmarking
Run benchmark methods.
sb <- bd %>% buildBench(data=geneExp, parallel = FALSE)
assayNames(sb) <- "qvalue"
sb <- addDefaultMetrics(sb)
Plot results.
rejections_scatter(sb, supplementary=FALSE)
plotFDRMethodsOverlap(sb, alpha=0.05, nsets=ncol(sb), order.by="freq", decreasing=TRUE, supplementary=FALSE)
Simulation set up
First we'll illustrate the logistic curve used to simulate the informative covariate.
In this study we'll sample from a logistic curve and use these values as probability
weights when selecting DE genes. A weakly informative covariate will be generated by
adding noise to the strongly informative covariate.
x <- seq(0,10, by = 0.1)
infcov <- 1 / (1+exp(-x + 5))
plot(x, infcov, type = "l", ylab = "informative covariate")
plot(infcov, pmin(1, abs(infcov + rnorm(length(x), 0, 0.25))),
ylab = "weak covariate", xlab = "true covariate")
Next, here's a function that will generate sampling weights for test statistics
of DE genes that downweights observations near zero, as well as extremely large values.
stat_wt <- function(z){
pos <- pmax(0,dgamma(z,4.5, 1-1e-4))
neg <- pmax(0,dgamma(-z,4.5, 1-1e-4))
res <- pos
res[z < 0] <- neg[z < 0]
return(res)
}
# Plot weights by test statistic value
df = data.frame(z = seq(-10,10, by = 0.1))
df$w = stat_wt(df$z)
ggplot(data=df, aes(x=z, y=w)) +
geom_line() +
ylab("Sampling Weight") +
xlab("Test statistic value")
Next, we'll analyze random splits of one condition, both with and without the
addition of simulated DE genes. Here we'll create a function that we can use to
run one replicate given sample size and number of DE gene settings. This will be
looped over many replications and results averaged over them.
#' @param X simnumber
#' @param seed random seed
#' @param sampleSize is the number of samples in each condition
#' @param nDE is the number of DE genes
#' @param bd is the bench design object
#' @param BPPARAM is the BiocParallel bpparam argument to pass to DESeq. Since
#' parallelization happens by forking multiple instances of this function, we
#' will specify only one worker per DEseq intance.
#' @param uninformativeCovariate logical indicating whether to replace the
#' informative covariate with a randomly generated (uninformative) covariate.
#' Default is FALSE.
#' @param pvalHists logical whether to return histograms of pvalues by covariate
#' instead of SB objects
simulateOneSplit <- function(X, rseed, nDE, sampleSize,
bd, uninformativeCovariate = FALSE, pvalHists = FALSE,
strongCovariate = TRUE){
# set random seed
set.seed(as.numeric(X)*as.numeric(rseed))
covar <- 1 / (1+exp(-runif(nrow(dds_full), 0, 10) + 5))
# select a random subset of 20 WT samples
dds_test <- dds_full[,sample(1:ncol(dds_full), sampleSize*2)]
# add a fake condition column to coldat
colData(dds_test)$fake <- factor(c(rep("A", sampleSize),
rep("B", sampleSize))[sample(1:(sampleSize*2),
sampleSize*2)])
design(dds_test) <- ~fake
pzero = rowSums(counts(dds_test)==0)/ncol(counts(dds_test))
dds_test <- dds_test[pzero < 0.5,]
covar <- covar[pzero < 0.5]
truth <- rep(FALSE, nrow(dds_test))
# make sure null comparison is truly null: if PC1 or PC2 sig different,
# or if test of join means of PCs 1:4 is significant,
# reshuffle sample labels
dds_test <- estimateSizeFactors(dds_test)
x <- t(counts(dds_test, normalize=TRUE))
pc <- prcomp(log(x + 0.5), scale.=TRUE)
a1 <- pc$x[colData(dds_test)$fake=="A",1]
b1 <- pc$x[colData(dds_test)$fake=="B",1]
a2 <- pc$x[colData(dds_test)$fake=="A",2]
b2 <- pc$x[colData(dds_test)$fake=="B",2]
p1 <- t.test(a1, b1)$p.value
p2 <- t.test(a2, b2)$p.value
tries <- 0
while(p1 < 0.10 || p2 < 0.10 && tries < 10){
colData(dds_test)$fake <- sample(colData(dds_test)$fake, ncol(dds_test))
x <- t(counts(dds_test, normalize=TRUE))
pc <- prcomp(log(x + 0.5), scale.=TRUE)
a1 <- pc$x[colData(dds_test)$fake=="A",1]
b1 <- pc$x[colData(dds_test)$fake=="B",1]
a2 <- pc$x[colData(dds_test)$fake=="A",2]
b2 <- pc$x[colData(dds_test)$fake=="B",2]
p1 <- t.test(a1, b1)$p.value
p2 <- t.test(a2, b2)$p.value
tries <- tries + 1
}
if(nDE > 0){
# pick random set of nDE genes to add signal to
DE <- sample(1:nrow(dds_test), nDE, prob = covar)
truth[DE] <- TRUE
# randomly sample a test statistic from original test statistics, with probability
# weights that downweight observations near zero
counts_new <- counts(dds_test)
stat <- rep(0, nrow(dds_test))
# assume fixed standard error (median observed) and calculate FC from stat
if(sampleSize == 5){
stat[DE] <- sample(res_5$stat, nDE, prob = stat_wt(res_5$stat))
log2FC <- stat * median(res_5$lfcSE)
}else if (sampleSize == 10){
stat[DE] <- sample(res_10$stat, nDE, prob = stat_wt(res_10$stat))
log2FC <- stat * median(res_10$lfcSE)
}else{
stop("Only sample sizes 5 and 10 are currently supported with pre-",
"computed fold changes for sampling from.")
}
# randomize which condition is shifted up or down
ran <- runif(nrow(dds_test))
refcond <- ifelse(ran < 0.5, "A", "B")
down <- which(ran < 0.5)
counts_new[down,colData(dds_test)$fake==unique(refcond[down])] <-
counts(dds_test)[down, colData(dds_test)$fake==unique(refcond[down])] *
2^log2FC[down]
counts_new[-down,colData(dds_test)$fake==unique(refcond[-down])] <-
counts(dds_test)[-down, colData(dds_test)$fake==unique(refcond[-down])] *
2^log2FC[-down]
counts_new <- apply(counts_new, 2, as.integer)
counts(dds_test) <- counts_new
}
# replace existing size factors
dds_test <- estimateSizeFactors(dds_test)
dds_test <- DESeq(dds_test, parallel = FALSE)
resTEST <- results(dds_test, name="fake_B_vs_A", independentFiltering = FALSE)
geneExp <- tbl_df(data.frame(geneName=rownames(resTEST),
pval=resTEST$pvalue,
SE=resTEST$lfcSE,
ind_covariate = covar,
effect_size = resTEST$log2FoldChange,
test_statistic = resTEST$stat,
qvalue = truth))
if (uninformativeCovariate){
geneExp <- mutate(geneExp, ind_covariate = runif(length(covar)))
}else if(!strongCovariate){
geneExp <- mutate(geneExp, ind_covariate = pmin(1, abs(covar + rnorm(length(covar), 0, 0.25))))
}
geneExp <- geneExp %>% dplyr::filter(!is.na(pval))
if(pvalHists){
return(strat_hist(geneExp, pvalue="pval", covariate="ind_covariate", maxy =10))
}
sb <- bd %>% buildBench(data=geneExp, parallel = FALSE,
truthCols = "qvalue",
ftCols = "ind_covariate")
assayNames(sb) <- "qvalue"
sb <- addDefaultMetrics(sb)
rowData(sb)$log2FC <- geneExp$effect_size
return(sb)
}
We'll also set some parameters that will be common to all simulations. These
include the number of replications, the bench design object, the set of
methods to exclude in the results plots, and the alpha cutoff level to
be used when plotting the aggregated Upset results.
B <- 100
excludeSet <- c("unadjusted", "bl-df02", "bl-df04", "bl-df05")
ualpha <- 0.05
# only keep one condition for subsetting in simulations that follow
dds_full <- dds_full[,colData(dds_full)$condition == "Snf2"]
Here's a helper function to return the number of methods with rejections at
a particular alpha level (this helps us determine whether or not to plot the
aggregated upset plot - if there aren't at least 2 methods it will throw an
error, which is a problem for the null simulations).
# To be included in the upset agg plot, method needs to have found on average
# at least one rejection per replicate. To create an upset plot, require that
# at least two methods rejected at this threshold.
#' @param res standardized metric data.table generated using
#' standardize_results.
#' @param alpha alpha cutoff
#' @param filterSet which methods to exclude from consideration
numberMethodsReject <- function(res, alphacutoff, filterSet){
res <- res %>%
filter(is.na(param.alpha) | (param.alpha == alphacutoff)) %>%
filter(!(blabel %in% filterSet)) %>%
filter(alpha == alphacutoff) %>%
filter(performanceMetric == "rejections") %>%
select(blabel, performanceMetric, value) %>%
group_by(blabel) %>%
summarize(mean_value = mean(value)) %>%
filter(mean_value > 1)
return(nrow(res))
}
Non-null Comparisons with Strongly Informative Covariate
D5S: DE 5v5 Strong
Here we'll repeat the above, but for a DE comparison 5 versus 5 Snf2 samples,
where the groups are selected randomly and 500 DE genes are added.
This will be done for 100 random splits.
Generate a list of SB results objects
rseed <- 198
sampleSize <- 5
nDE <- 500
if (!file.exists(resfile_d5)){
de5 <- mclapply(X=1:B, FUN=simulateOneSplit, rseed=rseed,
nDE=nDE, sampleSize=sampleSize, bd=bd, mc.cores=nCores)
saveRDS(de5, file = resfile_d5)
}else{
de5 <- readRDS(file = resfile_d5)
}
# pvalue histograms
plotfile <- file.path(datdir, "de5_pvalhists.pdf")
if (!file.exists(plotfile)){
hists <- mclapply(X=1:B, FUN=simulateOneSplit, rseed=rseed,
nDE=nDE, sampleSize=sampleSize, bd=bd,
pvalHists = TRUE, mc.cores = nCores)
pdf(plotfile, width=8, height=4)
for(i in 1:length(hists)){
print(hists[[i]])
}
dev.off()
}
Plot average results over replications
Plot results.
# Check for missing results (if any methods threw an error for relevant metrics).
rowSums(sapply(de5, function(x) colSums(is.na(assays(x)$qvalue)) > 0))
res5d <- plotsim_standardize(de5, alpha = seq(0.01, 0.10, 0.01))
plotsim_average(res5d, met="rejections",filter_set = excludeSet,
merge_ihw = TRUE, errorBars=TRUE)
plotsim_average(res5d, met="FDR",filter_set = excludeSet,
merge_ihw = TRUE, errorBars=TRUE)
plotsim_average(res5d, met="TPR",filter_set = excludeSet,
merge_ihw = TRUE, errorBars=TRUE)
plotsim_average(res5d, met=c("FPR", "TPR"), filter_set = excludeSet,
merge_ihw = TRUE)
covariateLinePlot(de5, alpha=0.05, covname="log2FC", nbins=25,
trans="log1p")
covariateLinePlot(de5, alpha=0.05, covname="ind_covariate", nbins=25,
trans="log1p")
if (numberMethodsReject(res5d, alphacutoff=ualpha, filterSet=excludeSet) >= 2){
aggupset(de5, alpha=ualpha, supplementary = FALSE, return_list = FALSE)
}else{
message("Not enough methods found rejections at alpha ", ualpha,
"; skipping upset plot")
}
D10S: DE 10v10 Strong
Here we'll repeat the above, but for a DE comparison 10 versus 10 Snf2 samples,
where the groups are selected randomly and 500 DE genes are added.
This will be done for 100 random splits.
Generate a list of SB results objects
rseed <- 961
sampleSize <- 10
nDE <- 500
if(!file.exists(resfile_d10)){
de10 <- mclapply(X=1:B, FUN=simulateOneSplit, rseed=rseed,
nDE=nDE, sampleSize=sampleSize, bd=bd, mc.cores=nCores)
saveRDS(de10, file = resfile_d10)
}else{
de10 <- readRDS(file = resfile_d10)
}
# pvalue histograms
plotfile <- file.path(datdir, "de10_pvalhists.pdf")
if (!file.exists(plotfile)){
hists <- mclapply(X=1:B, FUN=simulateOneSplit, rseed=rseed,
nDE=nDE, sampleSize=sampleSize, bd=bd,
pvalHists = TRUE, mc.cores = nCores)
pdf(plotfile, width=8, height=4)
for(i in 1:length(hists)){
print(hists[[i]])
}
dev.off()
}
Plot average results over replications
Plot results.
# Check for missing results (if any methods threw an error for relevant metrics).
rowSums(sapply(de10, function(x) colSums(is.na(assays(x)$qvalue)) > 0))
res10d <- plotsim_standardize(de10, alpha = seq(0.01, 0.10, 0.01))
plotsim_average(res10d, met="rejections",filter_set = excludeSet,
merge_ihw = TRUE, errorBars=TRUE)
plotsim_average(res10d, met="FDR",filter_set = excludeSet,
merge_ihw = TRUE, errorBars=TRUE)
plotsim_average(res10d, met="TPR",filter_set = excludeSet,
merge_ihw = TRUE, errorBars=TRUE)
plotsim_average(res10d, met=c("FPR", "TPR"), filter_set = excludeSet,
merge_ihw = TRUE)
covariateLinePlot(de10, alpha=0.05, covname="log2FC", nbins=25,
trans="log1p")
covariateLinePlot(de10, alpha=0.05, covname="ind_covariate", nbins=25,
trans="log1p")
if (numberMethodsReject(res10d, alphacutoff=ualpha, filterSet=excludeSet) >= 2){
aggupset(de10, alpha=ualpha, supplementary = FALSE, return_list = FALSE)
}else{
message("Not enough methods found rejections at alpha ", ualpha,
"; skipping upset plot")
}
Non-null Comparisons with Weakly Informative Covariate
D5W: DE 5v5 Weak
Here we'll repeat the above, but using a weaker informative covariate (that has noise
added to it).
Generate a list of SB results objects
rseed <- 198
sampleSize <- 5
nDE <- 500
if (!file.exists(resfile_d5_w)){
de5 <- mclapply(X=1:B, FUN=simulateOneSplit, rseed=rseed,
nDE=nDE, sampleSize=sampleSize, bd=bd, mc.cores=nCores,
strongCovariate = FALSE)
saveRDS(de5, file = resfile_d5_w)
}else{
de5 <- readRDS(file = resfile_d5_w)
}
Plot average results over replications
Plot results.
# Check for missing results (if any methods threw an error for relevant metrics).
rowSums(sapply(de5, function(x) colSums(is.na(assays(x)$qvalue)) > 0))
res5d <- plotsim_standardize(de5, alpha = seq(0.01, 0.10, 0.01))
plotsim_average(res5d, met="rejections",filter_set = excludeSet,
merge_ihw = TRUE, errorBars=TRUE)
plotsim_average(res5d, met="FDR",filter_set = excludeSet,
merge_ihw = TRUE, errorBars=TRUE)
plotsim_average(res5d, met="TPR",filter_set = excludeSet,
merge_ihw = TRUE, errorBars=TRUE)
plotsim_average(res5d, met=c("FPR", "TPR"), filter_set = excludeSet,
merge_ihw = TRUE)
covariateLinePlot(de5, alpha=0.05, covname="log2FC", nbins=25,
trans="log1p")
covariateLinePlot(de5, alpha=0.05, covname="ind_covariate", nbins=25,
trans="log1p")
if (numberMethodsReject(res5d, alphacutoff=ualpha, filterSet=excludeSet) >= 2){
aggupset(de5, alpha=ualpha, supplementary = FALSE, return_list = FALSE)
}else{
message("Not enough methods found rejections at alpha ", ualpha,
"; skipping upset plot")
}
D10W: DE 10v10 Weak
Here we'll repeat the above, but using a weaker informative covariate (that has noise
added to it).
Generate a list of SB results objects
rseed <- 961
sampleSize <- 10
nDE <- 500
if(!file.exists(resfile_d10_w)){
de10 <- mclapply(X=1:B, FUN=simulateOneSplit, rseed=rseed,
nDE=nDE, sampleSize=sampleSize, bd=bd, mc.cores=nCores,
strongCovariate = FALSE)
saveRDS(de10, file = resfile_d10_w)
}else{
de10 <- readRDS(file = resfile_d10_w)
}
Plot average results over replications
Plot results.
# Check for missing results (if any methods threw an error for relevant metrics).
rowSums(sapply(de10, function(x) colSums(is.na(assays(x)$qvalue)) > 0))
res10d <- plotsim_standardize(de10, alpha = seq(0.01, 0.10, 0.01))
plotsim_average(res10d, met="rejections",filter_set = excludeSet,
merge_ihw = TRUE, errorBars=TRUE)
plotsim_average(res10d, met="FDR",filter_set = excludeSet,
merge_ihw = TRUE, errorBars=TRUE)
plotsim_average(res10d, met="TPR",filter_set = excludeSet,
merge_ihw = TRUE, errorBars=TRUE)
plotsim_average(res10d, met=c("FPR", "TPR"), filter_set = excludeSet,
merge_ihw = TRUE)
covariateLinePlot(de10, alpha=0.05, covname="log2FC", nbins=25,
trans="log1p")
covariateLinePlot(de10, alpha=0.05, covname="ind_covariate", nbins=25,
trans="log1p")
if (numberMethodsReject(res10d, alphacutoff=ualpha, filterSet=excludeSet) >= 2){
aggupset(de10, alpha=ualpha, supplementary = FALSE, return_list = FALSE)
}else{
message("Not enough methods found rejections at alpha ", ualpha,
"; skipping upset plot")
}
Non-null Comparisons with Uninformative Covariate
Here we repeat the previous sections, using an uninformative covariate.
D5U: DE 5v5 Uninformative
Here we'll repeat the above, but for a DE comparison 5 versus 5 Snf2 samples,
where the groups are selected randomly and 500 DE genes are added.
This will be done for 100 random splits.
Generate a list of SB results objects
rseed <- 198
sampleSize <- 5
nDE <- 500
if (!file.exists(resfile_d5_uninfCov)){
de5 <- mclapply(X=1:B, FUN=simulateOneSplit, rseed=rseed,
nDE=nDE, sampleSize=sampleSize, bd=bd,
uninformativeCovariate = TRUE, mc.cores=nCores)
saveRDS(de5, file = resfile_d5_uninfCov)
}else{
de5 <- readRDS(file = resfile_d5_uninfCov)
}
Plot average results over replications
Plot results.
# Check for missing results (if any methods threw an error for relevant metrics).
rowSums(sapply(de5, function(x) colSums(is.na(assays(x)$qvalue)) > 0))
res5d <- plotsim_standardize(de5, alpha = seq(0.01, 0.10, 0.01))
plotsim_average(res5d, met="rejections",filter_set = excludeSet,
merge_ihw = TRUE, errorBars=TRUE)
plotsim_average(res5d, met="FDR",filter_set = excludeSet,
merge_ihw = TRUE, errorBars=TRUE)
plotsim_average(res5d, met="TPR",filter_set = excludeSet,
merge_ihw = TRUE, errorBars=TRUE)
plotsim_average(res5d, met=c("FPR", "TPR"), filter_set = excludeSet,
merge_ihw = TRUE)
covariateLinePlot(de5, alpha=0.05, covname="log2FC", nbins=25,
trans="log1p")
covariateLinePlot(de5, alpha=0.05, covname="ind_covariate", nbins=25,
trans="log1p")
if (numberMethodsReject(res5d, alphacutoff=ualpha, filterSet=excludeSet) >= 2){
aggupset(de5, alpha=ualpha, supplementary = FALSE, return_list = FALSE)
}else{
message("Not enough methods found rejections at alpha ", ualpha,
"; skipping upset plot")
}
D10U: DE 10v10 Uninformative
Here we'll repeat the above, but for a DE comparison 10 versus 10 Snf2 samples,
where the groups are selected randomly and 500 DE genes are added.
This will be done for 100 random splits.
Generate a list of SB results objects
rseed <- 961
sampleSize <- 10
nDE <- 500
if(!file.exists(resfile_d10_uninfCov)){
de10 <- mclapply(X=1:B, FUN=simulateOneSplit, rseed=rseed,
nDE=nDE, sampleSize=sampleSize, bd=bd,
uninformativeCovariate = TRUE, mc.cores=nCores)
saveRDS(de10, file = resfile_d10_uninfCov)
}else{
de10 <- readRDS(file = resfile_d10_uninfCov)
}
Plot average results over replications
Plot results.
# Check for missing results (if any methods threw an error for relevant metrics).
rowSums(sapply(de10, function(x) colSums(is.na(assays(x)$qvalue)) > 0))
res10d <- plotsim_standardize(de10, alpha = seq(0.01, 0.10, 0.01))
plotsim_average(res10d, met="rejections",filter_set = excludeSet,
merge_ihw = TRUE, errorBars=TRUE)
plotsim_average(res10d, met="FDR",filter_set = excludeSet,
merge_ihw = TRUE, errorBars=TRUE)
plotsim_average(res10d, met="TPR",filter_set = excludeSet,
merge_ihw = TRUE, errorBars=TRUE)
plotsim_average(res10d, met=c("FPR", "TPR"), filter_set = excludeSet,
merge_ihw = TRUE)
covariateLinePlot(de10, alpha=0.05, covname="log2FC", nbins=25,
trans="log1p")
covariateLinePlot(de10, alpha=0.05, covname="ind_covariate", nbins=25,
trans="log1p")
if (numberMethodsReject(res10d, alphacutoff=ualpha, filterSet=excludeSet) >= 2){
aggupset(de10, alpha=ualpha, supplementary = FALSE, return_list = FALSE)
}else{
message("Not enough methods found rejections at alpha ", ualpha,
"; skipping upset plot")
}
Simulation setting comparison
Here we compare the method ranks for the different sample sizes, weighting, and
covariate settings at alpha = 0.10.
plotMethodRanks(c(resfile_d5, resfile_d5_w, resfile_d10, resfile_d10_w,
resfile_d5_uninfCov, resfile_d10_uninfCov),
colLabels = c("5 S", "5 W", "10 S", "10 W",
"5 UI", "10 UI"),
alpha = 0.10, xlab = "Comparison",
excludeMethods = NULL)
Session information
sessionInfo()
stephaniehicks/benchmarkfdrData2019 documentation built on June 20, 2021, 10 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set(echo = TRUE)
Introduction
This is an analysis of the yeast data with 48 biological replicates in each of two conditions (analyzed in this publication). We chose this experiment because of the large number of biological replicates, which will allow us to (1) implement null comparisons on random subsets of samples within one condition, and (2) start with a null comparison and add in artificial differences to a subset of genes to define 'true positives'.
We will simulate a strongly informative and weakly informative independent covariate. We will also investigate the effects of an uninformative covariate by using a randomly generated covariate.
Processed count table is made available by the authors in their paper codebase GitHub repository.
In this Rmd we will carry out simulations for multiple
replicates, and plot results averaging over the replications. In contrast to
yeast-simulation.Rmd we will use an alternate simulation framework for generating
the effect sizes of the DE genes. Here we will use the fold changes of only genes
that were DE in a comparison of wild type versus knockout (with sample size 10).
Since this results in a distribution of effect sizes (log2 fold changes) under the
alternative that does not have a mode at zero, the assumptions for ashq are violated.
However, since this assumption is impossible to check, we still include this
method in the benchmark comparisons.
In this setting, we assume 500 DE genes and a bimodal alternative.
Set up workspace
# Load packages and source benchmark FDR library(SummarizedBenchmark) library(data.table) library(tidyr) library(dplyr) library(readr) library(ggplot2) library(magrittr) library(cowplot) library(purrr) library(DESeq2) library(tibble) library(ggthemes) library(R.utils) library(Hotelling) ## load helper functions for (f in list.files("../R", "\\.(r|R)$", full.names = TRUE)) { source(f) } # set up data / results directories datdir <- "yeast-data" resdir <- "../../results/RNAseq" dir.create(datdir, showWarnings = FALSE) dir.create(resdir, showWarnings = FALSE) # results files that will be generated # strong covariate resfile_d5 <- file.path(resdir, "yeastIIH-results-de5.rds") resfile_d10 <- file.path(resdir, "yeastIIH-results-de10.rds") # weak covariate resfile_d5_w <- file.path(resdir, "yeastIIHW-results-de5.rds") resfile_d10_w <- file.path(resdir, "yeastIIHW-results-de10.rds") # uninformative covariate resfile_d5_uninfCov <- file.path(resdir, "yeastIIH-results-de5-uninfCov.rds") resfile_d10_uninfCov <- file.path(resdir, "yeastIIH-results-de10-uninfCov.rds") # set up parallel backend library(parallel) nCores <- 20
Data download
First, we download the processed count data from GitHub. There is one file for the
Snf2 condition and one file for the wild type condition. The Snf2 condition is a yeast
strain that has the transcriptional regulator gene Snf2 knocked out. Each of these
files is a compressed tar.gz archive that contains a
single bam file for each replicate in each condition, which we'll place in a subdirectory
called r datdir.
download.file(url = "https://github.com/bartongroup/profDGE48/raw/master/Preprocessed_data/Snf2_countdata.tar.gz", destfile = file.path(datdir, "Snf2_countdata.tar.gz")) download.file(url = "https://github.com/bartongroup/profDGE48/raw/master/Preprocessed_data/WT_countdata.tar.gz", destfile = file.path(datdir, "WT_countdata.tar.gz")) gunzip(file.path(datdir, "Snf2_countdata.tar.gz")) gunzip(file.path(datdir, "WT_countdata.tar.gz")) untar(file.path(datdir, "Snf2_countdata.tar"), exdir = datdir) untar(file.path(datdir, "WT_countdata.tar"), exdir = datdir) file.remove(file.path(datdir, "Snf2_countdata.tar"), file.path(datdir, "WT_countdata.tar"))
Each of the data files contains two columns, one with a gene/feature name, and one with the count value.
We'll also download the list of 'bad' replicates which were especially poorly correlated to the others, as determined by the authors.
download.file(url = "https://github.com/bartongroup/profDGE48/raw/master/Bad_replicate_identification/exclude.lst", destfile = (file.path(datdir, "badreps.txt")))
Read data into R and create a count table
Here we make use of the map and map2 functions in the purrr package, to swiftly
apply the read_tsv function from readr to read in all of the 96 sample tables,
as well as add in the sample name (derived from the file name) to each subtable.
Finally, the reduce function is used to join all the replicates together.
We will remove samples that failed QC in the original study.
files <- dir(path = datdir, pattern = "*.bam.gbgout", full.names = TRUE) sample_names <- sapply(strsplit(dir(path = datdir, pattern = "*.bam.gbgout"), "_MID"), function(x) x[[1]]) badreps <- read_tsv(file.path(datdir, "badreps.txt"), col_names = FALSE)$X1 badreps <- unlist(lapply(strsplit(badreps, "_MID"), function(x) x[1])) counts <- files %>% purrr::map(read_tsv, col_names = FALSE) %>% # read in all the files individually purrr::map2(sample_names, ~ dplyr:::rename(.x, !! .y := X2, feature = X1) ) %>% # add sample names purrr::reduce(left_join, by = "feature") %>% # reduce with rbind into one dataframe dplyr::select(-badreps ) # remove badreps
Empirical analysis of 10 vs 10 (non-null)
Here we'll carry out an analysis of a subset of the full dataset, comparing the controls to the knockout samples. The fold changes observed in this comparison will be used when generating the non-null simulated data in the following sections. Here we use a subset of the samples (10 samples in each group) to mimic the sample sizes we use in simulation.
We're ready to construct a DESeq2 object. First we pull out the feature names and add them as rownames for the count table, and next we construct a column data object that houses the sample names, replicate numbers, and condition factor (WT versus Snf2 knockout).
Set up DESeq2 object
First we set up the DESeq2 object.
feats <- (counts %>% select(1))$feature counts <- as.matrix(counts %>% select(-1)) rownames(counts) <- feats coldat <- tibble(sample=colnames(counts)) %>% separate(sample, sep="_", into=c("condition", "replicate"), remove=FALSE) %>% mutate(condition = factor(condition)) # filter low count genes counts <- counts[rowMeans(counts) > 1,] dds_full <- DESeqDataSetFromMatrix(countData = counts, colData = coldat, design= ~ condition)
Run DESeq2
Next we run DESeq2, and subset results on genes with FDR < 0.05.
# results on full set set.seed(9384) subset.wt <- dds_full[,colData(dds_full)$condition == "WT"] subset.sn <- dds_full[,colData(dds_full)$condition == "Snf2"] subset <- cbind(subset.wt[, sample(1:ncol(subset.wt), 10)], subset.sn[, sample(1:ncol(subset.sn), 10)]) dds <- DESeq(subset) resultsNames(dds) # lists the coefficients res_10 <- results(dds, name="condition_WT_vs_Snf2", independentFiltering = FALSE) head(res_10) sum(res_10$padj < 0.05, na.rm = TRUE)
Check densities of the test statistic and effect sizes for significant DE genes.
data.frame(res_10[res_10$padj < 0.05 & !is.na(res_10$padj),]) %>% ggplot(aes(x = stat)) + geom_density() + xlab("Test Statistic") data.frame(res_10[res_10$padj < 0.05 & !is.na(res_10$padj),]) %>% ggplot(aes(x = log2FoldChange)) + geom_density() + xlab("Effect Size (log2 Fold Change)")
Check assumptions
There are several thousand genes detected using the full set. Next, we'll build input data frame for summarized benchmark.
geneExp <- tbl_df(data.frame(geneName=rownames(res_10), pval=res_10$pvalue, SE=res_10$lfcSE, ind_covariate = res_10$baseMean, effect_size=res_10$log2FoldChange, test_statistic=res_10$stat, pzero=rowSums(counts(subset)==0)/ncol(counts(subset)))) # filter NAs and those with less than 50% expressed geneExp <- geneExp %>% na.omit() %>% dplyr::filter(pzero < 0.5)
We'll create a plot to examine the distribution of effect sizes, since the ash method assumes that the distribution of true (unobserved) effect sizes is unimodal.
ggplot(data=geneExp, aes(effect_size)) + geom_histogram(bins=30)
We'll also explore how the standard error (used by ash) correlates with the independent covariate (used by methods that incorporate covariates), in order to get an idea of how these pieces of information relate to one another.
ggplot(data=geneExp, aes(x = ind_covariate, y = SE)) + geom_hex(bins = 100) + scale_x_continuous(trans="log10") + xlab("Covariate: Mean gene expression")
Look at covariate diagnostic plots.
strat_hist(geneExp, pvalue="pval", covariate="ind_covariate", maxy=30)
rank_scatter(geneExp, pvalue="pval", covariate="ind_covariate")
FDR benchmarking
Build common bench design. We also add in Scott's FDR Regression since our test statistics are approximately t-distributed. Note that the assumption of FDRreg-empirical that the distribution of non-null test statistics does not have a significant mass at zero is violated here, but we still include fdrreg-e in the evaluation since this assumption is impossible to check outside of a simulation setting (since it is not possible to know which tests are non-null). Thus, we'd like to evaluate how results change when the assumption is violated, since it's plausible that it is violated in the case studies.
bd <- initializeBenchDesign() bd <- addBMethod(bd, "fdrreg-t", FDRreg::FDRreg, function(x) { x$FDR }, z = test_statistic, features = model.matrix( ~ splines::bs(ind_covariate, df = 3) - 1), nulltype = 'theoretical', control = list(lambda = 0.01)) bd <- addBMethod(bd, "fdrreg-e", FDRreg::FDRreg, function(x) { x$FDR }, z = test_statistic, features = model.matrix( ~ splines::bs(ind_covariate, df = 3) - 1), nulltype = 'empirical', control = list(lambda = 0.01))
Run benchmark methods.
sb <- bd %>% buildBench(data=geneExp, parallel = FALSE) assayNames(sb) <- "qvalue" sb <- addDefaultMetrics(sb)
Plot results.
rejections_scatter(sb, supplementary=FALSE)
plotFDRMethodsOverlap(sb, alpha=0.05, nsets=ncol(sb), order.by="freq", decreasing=TRUE, supplementary=FALSE)
Empirical analysis of 5 vs 5 (non-null)
Here we'll carry out an analysis of a subset of the full dataset, comparing the controls to the knockout samples. The fold changes observed in this comparison will be used when generating the non-null simulated data in the following sections. Here we use a subset of the samples (5 samples in each group) to mimic the sample sizes we use in simulation.
Run DESeq2
First we run DESeq2, and subset results on genes with FDR < 0.05.
# results on full set set.seed(728) subset.wt <- dds_full[,colData(dds_full)$condition == "WT"] subset.sn <- dds_full[,colData(dds_full)$condition == "Snf2"] subset <- cbind(subset.wt[, sample(1:ncol(subset.wt), 5)], subset.sn[, sample(1:ncol(subset.sn), 5)]) dds <- DESeq(subset) resultsNames(dds) # lists the coefficients res_5 <- results(dds, name="condition_WT_vs_Snf2", independentFiltering = FALSE) head(res_5) sum(res_5$padj < 0.05, na.rm = TRUE)
Check densities of the test statistic and effect sizes for significant DE genes.
data.frame(res_5[res_5$padj < 0.05 & !is.na(res_5$padj),]) %>% ggplot(aes(x = stat)) + geom_density() + xlab("Test Statistic") data.frame(res_5[res_5$padj < 0.05 & !is.na(res_5$padj),]) %>% ggplot(aes(x = log2FoldChange)) + geom_density() + xlab("Effect Size (log2 Fold Change)")
Check assumptions
There are several thousand genes detected using the full set. Next, we'll build input data frame for summarized benchmark.
geneExp <- tbl_df(data.frame(geneName=rownames(res_5), pval=res_5$pvalue, SE=res_5$lfcSE, ind_covariate = res_5$baseMean, effect_size=res_5$log2FoldChange, test_statistic=res_5$stat, pzero=rowSums(counts(subset)==0)/ncol(counts(subset)))) # filter NAs and those with less than 50% expressed geneExp <- geneExp %>% na.omit() %>% dplyr::filter(pzero < 0.5)
We'll create a plot to examine the distribution of effect sizes, since the ash method assumes that the distribution of true (unobserved) effect sizes is unimodal.
ggplot(data=geneExp, aes(effect_size)) + geom_histogram(bins=30)
We'll also explore how the standard error (used by ash) correlates with the independent covariate (used by methods that incorporate covariates), in order to get an idea of how these pieces of information relate to one another.
ggplot(data=geneExp, aes(x = ind_covariate, y = SE)) + geom_hex(bins = 100) + scale_x_continuous(trans="log10") + xlab("Covariate: Mean gene expression")
Look at covariate diagnostic plots.
strat_hist(geneExp, pvalue="pval", covariate="ind_covariate", maxy=23)
rank_scatter(geneExp, pvalue="pval", covariate="ind_covariate")
FDR benchmarking
Run benchmark methods.
sb <- bd %>% buildBench(data=geneExp, parallel = FALSE) assayNames(sb) <- "qvalue" sb <- addDefaultMetrics(sb)
Plot results.
rejections_scatter(sb, supplementary=FALSE)
plotFDRMethodsOverlap(sb, alpha=0.05, nsets=ncol(sb), order.by="freq", decreasing=TRUE, supplementary=FALSE)
Simulation set up
First we'll illustrate the logistic curve used to simulate the informative covariate. In this study we'll sample from a logistic curve and use these values as probability weights when selecting DE genes. A weakly informative covariate will be generated by adding noise to the strongly informative covariate.
x <- seq(0,10, by = 0.1) infcov <- 1 / (1+exp(-x + 5)) plot(x, infcov, type = "l", ylab = "informative covariate") plot(infcov, pmin(1, abs(infcov + rnorm(length(x), 0, 0.25))), ylab = "weak covariate", xlab = "true covariate")
Next, here's a function that will generate sampling weights for test statistics of DE genes that downweights observations near zero, as well as extremely large values.
stat_wt <- function(z){ pos <- pmax(0,dgamma(z,4.5, 1-1e-4)) neg <- pmax(0,dgamma(-z,4.5, 1-1e-4)) res <- pos res[z < 0] <- neg[z < 0] return(res) } # Plot weights by test statistic value df = data.frame(z = seq(-10,10, by = 0.1)) df$w = stat_wt(df$z) ggplot(data=df, aes(x=z, y=w)) + geom_line() + ylab("Sampling Weight") + xlab("Test statistic value")
Next, we'll analyze random splits of one condition, both with and without the addition of simulated DE genes. Here we'll create a function that we can use to run one replicate given sample size and number of DE gene settings. This will be looped over many replications and results averaged over them.
#' @param X simnumber #' @param seed random seed #' @param sampleSize is the number of samples in each condition #' @param nDE is the number of DE genes #' @param bd is the bench design object #' @param BPPARAM is the BiocParallel bpparam argument to pass to DESeq. Since #' parallelization happens by forking multiple instances of this function, we #' will specify only one worker per DEseq intance. #' @param uninformativeCovariate logical indicating whether to replace the #' informative covariate with a randomly generated (uninformative) covariate. #' Default is FALSE. #' @param pvalHists logical whether to return histograms of pvalues by covariate #' instead of SB objects simulateOneSplit <- function(X, rseed, nDE, sampleSize, bd, uninformativeCovariate = FALSE, pvalHists = FALSE, strongCovariate = TRUE){ # set random seed set.seed(as.numeric(X)*as.numeric(rseed)) covar <- 1 / (1+exp(-runif(nrow(dds_full), 0, 10) + 5)) # select a random subset of 20 WT samples dds_test <- dds_full[,sample(1:ncol(dds_full), sampleSize*2)] # add a fake condition column to coldat colData(dds_test)$fake <- factor(c(rep("A", sampleSize), rep("B", sampleSize))[sample(1:(sampleSize*2), sampleSize*2)]) design(dds_test) <- ~fake pzero = rowSums(counts(dds_test)==0)/ncol(counts(dds_test)) dds_test <- dds_test[pzero < 0.5,] covar <- covar[pzero < 0.5] truth <- rep(FALSE, nrow(dds_test)) # make sure null comparison is truly null: if PC1 or PC2 sig different, # or if test of join means of PCs 1:4 is significant, # reshuffle sample labels dds_test <- estimateSizeFactors(dds_test) x <- t(counts(dds_test, normalize=TRUE)) pc <- prcomp(log(x + 0.5), scale.=TRUE) a1 <- pc$x[colData(dds_test)$fake=="A",1] b1 <- pc$x[colData(dds_test)$fake=="B",1] a2 <- pc$x[colData(dds_test)$fake=="A",2] b2 <- pc$x[colData(dds_test)$fake=="B",2] p1 <- t.test(a1, b1)$p.value p2 <- t.test(a2, b2)$p.value tries <- 0 while(p1 < 0.10 || p2 < 0.10 && tries < 10){ colData(dds_test)$fake <- sample(colData(dds_test)$fake, ncol(dds_test)) x <- t(counts(dds_test, normalize=TRUE)) pc <- prcomp(log(x + 0.5), scale.=TRUE) a1 <- pc$x[colData(dds_test)$fake=="A",1] b1 <- pc$x[colData(dds_test)$fake=="B",1] a2 <- pc$x[colData(dds_test)$fake=="A",2] b2 <- pc$x[colData(dds_test)$fake=="B",2] p1 <- t.test(a1, b1)$p.value p2 <- t.test(a2, b2)$p.value tries <- tries + 1 } if(nDE > 0){ # pick random set of nDE genes to add signal to DE <- sample(1:nrow(dds_test), nDE, prob = covar) truth[DE] <- TRUE # randomly sample a test statistic from original test statistics, with probability # weights that downweight observations near zero counts_new <- counts(dds_test) stat <- rep(0, nrow(dds_test)) # assume fixed standard error (median observed) and calculate FC from stat if(sampleSize == 5){ stat[DE] <- sample(res_5$stat, nDE, prob = stat_wt(res_5$stat)) log2FC <- stat * median(res_5$lfcSE) }else if (sampleSize == 10){ stat[DE] <- sample(res_10$stat, nDE, prob = stat_wt(res_10$stat)) log2FC <- stat * median(res_10$lfcSE) }else{ stop("Only sample sizes 5 and 10 are currently supported with pre-", "computed fold changes for sampling from.") } # randomize which condition is shifted up or down ran <- runif(nrow(dds_test)) refcond <- ifelse(ran < 0.5, "A", "B") down <- which(ran < 0.5) counts_new[down,colData(dds_test)$fake==unique(refcond[down])] <- counts(dds_test)[down, colData(dds_test)$fake==unique(refcond[down])] * 2^log2FC[down] counts_new[-down,colData(dds_test)$fake==unique(refcond[-down])] <- counts(dds_test)[-down, colData(dds_test)$fake==unique(refcond[-down])] * 2^log2FC[-down] counts_new <- apply(counts_new, 2, as.integer) counts(dds_test) <- counts_new } # replace existing size factors dds_test <- estimateSizeFactors(dds_test) dds_test <- DESeq(dds_test, parallel = FALSE) resTEST <- results(dds_test, name="fake_B_vs_A", independentFiltering = FALSE) geneExp <- tbl_df(data.frame(geneName=rownames(resTEST), pval=resTEST$pvalue, SE=resTEST$lfcSE, ind_covariate = covar, effect_size = resTEST$log2FoldChange, test_statistic = resTEST$stat, qvalue = truth)) if (uninformativeCovariate){ geneExp <- mutate(geneExp, ind_covariate = runif(length(covar))) }else if(!strongCovariate){ geneExp <- mutate(geneExp, ind_covariate = pmin(1, abs(covar + rnorm(length(covar), 0, 0.25)))) } geneExp <- geneExp %>% dplyr::filter(!is.na(pval)) if(pvalHists){ return(strat_hist(geneExp, pvalue="pval", covariate="ind_covariate", maxy =10)) } sb <- bd %>% buildBench(data=geneExp, parallel = FALSE, truthCols = "qvalue", ftCols = "ind_covariate") assayNames(sb) <- "qvalue" sb <- addDefaultMetrics(sb) rowData(sb)$log2FC <- geneExp$effect_size return(sb) }
We'll also set some parameters that will be common to all simulations. These include the number of replications, the bench design object, the set of methods to exclude in the results plots, and the alpha cutoff level to be used when plotting the aggregated Upset results.
B <- 100 excludeSet <- c("unadjusted", "bl-df02", "bl-df04", "bl-df05") ualpha <- 0.05 # only keep one condition for subsetting in simulations that follow dds_full <- dds_full[,colData(dds_full)$condition == "Snf2"]
Here's a helper function to return the number of methods with rejections at a particular alpha level (this helps us determine whether or not to plot the aggregated upset plot - if there aren't at least 2 methods it will throw an error, which is a problem for the null simulations).
# To be included in the upset agg plot, method needs to have found on average # at least one rejection per replicate. To create an upset plot, require that # at least two methods rejected at this threshold. #' @param res standardized metric data.table generated using #' standardize_results. #' @param alpha alpha cutoff #' @param filterSet which methods to exclude from consideration numberMethodsReject <- function(res, alphacutoff, filterSet){ res <- res %>% filter(is.na(param.alpha) | (param.alpha == alphacutoff)) %>% filter(!(blabel %in% filterSet)) %>% filter(alpha == alphacutoff) %>% filter(performanceMetric == "rejections") %>% select(blabel, performanceMetric, value) %>% group_by(blabel) %>% summarize(mean_value = mean(value)) %>% filter(mean_value > 1) return(nrow(res)) }
Non-null Comparisons with Strongly Informative Covariate
D5S: DE 5v5 Strong
Here we'll repeat the above, but for a DE comparison 5 versus 5 Snf2 samples, where the groups are selected randomly and 500 DE genes are added. This will be done for 100 random splits.
Generate a list of SB results objects
rseed <- 198 sampleSize <- 5 nDE <- 500 if (!file.exists(resfile_d5)){ de5 <- mclapply(X=1:B, FUN=simulateOneSplit, rseed=rseed, nDE=nDE, sampleSize=sampleSize, bd=bd, mc.cores=nCores) saveRDS(de5, file = resfile_d5) }else{ de5 <- readRDS(file = resfile_d5) } # pvalue histograms plotfile <- file.path(datdir, "de5_pvalhists.pdf") if (!file.exists(plotfile)){ hists <- mclapply(X=1:B, FUN=simulateOneSplit, rseed=rseed, nDE=nDE, sampleSize=sampleSize, bd=bd, pvalHists = TRUE, mc.cores = nCores) pdf(plotfile, width=8, height=4) for(i in 1:length(hists)){ print(hists[[i]]) } dev.off() }
Plot average results over replications
Plot results.
# Check for missing results (if any methods threw an error for relevant metrics). rowSums(sapply(de5, function(x) colSums(is.na(assays(x)$qvalue)) > 0)) res5d <- plotsim_standardize(de5, alpha = seq(0.01, 0.10, 0.01)) plotsim_average(res5d, met="rejections",filter_set = excludeSet, merge_ihw = TRUE, errorBars=TRUE) plotsim_average(res5d, met="FDR",filter_set = excludeSet, merge_ihw = TRUE, errorBars=TRUE) plotsim_average(res5d, met="TPR",filter_set = excludeSet, merge_ihw = TRUE, errorBars=TRUE) plotsim_average(res5d, met=c("FPR", "TPR"), filter_set = excludeSet, merge_ihw = TRUE) covariateLinePlot(de5, alpha=0.05, covname="log2FC", nbins=25, trans="log1p") covariateLinePlot(de5, alpha=0.05, covname="ind_covariate", nbins=25, trans="log1p") if (numberMethodsReject(res5d, alphacutoff=ualpha, filterSet=excludeSet) >= 2){ aggupset(de5, alpha=ualpha, supplementary = FALSE, return_list = FALSE) }else{ message("Not enough methods found rejections at alpha ", ualpha, "; skipping upset plot") }
D10S: DE 10v10 Strong
Here we'll repeat the above, but for a DE comparison 10 versus 10 Snf2 samples, where the groups are selected randomly and 500 DE genes are added. This will be done for 100 random splits.
Generate a list of SB results objects
rseed <- 961 sampleSize <- 10 nDE <- 500 if(!file.exists(resfile_d10)){ de10 <- mclapply(X=1:B, FUN=simulateOneSplit, rseed=rseed, nDE=nDE, sampleSize=sampleSize, bd=bd, mc.cores=nCores) saveRDS(de10, file = resfile_d10) }else{ de10 <- readRDS(file = resfile_d10) } # pvalue histograms plotfile <- file.path(datdir, "de10_pvalhists.pdf") if (!file.exists(plotfile)){ hists <- mclapply(X=1:B, FUN=simulateOneSplit, rseed=rseed, nDE=nDE, sampleSize=sampleSize, bd=bd, pvalHists = TRUE, mc.cores = nCores) pdf(plotfile, width=8, height=4) for(i in 1:length(hists)){ print(hists[[i]]) } dev.off() }
Plot average results over replications
Plot results.
# Check for missing results (if any methods threw an error for relevant metrics). rowSums(sapply(de10, function(x) colSums(is.na(assays(x)$qvalue)) > 0)) res10d <- plotsim_standardize(de10, alpha = seq(0.01, 0.10, 0.01)) plotsim_average(res10d, met="rejections",filter_set = excludeSet, merge_ihw = TRUE, errorBars=TRUE) plotsim_average(res10d, met="FDR",filter_set = excludeSet, merge_ihw = TRUE, errorBars=TRUE) plotsim_average(res10d, met="TPR",filter_set = excludeSet, merge_ihw = TRUE, errorBars=TRUE) plotsim_average(res10d, met=c("FPR", "TPR"), filter_set = excludeSet, merge_ihw = TRUE) covariateLinePlot(de10, alpha=0.05, covname="log2FC", nbins=25, trans="log1p") covariateLinePlot(de10, alpha=0.05, covname="ind_covariate", nbins=25, trans="log1p") if (numberMethodsReject(res10d, alphacutoff=ualpha, filterSet=excludeSet) >= 2){ aggupset(de10, alpha=ualpha, supplementary = FALSE, return_list = FALSE) }else{ message("Not enough methods found rejections at alpha ", ualpha, "; skipping upset plot") }
Non-null Comparisons with Weakly Informative Covariate
D5W: DE 5v5 Weak
Here we'll repeat the above, but using a weaker informative covariate (that has noise added to it).
Generate a list of SB results objects
rseed <- 198 sampleSize <- 5 nDE <- 500 if (!file.exists(resfile_d5_w)){ de5 <- mclapply(X=1:B, FUN=simulateOneSplit, rseed=rseed, nDE=nDE, sampleSize=sampleSize, bd=bd, mc.cores=nCores, strongCovariate = FALSE) saveRDS(de5, file = resfile_d5_w) }else{ de5 <- readRDS(file = resfile_d5_w) }
Plot average results over replications
Plot results.
# Check for missing results (if any methods threw an error for relevant metrics). rowSums(sapply(de5, function(x) colSums(is.na(assays(x)$qvalue)) > 0)) res5d <- plotsim_standardize(de5, alpha = seq(0.01, 0.10, 0.01)) plotsim_average(res5d, met="rejections",filter_set = excludeSet, merge_ihw = TRUE, errorBars=TRUE) plotsim_average(res5d, met="FDR",filter_set = excludeSet, merge_ihw = TRUE, errorBars=TRUE) plotsim_average(res5d, met="TPR",filter_set = excludeSet, merge_ihw = TRUE, errorBars=TRUE) plotsim_average(res5d, met=c("FPR", "TPR"), filter_set = excludeSet, merge_ihw = TRUE) covariateLinePlot(de5, alpha=0.05, covname="log2FC", nbins=25, trans="log1p") covariateLinePlot(de5, alpha=0.05, covname="ind_covariate", nbins=25, trans="log1p") if (numberMethodsReject(res5d, alphacutoff=ualpha, filterSet=excludeSet) >= 2){ aggupset(de5, alpha=ualpha, supplementary = FALSE, return_list = FALSE) }else{ message("Not enough methods found rejections at alpha ", ualpha, "; skipping upset plot") }
D10W: DE 10v10 Weak
Here we'll repeat the above, but using a weaker informative covariate (that has noise added to it).
Generate a list of SB results objects
rseed <- 961 sampleSize <- 10 nDE <- 500 if(!file.exists(resfile_d10_w)){ de10 <- mclapply(X=1:B, FUN=simulateOneSplit, rseed=rseed, nDE=nDE, sampleSize=sampleSize, bd=bd, mc.cores=nCores, strongCovariate = FALSE) saveRDS(de10, file = resfile_d10_w) }else{ de10 <- readRDS(file = resfile_d10_w) }
Plot average results over replications
Plot results.
# Check for missing results (if any methods threw an error for relevant metrics). rowSums(sapply(de10, function(x) colSums(is.na(assays(x)$qvalue)) > 0)) res10d <- plotsim_standardize(de10, alpha = seq(0.01, 0.10, 0.01)) plotsim_average(res10d, met="rejections",filter_set = excludeSet, merge_ihw = TRUE, errorBars=TRUE) plotsim_average(res10d, met="FDR",filter_set = excludeSet, merge_ihw = TRUE, errorBars=TRUE) plotsim_average(res10d, met="TPR",filter_set = excludeSet, merge_ihw = TRUE, errorBars=TRUE) plotsim_average(res10d, met=c("FPR", "TPR"), filter_set = excludeSet, merge_ihw = TRUE) covariateLinePlot(de10, alpha=0.05, covname="log2FC", nbins=25, trans="log1p") covariateLinePlot(de10, alpha=0.05, covname="ind_covariate", nbins=25, trans="log1p") if (numberMethodsReject(res10d, alphacutoff=ualpha, filterSet=excludeSet) >= 2){ aggupset(de10, alpha=ualpha, supplementary = FALSE, return_list = FALSE) }else{ message("Not enough methods found rejections at alpha ", ualpha, "; skipping upset plot") }
Non-null Comparisons with Uninformative Covariate
Here we repeat the previous sections, using an uninformative covariate.
D5U: DE 5v5 Uninformative
Here we'll repeat the above, but for a DE comparison 5 versus 5 Snf2 samples, where the groups are selected randomly and 500 DE genes are added. This will be done for 100 random splits.
Generate a list of SB results objects
rseed <- 198 sampleSize <- 5 nDE <- 500 if (!file.exists(resfile_d5_uninfCov)){ de5 <- mclapply(X=1:B, FUN=simulateOneSplit, rseed=rseed, nDE=nDE, sampleSize=sampleSize, bd=bd, uninformativeCovariate = TRUE, mc.cores=nCores) saveRDS(de5, file = resfile_d5_uninfCov) }else{ de5 <- readRDS(file = resfile_d5_uninfCov) }
Plot average results over replications
Plot results.
# Check for missing results (if any methods threw an error for relevant metrics). rowSums(sapply(de5, function(x) colSums(is.na(assays(x)$qvalue)) > 0)) res5d <- plotsim_standardize(de5, alpha = seq(0.01, 0.10, 0.01)) plotsim_average(res5d, met="rejections",filter_set = excludeSet, merge_ihw = TRUE, errorBars=TRUE) plotsim_average(res5d, met="FDR",filter_set = excludeSet, merge_ihw = TRUE, errorBars=TRUE) plotsim_average(res5d, met="TPR",filter_set = excludeSet, merge_ihw = TRUE, errorBars=TRUE) plotsim_average(res5d, met=c("FPR", "TPR"), filter_set = excludeSet, merge_ihw = TRUE) covariateLinePlot(de5, alpha=0.05, covname="log2FC", nbins=25, trans="log1p") covariateLinePlot(de5, alpha=0.05, covname="ind_covariate", nbins=25, trans="log1p") if (numberMethodsReject(res5d, alphacutoff=ualpha, filterSet=excludeSet) >= 2){ aggupset(de5, alpha=ualpha, supplementary = FALSE, return_list = FALSE) }else{ message("Not enough methods found rejections at alpha ", ualpha, "; skipping upset plot") }
D10U: DE 10v10 Uninformative
Here we'll repeat the above, but for a DE comparison 10 versus 10 Snf2 samples, where the groups are selected randomly and 500 DE genes are added. This will be done for 100 random splits.
Generate a list of SB results objects
rseed <- 961 sampleSize <- 10 nDE <- 500 if(!file.exists(resfile_d10_uninfCov)){ de10 <- mclapply(X=1:B, FUN=simulateOneSplit, rseed=rseed, nDE=nDE, sampleSize=sampleSize, bd=bd, uninformativeCovariate = TRUE, mc.cores=nCores) saveRDS(de10, file = resfile_d10_uninfCov) }else{ de10 <- readRDS(file = resfile_d10_uninfCov) }
Plot average results over replications
Plot results.
# Check for missing results (if any methods threw an error for relevant metrics). rowSums(sapply(de10, function(x) colSums(is.na(assays(x)$qvalue)) > 0)) res10d <- plotsim_standardize(de10, alpha = seq(0.01, 0.10, 0.01)) plotsim_average(res10d, met="rejections",filter_set = excludeSet, merge_ihw = TRUE, errorBars=TRUE) plotsim_average(res10d, met="FDR",filter_set = excludeSet, merge_ihw = TRUE, errorBars=TRUE) plotsim_average(res10d, met="TPR",filter_set = excludeSet, merge_ihw = TRUE, errorBars=TRUE) plotsim_average(res10d, met=c("FPR", "TPR"), filter_set = excludeSet, merge_ihw = TRUE) covariateLinePlot(de10, alpha=0.05, covname="log2FC", nbins=25, trans="log1p") covariateLinePlot(de10, alpha=0.05, covname="ind_covariate", nbins=25, trans="log1p") if (numberMethodsReject(res10d, alphacutoff=ualpha, filterSet=excludeSet) >= 2){ aggupset(de10, alpha=ualpha, supplementary = FALSE, return_list = FALSE) }else{ message("Not enough methods found rejections at alpha ", ualpha, "; skipping upset plot") }
Simulation setting comparison
Here we compare the method ranks for the different sample sizes, weighting, and covariate settings at alpha = 0.10.
plotMethodRanks(c(resfile_d5, resfile_d5_w, resfile_d10, resfile_d10_w, resfile_d5_uninfCov, resfile_d10_uninfCov), colLabels = c("5 S", "5 W", "10 S", "10 W", "5 UI", "10 UI"), alpha = 0.10, xlab = "Comparison", excludeMethods = NULL)
Session information
sessionInfo()
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.