probNonEquiv: `probNonEquiv` performs a Bayesian hypothesis test for equivalence between group means. It returns the posterior probability that |mu1-mu2|>logfc. `pvalTreat` is a wrapper to `treat` in package `limma`, which returns P-values for the same hypothesis test.

Description

probNonEquiv computes v_i=P(|theta_i| > logfc | data), where theta_i is the difference between group means for gene i. This posterior probability is based on the NNGCV model from package EBarrays, which has a formulation similar to limma in an empirical Bayes framework. Notice that the null hypothesis here is that |theta_i|

Subsequent differential expression calls are based on selecting large v_i. For instance, selecting v_i >= 0.95 guarantees that the posterior expected false discovery proportion (a Bayesian FDR analog) is below 0.05.

Usage

probNonEquiv(x, groups, logfc = log(2), minCount, method = "plugin", mc.cores=1)
pvalTreat(x, groups, logfc = log(2), minCount, p.adjust.method='none', mc.cores = 1)

Arguments

ExpressionSet containing expression levels, or list of ExpressionSets

groups

Variable in fData(x) indicating the two groups to compare (the case with more than 2 groups is not implemented).

logfc

Biologically relevant threshold for the log fold change, i.e. difference between groups means in log-scale

minCount

If specified, probabilities are only computed for rows with fData(x)$readCount >= minCount

method

Set to 'exact' for exact posterior probabilities (slower), 'plugin' for plug-in approximation (much faster). Typically both give very similar results.

mc.cores

Number of parallel processors to use. Ignored unless x is a list.

p.adjust.method

P-value adjustment method, passed on to p.adjust

Value

If x is a single ExpressionSet, probNonEquiv returns a vector with posterior probabilities (NA for rows with less than minCount reads). pvalTreat returns TREAT P-values instead.If x is a list of ExpressionSet, the function is applied to each element separately and results are returned as columns in the output matrix.

References

Rossell D, Stephan-Otto Attolini C, Kroiss M, Stocker A. Quantifying Alternative Splicing from Paired-End RNA-sequencing data. Annals of Applied Statistics, 8(1):309-330

McCarthy DJ, Smyth GK. Testing significance relative to a fold-change threshold is a TREAT. Bioinformatics, 25(6):765-771

Examples

Run this code

  #Simulate toy data
  p <- 50; n <- 10
  x <- matrix(rnorm(p*2*n),nrow=p)
  x[(p-10):p,1:n] <- x[(p-10):p,1:n] + 1.5
  x <- new("ExpressionSet",exprs=x)
  x$group <- rep(c('group1','group2'),each=n)

  #Posterior probabilities
  pp <- probNonEquiv(x, groups='group', logfc=0.5)
  d <- rowMeans(exprs(x[,1:n])) - rowMeans(exprs(x[,-1:-n]))
  plot(d,pp,xlab='Observed log-FC')
  abline(v=c(-.5,.5))

  #Check false positives
  truth <- rep(c(FALSE,TRUE),c(p-11,11))
  getRoc(truth, pp>.9)
  getRoc(truth, pp>.5)