Entering edit mode
Hi, I have run the SampledBlockwiseModules function for WGCNA in order to subsample my data and create modules from it iteratively. I then relabeled modules in each of the resampling runs so that full and resampled modules with best overlaps have the same labels. So, I was wondering how do I create a consensus network from the relabeled module labels? What I'm trying to do is combine the results from each of the subsamples to create the most optimal and robust network.
Code should be placed in three backticks as shown below
#Resampling analysis of module stability
# Number of resampling runs
nRuns = 10
power = 12
deepSplit = 2
minModuleSize = 30
networkType = "signed"
TOMType = "unsigned"
TOMDenom = "mean"
reassignThreshold = 0
mergeCutHeight = 0.20
# Proportion of missing data. Not needed for the calculations, but useful to know.
propNA = sum(is.na(datExpr))/length(datExpr)
propNA
tmf1 = system.time ( {
mods1 = sampledBlockwiseModules(
datExpr = datExpr,
nRuns = 10,
replace = TRUE,
fraction = 0.5,
maxBlockSize = 30000,
randomSeed = 12345,
checkSoftPower = TRUE,
nPowerCheckSamples = 2000,
skipUnsampledCalculation = FALSE,
corType = "pearson",
power = 12,
networkType = "signed",
TOMType = "unsigned",
TOMDenom = "mean",
deepSplit = 2,
mergeCutHeight = 0.20,
saveTOMs = FALSE,
saveTOMFileBase = "sampled.TOM" ) } )
# Print the timing results
print(tmf1)
# Save the resampled modules
save(tmf1, mods1, file = "sampledModuleExample-mods.RData")
#Run single iteration of network construction
slowAnalysis = function(datExpr)
{
cor = stats::cor(datExpr, use = "p")
cor[cor<0] = 0
adj = cor^12
dTOM = TOMdist(adj, TOMType = "unsigned", TOMDenom = "mean")
collectGarbage()
tree = stats::hclust(as.dist(dTOM), method = "a")
labels = cutreeDynamic(tree, minClusterSize = 30, distM = dTOM, deepSplit = 2)
mergedLabels = mergeCloseModules(datExpr, labels, cutHeight = 0.20)
mergedLabels
}
tms = system.time({slowLabels = slowAnalysis(datExpr)})
print(tms)
collectGarbage()
#Visualisation of results
nGenes = ncol(datExpr)
nSamples = nrow(datExpr)
# Define a matrix of labels for the original and all resampling runs
labels = matrix(0, nGenes, nRuns + 1)
labels[, 1] = mods0[[1]]$mods$colors
# Relabel modules in each of the resampling runs so that full and resampled modules with best overlaps have
# the same labels. This is achieved by the function matchLabels.
pind = initProgInd()
for (r in 2:(nRuns+1))
{
labels[, r] = matchLabels(mods0[[r]]$mods$colors, labels[, 1])
pind = updateProgInd((r-1)/nRuns, pind)
}
# Save the results
save(labels, file = "sampledModuleExample-matchedLabels.RData")
sessionInfo( )
R version 4.2.0 (2022-04-22) -- "Vigorous Calisthenics"
Copyright (C) 2022 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> sessionInfo()
R version 4.2.0 (2022-04-22)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 20.04.3 LTS
Matrix products: default
BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
locale:
[1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8 LC_COLLATE=C.UTF-8
[5] LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8 LC_PAPER=C.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
loaded via a namespace (and not attached):
[1] backports_1.4.1 Hmisc_4.7-0 scWGCNA_1.0.0 plyr_1.8.7
[5] igraph_1.3.2 lazyeval_0.2.2 sp_1.5-0 splines_4.2.0
[9] RApiSerialize_0.1.0 listenv_0.8.0 scattermore_0.8 GenomeInfoDb_1.32.2
[13] ggplot2_3.3.6 digest_0.6.29 foreach_1.5.2 htmltools_0.5.2
[17] GO.db_3.15.0 fansi_1.0.3 checkmate_2.1.0 magrittr_2.0.3
[21] memoise_2.0.1 tensor_1.5 cluster_2.1.3 doParallel_1.0.17
[25] ROCR_1.0-11 globals_0.15.1 fastcluster_1.2.3 Biostrings_2.64.0
[29] RcppParallel_5.1.5 matrixStats_0.62.0 spatstat.sparse_2.1-1 jpeg_0.1-9
[33] colorspace_2.0-3 blob_1.2.3 ggrepel_0.9.1 xfun_0.31
[37] dplyr_1.0.9 crayon_1.5.1 RCurl_1.98-1.7 jsonlite_1.8.0
[41] progressr_0.10.1 spatstat.data_2.2-0 impute_1.70.0 stringfish_0.15.7
[45] survival_3.2-13 zoo_1.8-10 iterators_1.0.14 glue_1.6.2
[49] polyclip_1.10-0 gtable_0.3.0 zlibbioc_1.42.0 XVector_0.36.0
[53] leiden_0.4.2 future.apply_1.9.0 BiocGenerics_0.42.0 abind_1.4-5
[57] scales_1.2.0 DBI_1.1.3 spatstat.random_2.2-0 miniUI_0.1.1.1
[61] Rcpp_1.0.8.3 htmlTable_2.4.0 viridisLite_0.4.0 xtable_1.8-4
[65] reticulate_1.25 spatstat.core_2.4-4 foreign_0.8-82 bit_4.0.4
[69] preprocessCore_1.58.0 Formula_1.2-4 stats4_4.2.0 htmlwidgets_1.5.4
[73] httr_1.4.3 RColorBrewer_1.1-3 ellipsis_0.3.2 Seurat_4.1.1
[77] ica_1.0-2 pkgconfig_2.0.3 nnet_7.3-17 uwot_0.1.11
[81] deldir_1.0-6 utf8_1.2.2 dynamicTreeCut_1.63-1 tidyselect_1.1.2
[85] rlang_1.0.3 reshape2_1.4.4 later_1.3.0 AnnotationDbi_1.58.0
[89] munsell_0.5.0 tools_4.2.0 cachem_1.0.6 cli_3.3.0
[93] generics_0.1.2 RSQLite_2.2.14 ggridges_0.5.3 stringr_1.4.0
[97] fastmap_1.1.0 goftest_1.2-3 knitr_1.39 bit64_4.0.5
[101] fitdistrplus_1.1-8 purrr_0.3.4 hdWGCNA_0.1.1.9002 RANN_2.6.1
[105] KEGGREST_1.36.2 pbapply_1.5-0 future_1.26.1 nlme_3.1-157
[109] mime_0.12 compiler_4.2.0 rstudioapi_0.13 plotly_4.10.0
[113] png_0.1-7 spatstat.utils_2.3-1 tibble_3.1.7 stringi_1.7.6
[117] rgeos_0.5-9 lattice_0.20-45 Matrix_1.4-1 vctrs_0.4.1
[121] pillar_1.7.0 lifecycle_1.0.1 spatstat.geom_2.4-0 lmtest_0.9-40
[125] RcppAnnoy_0.0.19 data.table_1.14.2 cowplot_1.1.1 bitops_1.0-7
[129] irlba_2.3.5 httpuv_1.6.5 patchwork_1.1.1 latticeExtra_0.6-29
[133] R6_2.5.1 qs_0.25.3 promises_1.2.0.1 KernSmooth_2.23-20
[137] gridExtra_2.3 IRanges_2.30.0 parallelly_1.32.0 codetools_0.2-18
[141] MASS_7.3-57 assertthat_0.2.1 SeuratObject_4.1.0 sctransform_0.3.3
[145] S4Vectors_0.34.0 GenomeInfoDbData_1.2.8 mgcv_1.8-40 parallel_4.2.0
[149] grid_4.2.0 rpart_4.1.16 tidyr_1.2.0 Rtsne_0.16
[153] base64enc_0.1-3 Biobase_2.56.0 shiny_1.7.1 WGCNA_1.71