probNonEquiv: 'probNonEquiv' performs a Bayesian hypothesis test for...
In casper: Characterization of Alternative Splicing based on Paired-End Reads

Description Usage Arguments Value Author(s) References See Also Examples

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|<logfc, e.g. isoforms with small fold changes are regarded as uninteresting.

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.

1
2
3

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)

`x`	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`

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.

Victor Pena, David Rossell

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

treat in package limma, p.adjust

  #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)