rotateStat {randRotation} | R Documentation |
This function generates rotations of data and calculates the provided statistic
on each rotation and the non-rotated (original) data.
This is the central function of the package.
rotateStat( initialised.obj, R = 10, statistic, ..., parallel = FALSE, BPPARAM = BiocParallel::bpparam() )
initialised.obj |
An initialised random rotation object as returned by |
R |
The number of resamples/rotations. Single |
statistic |
A function which takes a data matrix (same dimensions as |
... |
Further named arguments for |
parallel |
|
BPPARAM |
An optional |
The function takes an initialised randrot object
(initRandrot
) and a function that
calculates a statistic on the data. The statistic function thereby takes the
a matrix Y
as first argument. Any further arguments are passed to it
by ...
.
Together with pFdr
, this function implements
the workflow
described in (Hettegger et al. 2021).
Be aware that only data is rotated (see also
randrot
), so any additional information
including weights
, X
etc. need to be provided to
statistic
. See also package vignette and Examples
.
Parallel processing is implemented with the BiocParallel package of Bioconductor.
The default argument BiocParallel::bpparam()
for BPPARAM
returns the registered default backend.
See package documentation for further information and usage options.
If parallel = TRUE
the function calls in statistic
need to be called explicitly
with package name and "::". So e.g. calling lmFit
from the
limma
package is done with limma::lmFit(...)
, see also the
examples in the package vignette.
An object of class rotateStat
.
Peter Hettegger
Hettegger P, Vierlinger K, Weinhaeusel A (2021). “Random rotation for identifying differentially expressed genes with linear models following batch effect correction.” Bioinformatics. ISSN 1367-4803, doi: 10.1093/bioinformatics/btab063.
#set.seed(0) # Dataframe of phenotype data (sample information) # We simulate 2 sample classes processed in 3 batches pdata <- data.frame(batch = rep(1:3, c(10,10,10)), phenotype = rep(c("Control", "Cancer"), c(5,5))) features <- 100 # Matrix with random gene expression data edata <- matrix(rnorm(features * nrow(pdata)), features) rownames(edata) <- paste("feature", 1:nrow(edata)) mod1 <- model.matrix(~phenotype, pdata) # Initialisation of the random rotation class init1 <- initBatchRandrot(Y = edata, X = mod1, coef.h = 2, batch = pdata$batch) init1 # Definition of the batch effect correction procedure with subsequent calculation # of two-sided test statistics statistic <- function(., batch, mod, coef){ # The "capture.output" and "suppressMessages" simply suppress any output capture.output(suppressMessages( Y.tmp <- sva::ComBat(., batch = batch, mod) )) fit1 <- lm.fit(mod, t(Y.tmp)) abs(coef(fit1)[coef,]) } # We calculate test statistics for the second coefficient res1 <- rotateStat(initialised.obj = init1, R = 10, statistic = statistic, batch = pdata$batch, mod = mod1, coef = 2) hist(pFdr(res1))