probNonEquiv {casper} | R Documentation

## `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|<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.

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

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

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

### Author(s)

Victor Pena, David Rossell

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

### See Also

`treat` in package `limma`, `p.adjust`

### Examples

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

[Package casper version 2.24.2 Index]