varcompseq | R Documentation |
Wrapper function for gene-by-gene association testing of RNA-seq data
varcompseq( exprmat, covariates, variables2test, sample_group = NULL, weights_var2test_condi = TRUE, cov_variables2test_eff = diag(ncol(variables2test)), which_test = c("permutation", "asymptotic"), which_weights = c("loclin", "voom", "none"), n_perm = 1000, progressbar = TRUE, parallel_comp = TRUE, nb_cores = parallel::detectCores() - 1, preprocessed = FALSE, doPlot = TRUE, gene_based_weights = FALSE, bw = "nrd", kernel = c("gaussian", "epanechnikov", "rectangular", "triangular", "biweight", "tricube", "cosine", "optcosine"), exact = FALSE, transform = TRUE, padjust_methods = c("BH", "BY", "holm", "hochberg", "hommel", "bonferroni"), lowess_span = 0.5, na.rm_varcompseq = TRUE, homogen_traj = FALSE )
exprmat |
a numeric matrix of size |
covariates |
a numeric matrix of size |
variables2test |
a numeric design matrix of size |
sample_group |
a vector of length |
weights_var2test_condi |
a logical flag indicating whether heteroscedasticity
weights computation should be conditional on both the variables to be tested
|
cov_variables2test_eff |
a matrix of size |
which_test |
a character string indicating which method to use to approximate
the variance component score test, either |
which_weights |
a character string indicating which method to use to estimate
the mean-variance relationship weights. Possibilities are |
n_perm |
the number of perturbations. Default is |
progressbar |
logical indicating whether a progress bar should be displayed when computing permutations (only in interactive mode). |
parallel_comp |
a logical flag indicating whether parallel computation
should be enabled. Only Linux and MacOS are supported, this is ignored on Windows.
Default is |
nb_cores |
an integer indicating the number of cores to be used when
|
preprocessed |
a logical flag indicating whether the expression data have
already been preprocessed (e.g. log2 transformed). Default is |
doPlot |
a logical flag indicating whether the mean-variance plot should be drawn.
Default is |
gene_based_weights |
a logical flag used for |
bw |
a character string indicating the smoothing bandwidth selection method to use. See
|
kernel |
a character string indicating which kernel should be used.
Possibilities are |
exact |
a logical flag indicating whether the non-parametric weights accounting
for the mean-variance relationship should be computed exactly or extrapolated
from the interpolation of local regression of the mean against the
variance. Default is |
transform |
a logical flag used for |
padjust_methods |
multiple testing correction method used if |
lowess_span |
smoother span for the lowess function, between 0 and 1. This gives
the proportion of points in the plot which influence the smooth at each value.
Larger values give more smoothness. Only used if |
na.rm_varcompseq |
logical: should missing values in |
homogen_traj |
a logical flag indicating whether trajectories should be considered homogeneous.
Default is |
A list with the following elements:
which_test
: a character string carrying forward the value of the 'which_test
' argument
indicating which test was perform (either "asymptotic" or "permutation").
preprocessed
: a logical flag carrying forward the value of the 'preprocessed
' argument
indicating whether the expression data were already preprocessed, or were provided as raw counts and
transformed into log-counts per million.
n_perm
: an integer carrying forward the value of the 'n_perm
' argument indicating
the number of perturbations performed (NA
if asymptotic test was performed).
genesets
: carrying forward the value of the 'genesets
' argument defining the gene sets
of interest (NULL
for gene-wise testing).
pval
: computed p-values. A data.frame
with one raw for each each gene set, or
for each gene if genesets
argument is NULL
, and with 2 columns: the first one 'rawPval
'
contains the raw p-values, the second one contains the FDR adjusted p-values (according to
the 'padjust_methods
' argument) and is named 'adjPval
'.
Agniel D & Hejblum BP (2017). Variance component score test for time-course gene set analysis of longitudinal RNA-seq data, Biostatistics, 18(4):589-604. doi: 10.1093/biostatistics/kxx005. arXiv:1605.02351.
sp_weights
vc_test_perm
vc_test_asym
p.adjust
#rm(list=ls()) nsims <- 2 #100 res <- numeric(nsims) for(i in 1:nsims){ n <- 1000 nr=5 ni=50 r <- nr*ni t <- matrix(rep(1:nr), ni, ncol=1, nrow=r) sigma <- 0.5 b0 <- 1 #under the null: b1 <- 0 y.tilde <- b0 + b1*t + rnorm(r, sd = sigma) y <- t(matrix(rnorm(n*r, sd = sqrt(sigma*abs(y.tilde))), ncol=n, nrow=r) + matrix(rep(y.tilde, n), ncol=n, nrow=r)) x <- matrix(1, ncol=1, nrow=r) #run test res_genes <- varcompseq(exprmat=y, covariates=x, variables2test=t, sample_group=rep(1:ni, each=nr), which_test="asymptotic", which_weights="none", preprocessed=TRUE) mean(res_genes$pvals[, "rawPval"]>0.05) quantile(res_genes$pvals[, "rawPval"]) res[i] <- mean(res_genes$pvals[, "rawPval"]<0.05) cat(i,"\n") } mean(res) ## Not run: b0 <- 1 b1 <- 0 y.tilde <- b0 + b1*t + rnorm(r, sd = sigma) y <- t(matrix(rnorm(n*r, sd = sqrt(sigma*abs(y.tilde))), ncol=n, nrow=r) + matrix(rep(y.tilde, n), ncol=n, nrow=r)) res_genes <- varcompseq(exprmat=y, covariates=x, variables2test=t, sample_group=rep(1:ni, each=nr), which_weights="none", preprocessed=TRUE) summary(res_genes$pvals) ## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.