## ---- eval=FALSE-------------------------------------------------------------- # install.packages("BiocManager") ## ---- eval=FALSE-------------------------------------------------------------- # library("BiocManager") # BiocManager::install("swfdr") ## ----------------------------------------------------------------------------- library(swfdr) ## ---- eval=FALSE-------------------------------------------------------------- # BiocManager::install("qvalue") # library(qvalue) ## ---- echo=FALSE, eval=TRUE--------------------------------------------------- library(qvalue) ## ----------------------------------------------------------------------------- N <- 1000 ## ----D0_raw, cache=TRUE------------------------------------------------------- set.seed(12345) dataset0 <- matrix(rnorm(6*N), ncol=6) head(dataset0, 3) ## ----D0, cache=TRUE----------------------------------------------------------- # compute t-tests on each row of a matrix dataset_p <- function(D, group1=c(1,2,3), group2=c(4,5,6)) { t.test.groups <- function(x) { t.test(x[group1], x[group2])$p.value } round(apply(D, 1, t.test.groups), 5) # rounding for legibility } D0 <- data.frame(p=dataset_p(dataset0), truth=0) summary(D0$p) ## ----D1, cache=TRUE----------------------------------------------------------- dataset1 <- matrix(c(rnorm(3*N, 0, 1), rnorm(3*N, 2, 1)), ncol=6) head(dataset1, 3) D1 <- data.frame(p=dataset_p(dataset1), truth=1) summary(D1$p) ## ----------------------------------------------------------------------------- D2 <- rbind(D0, D1) dim(D2) table(D2$truth) ## ----Dhist, fig.width=7.5, fig.height=2.2------------------------------------- par(mfrow=c(1,3), las=1, mar=c(3,4,2,1)) p_hist <- function(p, breaks=20, freq=F, ylim=c(0, 7.2), ...) { hist(p, breaks=breaks, freq=freq, ylim=ylim, xlab="", ...) } p_hist(D0$p, col="#666666", main="p-values in D0") p_hist(D1$p, col="#cccccc", main="p-values in D1") p_hist(D2$p, col="#999999", main="p-values in D2") ## ----Dqs, cache=TRUE---------------------------------------------------------- library(qvalue) D0q <- qvalue(D0$p) D1q <- qvalue(D1$p) D2q <- qvalue(D2$p) c(D0=D0q$pi0, D1=D1q$pi0, D2=D2q$pi0) ## ----------------------------------------------------------------------------- summary(D0q) ## ----------------------------------------------------------------------------- # extract number of hits with values below a threshold hit_counts <- function(x, threshold=0.05) { counts <- c(sum(x$pvalues