Different results from blockwiseModules and blockwiseConsensusModules individual TOM
0
0
Entering edit mode
Levi • 0
@5c8480a0
Last seen 8 weeks ago
United States

I am trying to construct a consensus WGCNA network from RNAseq data on three brain regions in two lines of zebrafish. My main hitch is with different outputs for the individual networks. I have constructed the consensus network using blockwiseConsensusModules (as shown in the code below) and saved the individual TOMs to disk. I then constructed individual networks for each brain region (one example is shown below) from each of the individual TOMs.

However, when I construct the individual networks from blockwiseModules using the same settings as before, I get a different network. I am thinking that something I am doing after TOM calculation in the first part is to blame, but I can't figure out what. The tutorial pages for this have been taken down so it's difficult for me to reference what is going on. My rationale for using the individual TOMs from blockwiseConsensusModules is that I figure that way the consensus network matches the dynamics in the individual networks as closely as possible, but the differences in these code outputs is making me think otherwise.

Can I relate the individual networks calculated with blockwiseModules to the consensus network? What is going on? Any help would be so appreciated.

net_cons <- blockwiseConsensusModules(multiExpr, TOMType = "signed", networkType = "signed", power = softPower, minModuleSize = 30, deepSplit = 2, pamRespectsDendro = FALSE, mergeCutHeight = 0.25, numericLabels = TRUE, verbose = 5, maxBlockSize = 25000, corType = "bicor", minKMEtoStay = 0)
save(net_cons, file = "net_cons.RData")

#consensus modules
table(net_cons$colors)
length(table(net_cons$colors))

#Bringing in calculated TOMS

load("individualTOM-Set1-Block1.RData")
tomDS_Dm <- tomDS
head(tomDS_Dm)
dim(tomDS_Dm)

# We like large modules, so we set the minimum module size relatively high:
minModuleSize = 30

# Call the hierarchical clustering function
geneTreeDm = hclust(as.dist(1-tomDS_Dm), method = "average")
# Plot the resulting clustering tree (dendrogram)
sizeGrWindow(12,9)
plot(geneTreeDm, xlab="", sub="", main = "Dm Gene clustering on TOM-based dissimilarity",
labels = FALSE, hang = 0.04)

# Module identification using dynamic tree cut:
dynamicModsDm = cutreeDynamic(dendro = geneTreeDm, distM = as.matrix(1-tomDS_Dm),
deepSplit = 2, pamRespectsDendro = FALSE,
minClusterSize = minModuleSize)
table(dynamicModsDm)

# Convert numeric labels into colors
dynamicColorsDm = labels2colors(dynamicModsDm)
table(dynamicColorsDm)
# Plot the dendrogram and colors underneath
sizeGrWindow(8,6)
plotDendroAndColors(geneTreeDm, dynamicColorsDm, "Dynamic Tree Cut",
dendroLabels = FALSE, hang = 0.03,
addGuide = TRUE, guideHang = 0.05,
main = "Gene dendrogram and module colors Dm")

# Calculate eigengenes
MEListDm = moduleEigengenes(as.data.frame(multiExpr[[1]]), colors = dynamicColorsDm)
MEsDm = MEListDm$eigengenes
# Calculate dissimilarity of module eigengenes
MEDissDm= 1-cor(MEsDm);
# Cluster module eigengenes
METreeDm = hclust(as.dist(MEDissDm), method = "average");
# Plot the result
sizeGrWindow(7, 6)
plot(METreeDm, main = "Clustering of module eigengenes Dm",
xlab = "", sub = "")

MEDissThres = 0.25
# Plot the cut line into the dendrogram
abline(h=MEDissThres, col = "red")
# Call an automatic merging function
mergeDm = mergeCloseModules(as.data.frame(multiExpr[[1]]), dynamicColorsDm, cutHeight = MEDissThres, verbose = 3)
# The merged module colors
mergedColorsDm = mergeDm$colors;
# Eigengenes of the new merged modules
mergedMEsDm = mergeDm$newMEs

sizeGrWindow(12, 9)
pdf(file = "geneDendro-3_Dm.pdf", wi = 9, he = 6)
plotDendroAndColors(geneTreeDm, cbind(dynamicColorsDm, mergedColorsDm),
c("Dynamic Tree Cut", "Merged dynamic"),
dendroLabels = FALSE, hang = 0.03,
addGuide = TRUE, guideHang = 0.05)
dev.off()
table(mergedColorsDm)

# Rename to moduleColors
moduleColorsDm = mergedColorsDm
# Construct numerical labels corresponding to the colors
colorOrder = c("grey", standardColors(50));
moduleLabelsDm = match(moduleColorsDm, colorOrder)-1;
MEsDm = mergedMEsDm
length(table(moduleColorsDm)) #51
table(moduleColorsDm)
sum(table(moduleLabelsDm))

