tTestPerRow.plotVolcano: Volcano plot from a tTestPerRow result object.

In jvanheld/stats4bioinfo: Utilities for the book "Statistics for bioinformatics"

Description

Draw a volcano plot from the table produced by stats4bioinfo::tTestPerRow(). Abcissa represents the mean difference between groups. Ordinate the significance of the Welch test (sig=-log10(e-value)).

Usage

tTestPerRow.plotVolcano(ttpr.result, control.type = "fdr",
  alpha = ttpr.result$alpha, ylab = paste(sep = "", "-log10(", control.type,
  ")"), plot.ci = FALSE, ...)

Arguments

ttpr.result	Must be the result obtained with tTestPerRow()
...	Additional parameters are passed to VolcanoPlot()
control.type="fdr"	Type of control (supported: "fdr", "e.value", "p.value")
alpha=0.05	Alpha threshold for the control of false positives
plot.ci=FALSE	if TRUE, draw confidence intervals for each feature of the volcano plot.

Value

no return object

Author(s)

Jacques van Helden ([e-mail])

Examples

## Parameters to generate and analyse the simulated dataset
n.h1 <- 500 ## Number of rows under H1
n.h0 <- 500 ## Number of rows under H0
alpha <- 0.05
sample.labels <- c(rep("a", 50), rep("b", 50))

## Generate an artificial dataset with n.h1 truly positive and n.h0 truly negative rows
samples.per.groups <- unlist(table(sample.labels))
h1.data <- rnormPerGroup(n=samples.per.groups, mean=c(0, 0.5), sd=c(1,1), nrow=n.h1)
h0.data <- rnormPerGroup(n=samples.per.groups, mean=c(0, 0), sd=c(1,1), nrow=n.h0)
x <- rbind(h1.data$x, h0.data$x)
x.status <- rep(c("H1","H0"), times=c(n.h1, n.h0))
table(x.status)

## Run Student t-test on the simulated data set
x.student <- tTestPerRow(x, var.equal=TRUE, cl=sample.labels, alpha=alpha, test.group="b")

## Draw a classical Volcano plot with -log10(p-value) on the Y axis. 
tTestPerRow.plotVolcano(x.student, legend.corner="topleft", control.type = "p.value")

## Add a threshold on the effect size
tTestPerRow.plotVolcano(x.student, legend.corner="topleft", control.type = "p.value", effect.threshold = 0.6)

## Draw volcano plot with status- and density-based colors
tTestPerRow.plotVolcano(x.student, legend.corner="topleft", control.type = "p.value", effect.threshold = 0.6, density.colors=TRUE)

## Draw a Volcano plot with horizontal bars denoting the confidence intervals 
## around the difference between means
tTestPerRow.plotVolcano(x.student, legend.corner="topleft", 
   plot.ci=TRUE, control.type="p.value", alpha=)

## Draw a volcano plot where the colors represent the true status (TP, FP, TN, FN) 
## rather than the significance.
(confusion.table <- table(x.status, positive=x.student$table$p.value < alpha))
pred.status <- rep(NA, times=nrow(x.student$table))
pred.status[x.status=="H1" & x.student$table$p.value < alpha] <- "TP" ## True positives
pred.status[x.status=="H0" & x.student$table$p.value >= alpha] <- "TN" ## True negatives
pred.status[x.status=="H0" & x.student$table$p.value < alpha] <- "FP" ## False positive
pred.status[x.status=="H1" & x.student$table$p.value >= alpha] <- "FN" ## False negatives
table(pred.status)
pred.status.colors <- c("TP"="darkgreen", "FP"="red", "FN"="orange", "TN"="grey")
tTestPerRow.plotVolcano(x.student, legend.corner=NULL, 
   plot.ci=TRUE, plot.points=TRUE, tick.size=0.08, control.type="p.value", 
   col.points=pred.status.colors[pred.status])
 legend("topleft", legend=paste(table(pred.status), names(table(pred.status))),
        col=pred.status.colors[names(table(pred.status))], lwd=2)

## Draw two volcano plots separating positive and negative tests, to show that 
## the p-value cutoff is equivalent to a selection of the confidence interval 
## that do not cross the 0 value. Set the p-value threshold to 1/N, which is 
## equivalent to set the control on e-value <= 1.
par(mfrow=c(1,2))
corrected.alpha <- 1/nrow(x.student$table)
tTestPerRow.plotVolcano(x.student, col.points=NULL, col.positive="black",
   legend.corner=NULL, plot.ci=TRUE, tick.size=0, 
   control.type="p.value", alpha=corrected.alpha, main="Positive tests")
tTestPerRow.plotVolcano(x.student, col.points="grey", col.positive=NA,
   legend.corner=NULL, plot.ci=TRUE, tick.size=0, 
   control.type="p.value", alpha=corrected.alpha, main="Negative tests")
par(mfrow=c(1,1))

jvanheld/stats4bioinfo documentation built on May 20, 2019, 5:16 a.m.