WGCNA的入门和进阶：（二）WGCNA的R包实现
本文使用window本地Rstudio实现WGCNA的全面分析，附代码和实例供大家学习和交流。
WGCNA包的下载安装参考官网。https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/index.html
install.packages("BiocManager")
BiocManager::install("WGCNA")
femData <- read.csv("gene.csv", sep=',',header=TRUE) #加载数据
datExpr = as.data.frame(t(log2(femData[,-c(1,2)]+1))) #输入的基因表达量是TPM归一化的结果，此处对表达量再进行一次log2的转换。因为我输入的femData是样本名为列，基因名为行，所以此处使用t()进行转换，注意：得到的datExpr一定是基因名在列，样本名在行，一定不能搞混，要不然后面的所有分析都是错的。
names(datExpr) = femData$gene_id #列标签添加基因名称
rownames(datExpr) = names(femData[,-c(1,2)]) #行标签添加为样本名称
gsg = goodSamplesGenes(datExpr, verbose = 3) #检测缺失值
gsg$allOK #结果为TRUE，则所有选定基因都用于后续WGCNA
datExpr=datExpr[gsg$goodGenes] #如果gsg$allOK结果为FALSE，则后续选择好的基因用于WGCNA。
三．选择合适的软阈值β进行后续的基因共表达网络构建（参考官方说明书：https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/faq.html）
powers = c(c(1:10), seq(from = 12, to=30, by=2))
sft=pickSoftThreshold(datExpr,dataIsExpr = TRUE,
corOptions = list(use = 'p'),
networkType = 'unsigned') #目的是为了帮助用户选择一个合适的软阈值进行网络构建。
powers = c(c(1:10), seq(from = 12, to=30, by=2))
sft=pickSoftThreshold(datExpr,dataIsExpr = TRUE,
corOptions = list(use = 'p'),
networkType = 'unsigned') #目的是为了帮助用户选择一个合适的软阈值进行网络构建。
软阈值的选择标准（以下三个标准选择其中一个即可）
1、绘制软阈值和Scale Free Topology Model Fit, signed R^2的关系图，取使得Scale Free Topology Model Fit, signed R^2>0.8（或者自己卡一个阈值）的Soft Threshold作为软阈值。
plot(x=sft$fitIndices[,1],
y=-sign(sft$fitIndices[,3])*sft$fitIndices[,2],
xlab="Soft Threshold (power)",
ylab="Scale Free Topology Model Fit, signed R^2",type="n",
main = paste("Scale independence"))
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
labels=powers,cex=0.9,col="red")
2、绘制软阈值和 Mean Connectivity的关系图，取使得 Mean Connectivity <=100（或者自己卡一个阈值）的Soft Threshold作为软阈值。
plot(sft$fitIndices[,1], sft$fitIndices[,5],
xlab="Soft Threshold (power)", ylab="Mean Connectivity", type="n",
main = paste("Mean connectivity"))
text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=0.9,col="red")
如果无向网络在power小于15或有向网络power小于30内，没有一个power值可以使Scale Free Topology Model Fit, signed R^2达到0.8或Mean Connectivity降到100以下，可能是由于部分样品与其他样品差别太大造成的。这可能由批次效应、样品异质性或实验条件对表达影响太大等造成，可以通过绘制样品聚类查看分组信息、关联批次信息、处理信息和有无异常样品。如果这确实是由有意义的生物变化引起的，也可以使用经验power值。
softPower=3 #根据上述规则选择合适的softPower进行以下的加权相关性分析。
TOM = TOMsimilarityFromExpr(datExpr,
power = softPower) # TOM (Topological overlap matrix)：把基因矩阵转换为拓扑重叠矩阵，以降低噪音和假相关，获得新的基因距离矩阵，这个信息用于构建网络或绘制TOM图。
colnames(TOM) =rownames(TOM) = femData$gene_id #对构建好的TOM行列添加基因名称，得到的TOM矩阵则是两两基因加权的相关系数。
五．识别基因集：基于基因加权相关性得到的TOM矩阵，进行层级聚类分析，并根据设定标准切分聚类结果，获得不同的基因模块，用聚类树的分枝和不同颜色表示。
geneTree = hclust(as.dist(dissTOM),method='average') #对加权的基因矩阵进行聚类。
minModuleSize = 50 # 设置每个module中的最少基因个数为50。
dynamicMods = cutreeDynamic(dendro = geneTree,
pamRespectsDendro = FALSE,
minClusterSize = minModuleSize) #使用dynamic tree cut来识别基因集。
table(dynamicMods) #给出模块标签和每个模块的大小。label 0保留的为unasunsigned基因。
dynamicColors = labels2colors (dynamicMods) #在树状图下绘制模块颜色分配。注:灰色为unasunsigned基因。
六．计算Module Eigengenes（相关概念在上一专题进行具体阐述）并合并相似模块
MEList = moduleEigengenes (datExpr, colors = dynamicColors) #按照模块计算每个module的ME（也就是该模块的第一主成分）
MEDiss = 1-cor(MEs); # 计算根据模块特征向量基因计算模块相异度
METree = hclust(as.dist(MEDiss), method = "average"); #对不同的模块进行聚类
plotEigengeneNetworks(MEs, "Eigengene adjacency heatmap",
plot(METree, main = "Clustering of module eigengenes", xlab = "", sub = "")
abline(h=0.99, col="red")
merge_modules = mergeCloseModules(datExpr, dynamicColors, cutHeight = 0.1, verbose = 3)
mergedColors = merge_modules$colors; # 合并后的颜色：
mergedMEs = merge_modules$newMEs; # 新模块的特征向量MEs：
pdf("module cluster.pdf")
plotDendroAndColors(geneTree, cbind(dynamicColors, mergedColors),
c("Dynamic Tree Cut", "Merged dynamic"),
dendroLabels = FALSE, hang = 0.03,
addGuide = TRUE, guideHang = 0.05)
3、为了后续对指定模块中的基因进行下游分析，在此我们将每个模块的基因提取出来。
module_colors=unique(mergedColors)
for (color in module_colors){
module=gene.names[which(mergedColors==color)]
group1=rep(color,times=length(module))
module_gene=c(module_gene,module)
my_module=data.frame("colour"=cl,"gene"=module_gene)
write.csv(my_module,"Modules_gene.csv")
4、为了整体观察基因间加权的相关性，在此我们去掉那些grey module的基因（未聚类到任何一个module），并对剩余基因绘制TOM图。
hier1=flashClust(as.dist(diss1), method='average')
geneTree = hclust(as.dist(diss1), method = "average");
#set the diagonal of the dissimilarity to NA
pdf('Modules_heatmap.pdf')
TOMplot(diss1, hier1, as.character(mergedColors))
write.csv(diss1,"dissTOM.txt")
七．Eigengene significance 的计算（也可以当作Module significance，计算公式稍有区别，但表示的意思是一样的）
首先对样本进行分组，注意，在上面我们进行基因间的加权相关分析的时候是没有关注样本顺序的，这并不会影响对基因的聚类。但是，在后面的分析中，样本分组信息一定要对应样本名称，否则分析出来的结果就有偏差。
group=factor(rep(c("goup1","goup2","goup3"),each=15),levels = c("goup1","goup2","goup3")) #根据样本信息进行分组。
datTraits = data.frame(samples=rownames(datExpr),subtype=group)
design=model.matrix(~0 + datTraits$subtype)
colnames(design)=levels(factor(datTraits$subtype))
MEs0 = moduleEigengenes(datExpr, mergedColors)$eigengenes。#计算合并后的模块的第一主成分ME。
MEs = orderMEs(MEs0); #不同颜色的模块的ME值矩阵。
moduleTraitCor = cor(MEs, design , use = "p") #计算不同模块和样本性状的相关性。
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples) #计算不同模块和样本性状相关性的p值。
textMatrix = paste(signif(moduleTraitCor, 2), "\n(",signif(moduleTraitPvalue, 1), ")", sep = "");
dim(textMatrix) = dim(moduleTraitCor)
pdf('group_Modules_heatmap.pdf')
labeledHeatmap(Matrix = moduleTraitCor,
xLabels = colnames(design),
yLabels = row.names(moduleTraitCor),
ySymbols = row.names(moduleTraitCor),
colors = blueWhiteRed(50),
main = paste("Module-trait relationships"))
八．Module Membership MM的计算
#names (colors) of the modules
modNames = substring(names(MEs), 3)
geneModuleMembership = as.data.frame(cor(datExpr, MEs, use = "p")); #计算基因与module的相关性。
MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples)) #计算基因与module相关性的显著性
names(geneModuleMembership) = paste("MM", modNames, sep="");
names(MMPvalue) = paste("p.MM", modNames, sep="");
九．Gene significance GS 的计算
geneTraitSignificance = as.data.frame(cor(datExpr, design, use = "p"));
GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples));
names(geneTraitSignificance) = paste("GS.",names(geneTraitSignificance), sep="");
names(GSPvalue) = paste("p.GS.",names(geneTraitSignificance), sep="");
如果我们鉴定出某个module的基因与性状高度相关，可以进一步观察MM与GS的关系。以blue module举例说明。
column = match(module, modNames);
moduleGenes = mergedColors==module;
pdf("Module membership vs. gene significance.pdf")
verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
abs(geneTraitSignificance[moduleGenes, 1]),
xlab = paste("Module Membership in", module, "module"),
ylab = "Gene significance for time",
main = paste("Module membership in", module, "module vs. gene significance\n"),
cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
十．导出网络到edge and node list文件，可以使用Cytoscape读取，表现为同一module中基因直接的相关性。
for (color in module_colors){
gene=my_module$gene[which(my_module$colour==color)]
color_gene=as.matrix( TOM[gene,gene])
cyt = exportNetworkToCytoscape(color_gene,
edgeFile = paste('CytoscapeInput-edges-', color,'.txt', sep=''),
nodeFile = paste('CytoscapeInput-node-', color,'.txt', sep=''),
weighted = TRUE, threshold = 0, nodeNames = gene, altNodeNames = gene)
connectivity=abs(cor(datExpr,use="p"))^softpower #跟TOM一样
Alldegrees=intramodularConnectivity(connectivity, mergedColors) #计算每个基因module内的连接度和module外的连接度以及total连接度。
#Alldegrees$gene=rownames(Alldegrees)
datKME=signedKME(datExpr, MEs, outputColumnName="kME_MM.")#计算MM。
GS=cor(datExpr,design,use="p") #注意，如果是design是离散变量，一定要将design放在后面
combin_data=cbind(Alldegrees,datKME,GS1)
write.csv(combin_data,"combine_hub.csv")
一般想要找某个module(以blue为例)的hub基因的话，需要对kwithin，kME_MM_blue，GS同时卡阈值，如果某个基因kwithin，kME_MM_blue，ABS(GS)的值越大，则成为hub基因的可能性越大。
(责任编辑：伏泽 微信：18520221056)
|