Installation
suppressPackageStartupMessages(library(reconsi))
cat("reconsi package version", as.character(packageVersion("reconsi")), "\n")## reconsi package version 1.0.0General use
We illustrate the general use of the package on a synthetic dataset. The default Wilcoxon rank-sum test is used.
#Create some synthetic data with 90% true null hypothesis
p = 200; n = 50
x = rep(c(0,1), each = n/2)
mat = cbind(
matrix(rnorm(n*p/10, mean = 5+x),n,p/10), #DA
matrix(rnorm(n*p*9/10, mean = 5),n,p*9/10) #Non DA
)
#Provide just the matrix and grouping factor, and test using the collapsed null
fdrRes = reconsi(mat, x)
#The estimated tail-area false discovery rates.
estFdr = fdrRes$FdrThe method provides an estimate of the proportion of true null hypothesis, which is close to the true 90%.
## [1] 0.8513562The result of the procedure can be represented graphically as follows:
## Warning: Removed 121 row(s) containing missing values (geom_path).
The approximate correlation matrix of test statistic can be visualized as follows:
It is also possible to provide a custom test function, which must accept at least a ‘y’ response variable and a ‘x’ grouping factor. Additionally, quantile, distribution and density functions should be supplied for transformation through quantiles to z-values.
#With a custom function, here linera regression
fdrResLm = reconsi(mat, x, B = 5e1,
test = function(x, y){
fit = lm(y~x)
c(summary(fit)$coef["x","t value"], fit$df.residual)},
quantileFun = function(p){qt(p = p[1], df = p[2])},
distFun = function(q){pt(q = q[1], df = q[2])},
densFun = function(x){dt(x = x[1], df = x[2])})This framework also accepts more than 2 groups, and additional covariates through the “argList” argument.
#3 groups
p = 100; n = 60
x = rep(c(0,1,2), each = n/3)
mu0 = 5
mat = cbind(
matrix(rnorm(n*p/10, mean = mu0+x),n,p/10), #DA
matrix(rnorm(n*p*9/10, mean = mu0),n,p*9/10) #Non DA
)
#Provide an additional covariate through the 'argList' argument
z = rpois(n , lambda = 2)
fdrResLmZ = reconsi(mat, x, B = 5e1,
test = function(x, y, z){
fit = lm(y~x+z)
c(summary(fit)$coef["x","t value"], fit$df.residual)},
quantileFun = function(p){qt(p = p[1], df = p[2])},
distFun = function(q){pt(q = q[1], df = q[2])},
argList = list(z = z))If the null distribution of the test statistic is not known, it is also possbile to execute the procedure on the scale of the original test statistics, rather than z-values by setting zValues = FALSE. This may be numerically less stable.
fdrResKruskal = reconsi(mat, x, B = 5e1,
test = function(x, y){kruskal.test(y~x)$statistic}, zValues = FALSE)Alternatively, the same resampling instances may be used to determine the marginal null distributions as to estimate the collapsed null distribution, by setting the “permZvals” flag to true.
fdrResKruskalPerm = reconsi(mat, x, B = 5e1,
test = function(x, y){
kruskal.test(y~x)$statistic}, permZvals = TRUE)When no grouping variable is available, one can perform a bootstrap as resampling procedure. This is achieved by simply not supplying a grouping variable “x”. Here we perform a one sample Wilcoxon rank sum test for equality of the means to 0.
Case study
We illustrate the package using an application from microbiology. The species composition of a community of microorganisms can be determined through sequencing. However, this only yields compositional information, and knowledge of the population size can be acquired by cell counting through flow cytometry. Next, the obtained species compositions can multiplied by the total population size to yield approximate absolute cell counts per species. Evidently, this introduces strong correlation between the tests due to the common factor. In other words: random noise in the estimation of the total cell counts will affect all hypotheses. Therefore, we employ resampling to estimate the collapsed null distribution, that will account for this dependence.
The dataset used is taken from Vandeputte et al., 2017 (see manuscript), a study on gut microbiome in healthy and Crohn’s disease patients. The test looks for differences in absolute abundance between healthy and diseased patients. It relies on the phyloseq package, which is the preferred way to interact with our machinery for microbiome data.
#The grouping and flow cytometry variables are present in the phyloseq object, they only need to be called by their name.
testVanDePutte = testDAA(Vandeputte, groupName = "Health.status", FCname = "absCountFrozen", B = 1e2L)The estimated tail-area false discovery rates can then simply be extracted as
## 0% 25% 50% 75% 100%
## 1.549927e-17 1.884717e-02 3.217944e-01 7.424020e-01 8.769332e-01