| roast.DGEList {edgeR} | R Documentation |
Rotation gene set testing for Negative Binomial generalized linear models.
## S3 method for class 'DGEList'
roast(y, index=NULL, design=NULL, contrast=ncol(design), set.statistic="mean",
gene.weights=NULL, array.weights=NULL, weights=NULL, block=NULL, correlation,
var.prior=NULL, df.prior=NULL, trend.var=FALSE, nrot=999)
## S3 method for class 'DGEList'
mroast(y, index=NULL, design=NULL, contrast=ncol(design), set.statistic="mean",
gene.weights=NULL, array.weights=NULL, weights=NULL, block=NULL, correlation,
var.prior=NULL, df.prior=NULL, trend.var=FALSE, nrot=999, adjust.method="BH", midp=TRUE, sort="directional")
y |
|
index |
index vector specifying which rows (genes) of |
design |
design matrix |
contrast |
contrast for which the test is required. Can be an integer specifying a column of |
set.statistic |
summary set statistic. Possibilities are |
gene.weights |
optional numeric vector of weights for genes in the set. Can be positive or negative. For |
array.weights |
optional numeric vector of array weights. |
weights |
optional matrix of observation weights. If supplied, should be of same dimensions as |
block |
optional vector of blocks. |
correlation |
correlation between blocks. |
var.prior |
prior value for residual variances. If not provided, this is estimated from all the data using |
df.prior |
prior degrees of freedom for residual variances. If not provided, this is estimated using |
trend.var |
logical, should a trend be estimated for |
nrot |
number of rotations used to estimate the p-values. |
adjust.method |
method used to adjust the p-values for multiple testing. See |
midp |
logical, should mid-p-values be used in instead of ordinary p-values when adjusting for multiple testing? |
sort |
character, whether to sort output table by directional p-values ( |
This function implements a method for the ROAST gene set test from Wu et al (2010) for the digital gene expression data, eg. RNA-Seq data.
Basically, the Negative Binomial generalized linear models are fitted for count data. The fitted values are converted into z-scores, and then it calls the roast function in limma package to conduct the gene set test.
It tests whether any of the genes in the set are differentially expressed.
This allows users to focus on differential expression for any coefficient or contrast in a generalized linear model.
If contrast is not specified, the last coefficient in the model will be tested.
The arguments array.weights, block and correlation have the same meaning as they for for the lmFit function.
The arguments df.prior and var.prior have the same meaning as in the output of the eBayes function.
If these arguments are not supplied, they are estimated exactly as is done by eBayes.
The argument gene.weights allows directions or weights to be set for individual genes in the set.
The gene set statistics "mean", "floormean", "mean50" and msq are defined by Wu et al (2010).
The different gene set statistics have different sensitivities to small number of genes.
If set.statistic="mean" then the set will be statistically significantly only when the majority of the genes are differentially expressed.
"floormean" and "mean50" will detect as few as 25% differentially expressed.
"msq" is sensitive to even smaller proportions of differentially expressed genes, if the effects are reasonably large.
The output gives p-values three possible alternative hypotheses,
"Up" to test whether the genes in the set tend to be up-regulated, with positive t-statistics,
"Down" to test whether the genes in the set tend to be down-regulated, with negative t-statistics,
and "Mixed" to test whether the genes in the set tend to be differentially expressed, without regard for direction.
roast estimates p-values by simulation, specifically by random rotations of the orthogonalized residuals (Langsrud, 2005), so p-values will vary slightly from run to run.
To get more precise p-values, increase the number of rotations nrot.
The p-value is computed as (b+1)/(nrot+1) where b is the number of rotations giving a more extreme statistic than that observed (Phipson and Smyth, 2010).
This means that the smallest possible p-value is 1/(nrot+1).
mroast does roast tests for multiple sets, including adjustment for multiple testing.
By default, mroast reports ordinary p-values but uses mid-p-values (Routledge, 1994) at the multiple testing stage.
Mid-p-values are probably a good choice when using false discovery rates (adjust.method="BH") but not when controlling the family-wise type I error rate (adjust.method="holm").
roast performs a self-contained test in the sense defined by Goeman and Buhlmann (2007).
For a competitive gene set test, see camera.DGEList.
roast produces an object of class Roast. See roast for details.
mroast produces a data.frame. See mroast for details.
Yunshun Chen and Gordon Smyth
Goeman, JJ, and Buhlmann, P (2007). Analyzing gene expression data in terms of gene sets: methodological issues. Bioinformatics 23, 980-987.
Langsrud, O (2005). Rotation tests. Statistics and Computing 15, 53-60.
Phipson B, and Smyth GK (2010). Permutation P-values should never be zero: calculating exact P-values when permutations are randomly drawn. Statistical Applications in Genetics and Molecular Biology, Volume 9, Article 39.
Routledge, RD (1994). Practicing safe statistics with the mid-p. Canadian Journal of Statistics 22, 103-110.
Wu, D, Lim, E, Francois Vaillant, F, Asselin-Labat, M-L, Visvader, JE, and Smyth, GK (2010). ROAST: rotation gene set tests for complex microarray experiments. Bioinformatics 26, 2176-2182. http://bioinformatics.oxfordjournals.org/content/26/17/2176
mu <- matrix(10, 100, 4) group <- factor(c(0,0,1,1)) design <- model.matrix(~group) # First set of 10 genes that are genuinely differentially expressed iset1 <- 1:10 mu[iset1,3:4] <- mu[iset1,3:4]+10 # Second set of 10 genes are not DE iset2 <- 11:20 # Generate counts and create a DGEList object y <- matrix(rnbinom(100*4, mu=mu, size=10),100,4) y <- DGEList(counts=y, group=group) # Estimate dispersions y <- estimateDisp(y, design) roast(y, iset1, design, contrast=2) mroast(y, iset1, design, contrast=2) mroast(y, list(set1=iset1, set2=iset2), design, contrast=2)