test_modules: Test modules for differential co-expression

In arbet003/discoMod: Identifying differentially co-expressed modules between two groups

Description

For a given module, the null hypothesis is that the network structure is the same for both phenotype groups, while the alternative hypothesis is that the network structure differs between the two groups. The following test statistics are calculated (all defined in a forthcoming manuscript), using permutations to calculate p-values. PND: p-norm difference test (exponents 4,6,8, and 20), GHD: Generalised Hamming Distance, DispIndex: Dispersion Index, MAD: mean absolute deviation, pairedT: paired t-test, and wilcoxSRT: Wilcoxon signed rank test. PND6 is the recommended default test based on simulation results from a forthcoming manuscript.

Usage

test_modules(
  group1_modules,
  group2_modules,
  cortype = "spearman",
  othertype = NULL,
  adjpower = 1,
  numperms = 500,
  perm_stats = FALSE,
  parallel = FALSE,
  cores = 2
)

Arguments

group1_modules	list of modules in phenotype group 1. Each module (element of the list) is a matrix with samples in rows, and genes in columns. You could use find_modules() to create this list, otherwise you can explore other clustering methods (e.g. hierarchical clustering) and create your own user-defined list
group2_modules	list of modules in phenotype group 2.
cortype	"spearman" (default), "pearson", or "kendall"; type of correlation used to measure the network structure within a module.
othertype	If set to NULL (default), then correlation is used to represent the network. "Adjacency" and "TOM" (topological overlap measure) are other options (only "signed" versions are considered) see WGCNA::adjacency() and WGCNA::TOMsimilarity() for details. Option "TOM" is more computationally intensive.
adjpower	if othertype="Adjacency" (or "TOM"), then adjpower represents the exponent used on all correlations (the higher the value, the more the correlations are shrunk towards zero).
numperms	number of permutations used to calculate p-values
perm_stats	TRUE or FALSE, whether or not to store all permutation test statistics (can take up a lot of memory if there are many modules and numperms is large)
parallel	TRUE or FALSE, whether or not to create a doSNOW cluster used for parallel computing (the modules are evaluated in parallel, not the permutations)
cores	number of cores to use if parallel=T. Use parallel::detectCores() to check how many cores you have available.

Value

• test_stats : observed test statistics for each module (rows correspond to the order of modules in group1_modules and group2_modules).

• pvalues : p-values calculated using permutations.

• qvalues : false discovery rate multiple testing adjusted p-values (i.e. stats::p.adjust() with method="fdr")

• perm_results : if perm_stats=TRUE then perm_results will be a list of the permutation distribution of each test statistic for all modules.

Examples

######## load example dataset
library(multtest)
data(golub)
X1 = golub[,which(golub.cl==0)]
X2 = golub[,which(golub.cl==1)]
rownames(X1) = golub.gnames[,3]
rownames(X2) = golub.gnames[,3]

# use subset of 200 genes for example
set.seed(1234)
ind = sample(1:nrow(X1),200)
X1 = X1[ind,]
X2 = X2[ind,]

######## Derive modules in group 1
modules = find_modules(X1,X2,cluster_group=1)
modules$num_modules # number of modules estimated by BIC (modules with < 3 genes are excluded)
ngm = unlist(lapply(modules$group1_modules, ncol)) # number of genes per module
summary(ngm)


######## test modules for differential co-expression
testmods = test_modules(group1_modules = modules$group1_modules,group2_modules = modules$group2_modules)
# View(testmods$pvalues)
# View(testmods$qvalues)

which(testmods$pvalues$PND6 <= 0.05)
which(testmods$qvalues$PND6 <= 0.05)

# use parallel computing:
# testmods = test_modules(group1_modules = modules$group1_modules,
#                        group2_modules = modules$group2_modules,
#                        parallel=TRUE,
#                        cores=4)


######## plot module 5
heat = corrheatmap(modules$group1_modules[[5]],modules$group2_modules[[5]])
# plot(heat$both_plots)
# plot(heat$group1_plot)
# plot(heat$group2_plot)

arbet003/discoMod documentation built on Dec. 19, 2021, 4:36 a.m.