inst/doc/independent_filtering_plots.R

In genefilter: genefilter: methods for filtering genes from high-throughput experiments

## ----knitr, echo=FALSE, results="hide"----------------------------------------
library("knitr")
opts_chunk$set(tidy=FALSE,dev="png",fig.show="hide",
               fig.width=4,fig.height=4.5,dpi=240,
               message=FALSE,error=FALSE,warning=FALSE)

## ----style, eval=TRUE, echo=FALSE, results="asis"--------------------------
BiocStyle:::latex()

## ----setup, echo=FALSE--------------------------------------------------------
options( width = 80 )

## ----libraries----------------------------------------------------------------
library("genefilter")
library("ALL")
data("ALL")

## ----sample_data, cache=TRUE--------------------------------------------------
bcell <- grep("^B", as.character(ALL$BT))
moltyp <- which(as.character(ALL$mol.biol) %in% 
                c("NEG", "BCR/ABL"))
ALL_bcrneg <- ALL[, intersect(bcell, moltyp)]
ALL_bcrneg$mol.biol <- factor(ALL_bcrneg$mol.biol)
n1 <- n2 <- 3
set.seed(1969)
use <- unlist(tapply(1:ncol(ALL_bcrneg), 
                     ALL_bcrneg$mol.biol, sample, n1))
subsample <- ALL_bcrneg[,use]

## ----stats, cache=TRUE--------------------------------------------------------
S <- rowSds( exprs( subsample ) )
temp <- rowttests( subsample, subsample$mol.biol )
d <- temp$dm
p <- temp$p.value
t <- temp$statistic

## ----filter_volcano, include=FALSE--------------------------------------------
S_cutoff <- quantile(S, .50)
filter_volcano(d, p, S, n1, n2, alpha=.01, S_cutoff)

## ----kappa, include=FALSE-----------------------------------------------------
t <- seq(0, 5, length=100)
plot(t, kappa_t(t, n1, n2) * S_cutoff, 
     xlab="|T|", ylab="Fold change bound", type="l")

## ----table--------------------------------------------------------------------
table(ALL_bcrneg$mol.biol)

## ----filtered_p---------------------------------------------------------------
S2 <- rowVars(exprs(ALL_bcrneg))
p2 <- rowttests(ALL_bcrneg, "mol.biol")$p.value
theta <- seq(0, .5, .1)
p_bh <- filtered_p(S2, p2, theta, method="BH")

## ----p_bh---------------------------------------------------------------------
head(p_bh)

## ----rejection_plot-----------------------------------------------------------
rejection_plot(p_bh, at="sample",
               xlim=c(0,.3), ylim=c(0,1000),
               main="Benjamini & Hochberg adjustment")

## ----filtered_R---------------------------------------------------------------
theta <- seq(0, .80, .01)
R_BH <- filtered_R(alpha=.10, S2, p2, theta, method="BH")

## ----R_BH---------------------------------------------------------------------
head(R_BH)

## ----filtered_R_plot----------------------------------------------------------
plot(theta, R_BH, type="l",
     xlab=expression(theta), ylab="Rejections",
     main="BH cutoff = .10"
     )

## ----sessionInfo, results='asis', echo=FALSE----------------------------------
si <- as.character( toLatex( sessionInfo() ) )
cat( si[ -grep( "Locale", si ) ], sep = "\n" )