This is the module output I get from this approach: enter image description here


net_Dm <- blockwiseModules(multiExpr[[1]]$data, TOMType = "signed", networkType = "signed", power = softPower, minModuleSize = 30, deepSplit = 2, pamRespectsDendro = FALSE, mergeCutHeight = 0.25, numericLabels = TRUE, verbose = 5, maxBlockSize = 25000, corType = "bicor", minKMEtoStay = 0)
save(net_Dm, file = "net_Dm.RData")
table(net_Dm$colors)

This is the output I get from this approach:

enter image description here

Session Info below:

R version 4.3.0 (2023-04-21 ucrt) Platform: x86_64-w64-mingw32/x64 (64-bit) Running under: Windows 10 x64 (build 19045)

Matrix products: default

locale: 1 LC_COLLATE=English_United States.utf8 LC_CTYPE=English_United States.utf8 LC_MONETARY=English_United States.utf8 [4] LC_NUMERIC=C LC_TIME=English_United States.utf8

time zone: America/Chicago tzcode source: internal

attached base packages: 1 stats graphics grDevices utils datasets methods base

other attached packages: 1 lubridate_1.9.2 forcats_1.0.0 stringr_1.5.0 dplyr_1.1.2 purrr_1.0.1
[6] readr_2.1.4 tidyr_1.3.0 tibble_3.2.1 ggplot2_3.4.2 tidyverse_2.0.0
[11] WGCNA_1.72-1 fastcluster_1.2.3 dynamicTreeCut_1.63-1

loaded via a namespace (and not attached): 1 tidyselect_1.2.0 blob_1.2.4 Biostrings_2.68.1 bitops_1.0-7
[5] fastmap_1.1.1 RCurl_1.98-1.12 digest_0.6.31 rpart_4.1.19
[9] timechange_0.2.0 lifecycle_1.0.3 cluster_2.1.4 survival_3.5-5
[13] KEGGREST_1.40.0 RSQLite_2.3.1 magrittr_2.0.3 compiler_4.3.0
[17] rlang_1.1.1 Hmisc_5.1-0 tools_4.3.0 utf8_1.2.3
[21] yaml_2.3.7 data.table_1.14.8 knitr_1.43 htmlwidgets_1.6.2
[25] bit_4.0.5 withr_2.5.0 foreign_0.8-84 BiocGenerics_0.46.0
[29] nnet_7.3-19 grid_4.3.0 stats4_4.3.0 preprocessCore_1.62.1
[33] fansi_1.0.4 colorspace_2.1-0 GO.db_3.17.0 scales_1.2.1
[37] iterators_1.0.14 cli_3.6.1 rmarkdown_2.22 crayon_1.5.2
[41] generics_0.1.3 rstudioapi_0.14 tzdb_0.4.0 httr_1.4.6
[45] DBI_1.1.3 cachem_1.0.8 zlibbioc_1.46.0 splines_4.3.0
[49] parallel_4.3.0 impute_1.74.1 AnnotationDbi_1.62.1 XVector_0.40.0
[53] matrixStats_0.63.0 base64enc_0.1-3 vctrs_0.6.2 Matrix_1.5-4.1
[57] hms_1.1.3 IRanges_2.34.0 S4Vectors_0.38.1 bit64_4.0.5
[61] Formula_1.2-5 htmlTable_2.4.1 foreach_1.5.2 glue_1.6.2
[65] codetools_0.2-19 stringi_1.7.12 gtable_0.3.3 GenomeInfoDb_1.36.0
[69] munsell_0.5.0 pillar_1.9.0 htmltools_0.5.5 GenomeInfoDbData_1.2.10 [73] R6_2.5.1 doParallel_1.0.17 evaluate_0.21 lattice_0.21-8
[77] Biobase_2.60.0 png_0.1-8 backports_1.4.1 memoise_2.0.1
[81] Rcpp_1.0.10 gridExtra_2.3 checkmate_2.2.0 xfun_0.39
[85] pkgconfig_2.0.3

RNASeqData WGCNA consensusnetworks • 482 views
ADD COMMENT
0
Entering edit mode

Why do you think the modules for an individual block should be remotely similar to the consensus modules generated using all the blocks? Did you read something somewhere that lead you to think that was likely?

ADD REPLY
0
Entering edit mode

Hi James,

I do not think this - my understanding is that the consensus modules are calculated from the quantile of individual topological overlap matrices in each data set. Thus, my expectation is that the individual topological overlap matrices calculated (and stored) from my blockwiseConsensusModules for each of my individual data sets should be the same as the individual topological overlap matrices from blockwiseModules for each individual data set. Am i misunderstanding something in the calculation of these networks? My primary problem here is that the individual TOMs my consensus network is calculated from differ from the individual TOMs I calculated wit ha separate function.

ADD REPLY

Login before adding your answer.

Traffic: 382 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6