WGCNA 分析确定重要模块之后，如何找模块内的hubgene？
想找hubgene,必须先理解其概念：网络中处于核心位置的点（基因）。
在WGCNA中，由于进行了无尺度化处理，具体表现为少数基因（图中红点）与众多基因建立联系，其余大多数基因则没有联系，数值化即 hubgene 为同一网络中具有更高的 connectivity或 degree。
然而抛开网络本身，假如已有先验信息例如知道某个基因具有重要功能，这种也应该作为hubgene在网络中重点关注,即使它并没有太高的connectivity或degree。
power=6
moduleColors = labels2colors(net_power$colors) #net_power为一步法构建出的对象
1. 基于连通度connectivity，等同于常规网络中的degree。
由于加权共表达网络用数值度量了点点相关性，所以网络中某点的连通度，相当于与其余所有点相关性的加和。
WGCNA中又将网络划分成了不同modules,所以连通度又可以分为整个网络的 kTotal ,和某个模块的 kWithin 。
kTotal：某基因与网络中其余所有有边相连基因的 degree 之和。
kWithin：某基因与同一模块中其余所有有边相连基因的 degree 之和。对于已知模块，kWithin越大，越可能为 hubgene 。
kOut=kTotal-kWithin, 此基因假设在黄色模块中，则 kOut为 与黄色模块之外所有模块的连通度之和。
kDiff=kIn-kOut=2*kIN-kTotal，评估此基因在某模块中连通性的强弱，不太常用。
计算上有两个策略
参考：https://www.rdocumentation.org/packages/WGCNA/versions/1.70-3/topics/intramodularConnectivity
#第一种： 先计算邻接矩阵，再计算连通度，推荐此种,邻接矩阵(Adjacency matrix)指基因和基因之间的加权相关性值取power次方即无尺度化之后构成的矩阵。
adjacency = abs(cor(datExpr,use="p"))^power #datExpr为表达矩阵，power请与一步法中的软阈值保持一致
kIM <- intramodularConnectivity(adjacency, moduleColors)
#第二种：直接输入表达矩阵计算
kIM <- intramodularConnectivity.fromExpr(datExpr, colors, power = power)
#查看
head(kIM)
# kTotal kWithin kOut kDiff
# Gene1 90 80 10 70
# Gene2 9 8 1 7
2. WGCNA自带一个函数，可以提取每个模块中连通度最高的基因
参考：https://www.rdocumentation.org/packages/WGCNA/versions/1.70-3/topics/chooseTopHubInEachModule
hub = chooseTopHubInEachModule(datExpr,moduleColors, omitColors = "grey", power = power)
3. 有的文献中认为，hubgene 应该在某一模块中，与核心性状关联且与此模块也同样关联,比如：|GS|>.2&|MM|>.8 。
GS 为 gene significance 缩写,反映某基因表达量与对应表型数值的相关性，调用 cor 函数计算相关性，绝对值越大，相关性越大，越趋近0，越不相关。
MM 为 module membership 缩写，反映某基因表达量与模块特征值(module eigengenes)的相关。同样调用 cor 函数计算相关性，绝对值越大，相关性越大，越趋近0，越不相关。
ps: 模块特征值 (module eigengenes) 是指基因和样本矩阵进行PCA分析后，每个模块的第一主成分即PC1,是此模块特性的高度代表,可以理解为一个模块就是一个超级代表基因。
#3.1.首先计算模块特征值(module eigengenes)
MEs = moduleEigengenes(datExpr, moduleColors)$eigengenes
#3.2.计算module membership即MM
datKME =signedKME(datExpr, MEs, outputColumnName="MM.")
#3.3.gene significance即GS
GS = as.data.frame(cor(datExpr, datTraits, use = "p")) #datTraits为目标性状的矩阵文件
#3.4.循环提取每个模块的hubgene
MMname = colnames(datKME)
for (mm in MMname){
FilterGenes =abs(GS)> .2 & abs(datKME$mm)>.8
hubgenes = dimnames(data.frame(datExpr))[[2]][FilterGenes]
}
需要注意的是这里并没有添加显著检验，正确的做法应该是同时针对pvalue进行控制，比如：|GS|>.2 & |MM|>.8 & pvalue<0.05
4. 同时，WGCNA官方建议结合 GS 与 kWithin 进行筛选
Finding genes with high gene significance and high intramodular connectivity in interesting modules,比如: kWithin > 30 & |GS| >0.5
5. 实际操作中，由于导出到 cytoscape 等可视化工具时边的数量格外多，就需要额外进行一个过滤步骤。此时常用的是针对 TOM 值做一个阈值的筛选，例如只保留 TOM>0.8 的基因关系对，此时计算模块内基因的 degree ，值越高，则越可能为hubgene 。
需要注意的时，此时的网络中，由于设定了阈值0.8，则默认剩余基因只要存在边的相连即为存在相关性，则这里的点与点之间数值价值就不大了（不是没有价值），此时degree也可计算为某基因，与之相连边的数目总和，同样也是相连点的数目总和。
ps: TOM (Topological overlap matrix)：把邻接矩阵转换为拓扑重叠矩阵，是另一种距离矩阵。TOM值不止考虑两个基因的关系，还考虑依赖于这一对基因对存在的其它基因的影响，并将之量化，充分考虑了基因表达调控网络的复杂性，也更能代表两个基因之间的实际关系。
#5.1. 载入TOM值
# 如果采用分步计算，或设置的blocksize>=总基因数，直接load计算好的TOM结果
# 否则需要再计算一遍，比较耗费时间
TOM = TOMsimilarityFromExpr(datExpr, power = power)
#假如之前一步法中设定一个block计算，则可以直接载入节省时间
# load(net_power$TOMFiles[1], verbose=T)
# TOM <- as.matrix(TOM)
#5.2. 设定阈值并输出文件
geneid_allnet <- names(datExpr)
MEs = moduleEigengenes(datExpr, moduleColors)$eigengenes
modNames <- substring(names(MEs),3)
TOMcutoff <- 0.8 #实际操作中我建议设定一个小的值，到cytoscape中再进行个性化调整
#循环输出所有模块
for (mod in 1:nrow(table(moduleColors)))
{
modules = names(table(moduleColors))[mod]
# Select module probes
probes = names(datExpr)
inModule = (moduleColors == modules)
modProbes = probes[inModule]
modGenes = modProbes
# Select the corresponding Topological Overlap
modTOM = TOM[inModule, inModule]
dimnames(modTOM) = list(modProbes, modProbes)
# Export the network into edge and node list files Cytoscape can read
cyt = exportNetworkToCytoscape(modTOM,
edgeFile = paste("CytoscapeInput-edges-", modules , ".txt", sep=""),
nodeFile = paste("CytoscapeInput-nodes-", modules, ".txt", sep=""),
weighted = TRUE,
threshold = TOMcutoff,
nodeNames = modProbes,
altNodeNames = modGenes,
nodeAttr = moduleColors[inModule])
}
#5.3，接下来针对CytoscapeInput-edges-*txt文件，进行统计，看哪个基因degree（边的数目总和）高，则此基因为hubgene
6. 除开以上连通度，度的衡量和过滤，hubgene的筛选还要采取灵活的方法，结合自己的研究目的及更多的数据进行。
例如加入转录因子信息，由于转录因子的特性，往往成为hubgene。
或者某合成通路中的限速酶，也可能会成为模块内的hubgene。
假如这个限速酶连通度并不高，那么跟它有关系的点里有没有连通度高的呢？如果这个点又是一个转录因子，那可就太好了。
