I get the same modules but different number of genes WGCNA

Question

[debitboro]

Hi Biostars,

I've a question about WGCNA for co-expression network analysis. When I try to merge the modules basis on their expression profile using the following code:

merge = mergeCloseModules(datExpr, dynamicColors, cutHeight = 0.2, verbose = 3)

I get the merged module with the number of genes inside each module using the code to display those modules:

merge$colors

I get the same modules like above, BUT with different number of genes in each module, when I calculate the gene-module membership using:

geneModuleMembership = as.data.frame(cor(datExpr, MEs, use = "p"));
listModule <- c() #new list of modules
for (i in 1: numberGenes){
listModule <- c(listModule, names(which.max(abs(geneModuleMembership[i,]))))
}
table(listModule)

Do you have an idea why I get different number of genes ?

Thanks in advance

[GenoMax]

Different by how much? Many programs produce non-deterministic output unless they allow for a way to produce deterministic output.

[debitboro]

Thank you genomax2 for you response.

I post here the results for the two sets:

 table(merge$color)

    black          blue         brown          cyan     darkgreen
      345          1509           679           580           236
 darkgrey    darkorange       darkred darkturquoise         green
      103            92           294           418           656
greenyellow          grey   lightyellow          pink        purple
          689           161          2330           691           294
    royalblue        salmon
          143           228

and

table(listModule)

    MEblack          MEblue         MEbrown          MEcyan     MEdarkgreen
        546             909             775             562             307
 MEdarkgrey    MEdarkorange       MEdarkred MEdarkturquoise         MEgreen
        381             318             446             435             764
MEgreenyellow          MEgrey   MElightyellow          MEpink        MEpurple
            668             335            1389             550             333
    MEroyalblue        MEsalmon
            338             392

Answer 1

[John Ma]

The problem is why you used geneModuleMembership = as.data.frame(cor(datExpr, MEs, use = "p")) in the first place. WGCNA doesn't use Pearson correlation to assign genes to modules, so for all facts and purposes you're doing a different method of clustering. Certainly mergeCloseModules uses correlation, but that involves the correlation between module eigengenes. Because of this, there's no guarantee that any individual gene's correlation with its own ME is the smallest.

[debitboro]

Thank you MMa,

Yeah, I think I'll use the gene-module membership measure calculated from mergeCloseModules() instead to use the Pearson's correlation

[GenoMax]

If this answer helped you please up-vote/accept (green check mark) to provide closure for the thread.