Question: Wgcna blockwisemodules: lower power results in modules with positive and negative kme values.
0
gravatar for ductylerlv
23 days ago by
ductylerlv0
ductylerlv0 wrote:

Hello All,

Our lab is getting some interesting results whereby we have modules that have both positive and negative kme values in a signed network. If a higher power is used this goes away but the lower power already surpasses the R2 and k.means thresholds. Any insight would be helpful.

Best,
Duc

wgcna • 155 views
ADD COMMENTlink modified 6 days ago by edammer0 • written 23 days ago by ductylerlv0
2

can you post the code ?what you had done?

ADD REPLYlink written 23 days ago by krushnach80390

I was the one who encountered negative and non-correlated kME values for specific network nodes in WGCNA with the following input and code. The soft threshold criterion for scale free topology was met at power, beta, 7.6. We chose 8.5 to reduce noise further. SFT power vs. R² for this dataset

And the resulting protein network has module(s) like the below violet module, where we never usually see kMEintramodule<0.20, but here we do: lower kME violet module members

Clearly some of the lower violet module members should go to orangered4 with kMEorangered4>0.7.

Interestingly, some orangered non-hub members suffer from anticorrelation to their own ME, precisely what should not be seen in signed network modules: lower kME orangered4 module members

Full code to reproduce the problem follows:

rm(list=ls())
options(stringsAsFactors=FALSE)
rootdir <- "drive:/folder/location/" # This is the folder containing all of the analysis scripts, input, and output for this project

library(WGCNA)
allowWGCNAThreads()

setwd(rootdir)
cleanDat<-read.table(file=paste0(rootdir,"/final_cleanDat_11240x243_BootAgeSexPMIregr_DxNumericProtected_regr_Sept20_2018.txt"),sep="\t")

allowWGCNAThreads()
powers <- seq(4,15,by=0.5)
tableSFT<-pickSoftThreshold(t(cleanDat),powerVector=powers,corFnc="bicor",blockSize=15000,verbose=3,networkType="signed")[[2]] #,corFnc="bicor"
plot(tableSFT[,1],tableSFT[,2],xlab="Power (Beta)",ylab="SFT R²")

# pickSoftThreshold: calculating connectivity for given powers...
#   ..working on genes 1 through 11240 of 11240
#   Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k.
#1    4.0    0.715 -8.28          0.959  810.00    790.00 1150.0
...
#8    7.5    0.903 -4.52          0.991  104.00     94.50  254.0
#9    8.0    0.909 -4.23          0.994   79.00     70.90  211.0
...

powers <- seq(7.6,7.9,by=0.1)
tableSFT<-pickSoftThreshold(t(cleanDat),powerVector=powers,corFnc="bicor",networkType="signed")
#plot(tableSFT[,1],tableSFT[,2],xlab="Power (Beta)",ylab="SFT R²")
tableSFT$powerEstiamte
#[1] 7.6

#  Power SFT.R.sq slope truncated.R.sq mean.k. median.k. max.k.
#1   7.6    0.905 -4.46          0.992    98.3      89.2    245
#2   7.7    0.904 -4.41          0.992    93.1      84.2    236
#3   7.8    0.905 -4.35          0.993    88.1      79.5    227
#4   7.9    0.907 -4.29          0.993    83.5      75.0    219

## Run an automated network analysis (WGCNA 1.64 Rx64 v3.5.1 WINDOWS)

power=8.5; #first tried: 7.6
minModSize=30
enforceMMS=FALSE

net <- blockwiseModules(t(cleanDat),power=power,deepSplit=4,minModuleSize=minModSize,
                        mergeCutHeight=0.07, TOMDenom="mean", #detectCutHeight=0.9999,
corType="bicor",networkType="signed",pamStage=TRUE,pamRespectsDendro=TRUE,reassignThresh=0.05,
                        verbose=3,saveTOMs=FALSE,maxBlockSize=12000)

# Summary Module Table
table(net$colors)["grey"] #3179/11,240 are grey/unassigned.

nModules<-length(table(net$colors))-1
modules<-cbind(colnames(as.matrix(table(net$colors))),table(net$colors))
orderedModules<-cbind(Mnum=paste("M",seq(1:nModules),sep=""),Color=labels2colors(c(1:nModules)))
modules<-modules[match(as.character(orderedModules[,2]),rownames(modules)),]
as.data.frame(cbind(orderedModules,Size=modules))

