Description Usage Arguments Value References See Also Examples
Wrapper function for performing gene set analysis of (potentially longitudinal) RNAseq data
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29  tcgsa_seq(
y,
x,
phi,
weights_phi_condi = TRUE,
genesets,
indiv = NULL,
Sigma_xi = diag(ncol(phi)),
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 = TRUE,
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,
R = NULL,
homogen_traj = FALSE,
na.rm_tcgsaseq = TRUE,
verbose = TRUE
)

y 
a numeric matrix of size 
x 
a numeric matrix of size 
phi 
a numeric design matrix of size 
weights_phi_condi 
a logical flag indicating whether heteroscedasticity
weights computation should be conditional on both the variable(s) to be tested

genesets 
either a vector of index or subscripts that defines which rows of 
indiv 
a vector of length 
Sigma_xi 
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 meanvariance 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 meanvariance 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 nonparametric weights accounting
for the meanvariance 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 
R 
library.size (optional, important to provide if 
homogen_traj 
a logical flag indicating whether trajectories should be considered homogeneous.
Default is 
na.rm_tcgsaseq 
logical: should missing values in 
verbose 
logical: should informative messages be printed during the
computation? 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 logcounts 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 genewise testing).
pval
: computed pvalues. 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 pvalues, the second one contains the FDR adjusted pvalues (according to
the 'padjust_methods
' argument) and is named 'adjPval
'.
Agniel D & Hejblum BP (2017). Variance component score test for timecourse gene set analysis of longitudinal RNAseq data, Biostatistics, 18(4):589604. 10.1093/biostatistics/kxx005. arXiv:1605.02351.
Law, C. W., Chen, Y., Shi, W., & Smyth, G. K. (2014). voom: Precision weights unlock linear model analysis tools for RNAseq read counts. Genome Biology, 15(2), R29.
sp_weights
vc_test_perm
vc_test_asym
p.adjust
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44  #rm(list=ls())
nsims < 2 #100
res_quant < list()
for(i in 1:2){
n < 2000#0
nr < 3
r < nr*20#4*nr#100*nr
t < matrix(rep(1:nr), r/nr, ncol=1, nrow=r)
sigma < 0.4
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 < tcgsa_seq(y, x, phi=t, genesets=lapply(0:9, function(x){x*10+(1:10)}),
Sigma_xi=matrix(1), indiv=rep(1:(r/nr), each=nr), which_test="asymptotic",
which_weights="none", preprocessed=TRUE)
res_genes < tcgsa_seq(y, x, phi=cbind(t),#, rnorm(r)), #t^2
genesets=NULL,
Sigma_xi=diag(1), indiv=rep(1:(r/nr), each=nr), which_test="asymptotic",
which_weights="none", preprocessed=TRUE)
length(res_genes$pvals[, "rawPval"])
quantile(res_genes$pvals[, "rawPval"])
res_quant[[i]] < res_genes$pvals[, "rawPval"]
}
#round(rowMeans(sapply(res_quant, quantile)), 3)
#plot(density(unlist(res_quant)))
#mean(unlist(res_quant)<0.05)
## Not run:
res_genes < tcgsa_seq(y, x, phi=t, genesets=NULL,
Sigma_xi=matrix(1), indiv=rep(1:(r/nr), each=nr), which_test="permutation",
which_weights="none", preprocessed=TRUE, n_perm=1000, parallel_comp = FALSE)
mean(res_genes$pvals$rawPval < 0.05)
summary(res_genes$pvals$FDR)
## End(Not run)

Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.