Methods in Biostatistics 2 (140.652)

Multiple comparisons

For this part, you need to install the ‘qvalue’ package from Bioconductor. Uncomment the first two lines for the installation. Note this can take a wee bit.

# source("http://bioconductor.org/biocLite.R")
# biocLite("qvalue")
library(qvalue)

Null results

set.seed(1)
p <- runif(1e5)
hist(p,breaks=101,col="lightgrey",prob=T)
abline(h=1,lwd=2,col="red")

Example 1: enrichment for low p-values

p <- c(runif(1e5),rbeta(1e3,1,1.5)/10)
hist(p,breaks=101,col="lightgrey",prob=T)
abline(h=1,lwd=2,col="red")

# Bonferroni
0.05/length(p)

## [1] 4.950495e-07

min(p)

## [1] 2.781628e-06

# q-values
q <- qvalue(p)
wh <- order(q$pvalues)
cbind(q$pvalues,q$qvalues)[wh,][1:20,]

##               [,1]      [,2]
##  [1,] 2.781628e-06 0.2751326
##  [2,] 8.556060e-06 0.3871881
##  [3,] 1.303945e-05 0.3871881
##  [4,] 1.565809e-05 0.3871881
##  [5,] 1.964995e-05 0.3887179
##  [6,] 8.649053e-05 0.6954618
##  [7,] 8.902334e-05 0.6954618
##  [8,] 8.986238e-05 0.6954618
##  [9,] 9.474996e-05 0.6954618
## [10,] 9.888271e-05 0.6954618
## [11,] 1.093899e-04 0.6954618
## [12,] 1.193085e-04 0.6954618
## [13,] 1.254713e-04 0.6954618
## [14,] 1.312552e-04 0.6954618
## [15,] 1.363505e-04 0.6954618
## [16,] 1.445739e-04 0.6954618
## [17,] 1.500996e-04 0.6954618
## [18,] 1.502968e-04 0.6954618
## [19,] 1.536754e-04 0.6954618
## [20,] 1.563255e-04 0.6954618

Example 2: enrichment for low p-values

p <- c(runif(1e5),runif(10,0,1e-6))
hist(p,breaks=101,col="lightgrey",prob=T)
abline(h=1,lwd=2,col="red")

0.05/length(p)

## [1] 4.9995e-07

min(p)

## [1] 1.929471e-07

q <- qvalue(p)
wh <- order(q$pvalues)
cbind(q$pvalues,q$qvalues)[wh,][1:20,]

##               [,1]        [,2]
##  [1,] 1.929471e-07 0.009320656
##  [2,] 2.463036e-07 0.009320656
##  [3,] 4.514845e-07 0.009320656
##  [4,] 6.301916e-07 0.009320656
##  [5,] 7.234962e-07 0.009320656
##  [6,] 7.257988e-07 0.009320656
##  [7,] 7.650754e-07 0.009320656
##  [8,] 7.839321e-07 0.009320656
##  [9,] 9.010656e-07 0.009320656
## [10,] 9.319724e-07 0.009320656
## [11,] 1.172884e-05 0.106636514
## [12,] 2.192124e-05 0.182695250
## [13,] 3.785221e-05 0.291199956
## [14,] 4.770258e-05 0.340766753
## [15,] 6.872625e-05 0.458220824
## [16,] 8.272938e-05 0.515977135
## [17,] 9.252643e-05 0.515977135
## [18,] 9.286660e-05 0.515977135
## [19,] 9.864545e-05 0.519238520
## [20,] 1.115340e-04 0.557725904

Also, check out the function p.adjust() !

End of code