#                Mnum           Color Size
#turquoise         M1       turquoise  643
#blue              M2            blue  494
#brown             M3           brown  423
...

tmpMEs <- MEs <- net$MEs
colnames(tmpMEs) <- paste("ME",colnames(MEs),sep="")
kMEdat <- signedKME(t(cleanDat), tmpMEs, corFnc="bicor")

write.table(cbind(rownames(cleanDat),net$colors,kMEdat),file=paste(rootdir,"/ModuleAssignments_TOMdenomMEAN_netPwr",power,"_MergeHeight0.07_PAMstageTRUE.txt",sep=""),sep="\t")

Formatted ModuleAssignments (kME table) can be found here: Excel Binary kME Table for above network output

Input for cleanDat log2(abundance) matrix loaded in above code is here: ZIP with log2(abundance) matrix text tab-separated file

...explicit use of minKMEtoStay=0.30 has never been necessary before.

ADD REPLYlink written 6 days ago by edammer0

I should add that minKMEtoStay=0.30 is a default parameter/setting: blockwiseModules in WGCNA v1.64 default parameters

Rerunning the above blockwiseModules() function with this parameter added to the others gives identical modules.

My sessionInfo follows:

> sessionInfo()
R version 3.5.1 (2018-07-02)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 17134)

Matrix products: default

locale:
[1] LC_COLLATE=English_United States.1252  LC_CTYPE=English_United States.1252   
[3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C                          
[5] LC_TIME=English_United States.1252    

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

other attached packages:
[1] WGCNA_1.64-1          fastcluster_1.1.25    dynamicTreeCut_1.63-1

loaded via a namespace (and not attached):
 [1] Biobase_2.40.0        bit64_0.9-7           splines_3.5.1         foreach_1.4.4        
 [5] Formula_1.2-3         assertthat_0.2.0      stats4_3.5.1          latticeExtra_0.6-28  
 [9] blob_1.1.1            fit.models_0.5-14     robustbase_0.93-2     impute_1.54.0        
[13] pillar_1.3.0          RSQLite_2.1.1         backports_1.1.2       lattice_0.20-35      
[17] glue_1.3.0            digest_0.6.17         RColorBrewer_1.1-2    checkmate_1.8.5      
[21] colorspace_1.3-2      htmltools_0.3.6       preprocessCore_1.42.0 Matrix_1.2-14        
[25] plyr_1.8.4            pcaPP_1.9-73          pkgconfig_2.0.2       purrr_0.2.5          
[29] GO.db_3.6.0           mvtnorm_1.0-8         scales_1.0.0          htmlTable_1.12       
[33] tibble_1.4.2          IRanges_2.14.11       ggplot2_3.0.0         nnet_7.3-12          
[37] BiocGenerics_0.26.0   lazyeval_0.2.1        survival_2.42-6       magrittr_1.5         
[41] crayon_1.3.4          memoise_1.1.0         doParallel_1.0.11     MASS_7.3-50          
[45] foreign_0.8-71        tools_3.5.1           data.table_1.11.4     matrixStats_0.54.0   
[49] stringr_1.3.1         S4Vectors_0.18.3      munsell_0.5.0         cluster_2.0.7-1      
[53] AnnotationDbi_1.42.1  bindrcpp_0.2.2        compiler_3.5.1        rlang_0.2.2          
[57] grid_3.5.1            iterators_1.0.10      rstudioapi_0.7        htmlwidgets_1.2      
[61] robust_0.4-18         base64enc_0.1-3       gtable_0.2.0          codetools_0.2-15     
[65] DBI_1.0.0             rrcov_1.4-4           R6_2.2.2              gridExtra_2.3        
[69] knitr_1.20            dplyr_0.7.6           bit_1.1-14            bindr_0.1.1          
[73] Hmisc_4.1-1           stringi_1.2.4         parallel_3.5.1        Rcpp_0.12.18         
[77] rpart_4.1-13          acepack_1.4.1         DEoptimR_1.0-8        tidyselect_0.2.4     
>
ADD REPLYlink written 6 days ago by edammer0
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1706 users visited in the last hour