I found a related but not quite the same quesiton in this post (Using WGCNA egiengenes for modules in correlated networks) as mine is specific about findind the hub genes in modules of interest. The hub gene I refer to in this post follows the definition in 7d of this turorial (https://labs.genetics.ucla.edu/horvath/CoexpressionNetwork/Rpackages/WGCNA/Tutorials/Simulated-07-Membership.pdf): abs(GS) and abs(MM) above a certain threshold.
Following this criteria I found some of the hubs turned out coming from modules outside the module of interest. For example, if I want to find hub genes in module brown by requesting hub=abs(GS) > .75 & abs(datKME$MM.brown)> .75, a few hubs that met the requirement belong to non-brown module. This is due to the fact that their MMbrown was smaller than -0.75.
My understanding for MM is the corelation of the gene expression and the module eigengene. Finding hubs with negative corelation with eigengene seems counter-intuitive to me. Maybe I miss some part about this concept. Any input/clarification would be greatly appreciated.