如果有一点R语言基础会发现上面代码就做了三件事 ， 选择基因 ， 过滤基因 ， 基因聚类 。为什么要做这些事情 ， 下面来一件一件说 。
选择基因中有三种方法 。其实这个部分是很灵活的 ， 三种方法都可以 ， 如果计算机资源可以 ， 那完全可以采用所有基因进行后续分析 。
过滤基因 。这相当于数据的清洗 。
基因聚类 。这一部分特别重要 。有人说可以不做 ， 但自我感觉重要 ， 因为我们后面做的都是与聚类相关的 ， 如果不做这一步直接做后面相关的聚类可信度会少很多 。
筛选软阈值 。为什么要筛选软阈值？那就需要了解一下软阈值的意义 。这是WGCNA的底层思维 。筛选到软阈值之后就能构建模块了 ， 也就是构建共表达网络 。有没有发现 ， 这是一步一步递进的 。
主要运用的函数是 ， 这个函数有必要好好研究一下 ， 源代码我会放在最后 ， 有时间我再来写篇文章 。
这个函数运行的过程我大致讲一下 ， 要好好理解才行 。
第一步
在得到一个m × n 表达矩阵（m个样本 ， n个基因）后 ， 第一步是计算两基因在多样本的表达水平相关性（例如  ， 值一般在-1 ~ 1之间） ， 从而得到表达相似度矩阵 matix ， 简记为S 。（在这里要思考一下为什么要计算表达水平相关性 ， 举一反三 ， 计算别的相关性可以吗？）
根据是否考虑相关性的正负性 ， 有两种处理方式 。
不要把sign和弄混了 。前者表示相关的正负性 ， 后者表示方向性 。使用WGCNA建立的都是网络 。
第二步
使用 将基因表达相似矩阵转为共表达邻接矩阵 。根据 的特点分为hard与soft两种 。
第三步
鉴定模块及相关分析
# 清空所有变量rm(list = ls())# 加载包library(WGCNA)# 允许多线程运行enableWGCNAThreads()# 加载表达矩阵load("data/Step01-WGCNA_input.Rda")### 选择软阈值 ###powers = c(c(1:10), seq(from = 12, to=20, by=2))# 进行网络拓扑分析sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)#β=power,就是软阈值# 可视化结果sizeGrWindow(9, 5)pdf(file = "figures/Step02-SoftThreshold.pdf", width = 9, height = 5);par(mfrow = c(1,2))#一个画板上 ， 画两个图 ， 一行两列cex1 = 0.9;# 无尺度网络阈值得选择plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],xlab="Soft Threshold (power)",#x轴ylab="Scale Free Topology Model Fit,signed R^2",type="n",#y轴main = paste("Scale independence"));#标题text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],labels=powers,cex=cex1,col="red");# 用红线标出R^2的参考值abline(h=0.90,col="red")# 平均连接度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=cex1,col="red")dev.off()# 无尺度网络检验 ， 验证构建的网络是否是无尺度网络softpower=sft$powerEstimateADJ = abs(cor(datExpr,use="p"))^softpower#相关性取绝对值再幂次k = as.vector(apply(ADJ,2,sum,na.rm=T))#对ADJ的每一列取和 ， 也就是频次pdf(file = "figures/Step02-scaleFree.pdf",width = 14)par(mfrow = c(1,2))hist(k)#直方图scaleFreePlot(k,main="Check scale free topology")dev.off()### 一步构建网络 ###net = blockwiseModules(datExpr, #处理好的表达矩阵power = sft$powerEstimate,#选择的软阈值TOMType = "unsigned", #拓扑矩阵类型 ， none表示邻接矩阵聚类 ， unsigned最常用 ， 构建无方向minModuleSize = 30,#网络模块包含的最少基因数reassignThreshold = 0, #模块间基因重分类的阈值mergeCutHeight = 0.25,#合并相异度低于0.25的模块numericLabels = TRUE, #true ， 返回模块的数字标签 false返回模块的颜色标签pamRespectsDendro = FALSE,#调用动态剪切树算法识别网络模块后 ， 进行第二次的模块比较 ， 合并相关性高的模块saveTOMs = TRUE,#保存拓扑矩阵saveTOMFileBase = "data/Step02-fpkmTOM", verbose = 3)#0 ， 不反回任何信息 ， ＞0返回计算过程# 保存网络构建结果save(net, file = "data/Step02-One_step_net.Rda")# 加载网络构建结果load(file = "data/Step02-One_step_net.Rda")# 打开绘图窗口sizeGrWindow(12, 9)pdf(file = "figures/Step02-moduleCluster.pdf", width = 12, height = 9);# 将标签转化为颜色mergedColors = labels2colors(net$colors)# 绘制聚类和网络模块对应图plotDendroAndColors(dendro = net$dendrograms[[1]], #hclust函数生成的聚类结果colors = mergedColors[net$blockGenes[[1]]],#基因对应的模块颜色groupLabels = "Module colors",#分组标签dendroLabels = FALSE, #false,不显示聚类图的每个分支名称hang = 0.03,#调整聚类图分支所占的高度addGuide = TRUE, #为聚类图添加辅助线guideHang = 0.05,#辅助线所在高度main = "Gene dendrogram and module colors")dev.off()# 加载TOM矩阵load("data/Step02-fpkmTOM-block.1.RData")# 网络特征向量MEs = moduleEigengenes(datExpr, mergedColors)$eigengenes# 对特征向量排序MEs = orderMEs(MEs)# 可视化模块间的相关性sizeGrWindow(5,7.5);pdf(file = "figures/Step02-moduleCor.pdf", width = 5, height = 7.5);par(cex = 0.9)plotEigengeneNetworks(MEs, "", marDendro = c(0,4,1,2), marHeatmap = c(3,4,1,2), cex.lab = 0.8, xLabelsAngle = 90)dev.off()## TOMplotdissTOM = 1-TOMsimilarityFromExpr(datExpr, power = sft$powerEstimate); #1-相关性=相异性nSelect = 400# 随机选取400个基因进行可视化 ， 设置seed值 ， 保证结果的可重复性set.seed(10);#设置随机种子数select = sample(nGenes, size = nSelect);#从5000个基因选择400个selectTOM = dissTOM[select, select];#选择这400*400的矩阵# 对选取的基因进行重新聚类selectTree = hclust(as.dist(selectTOM), method = "average")#用hclust重新聚类selectColors = mergedColors[select];#提取相应的颜色模块# 打开一个绘图窗口sizeGrWindow(9,9)pdf(file = "figures/Step02-TOMplot.pdf", width = 9, height = 9);# 美化图形的设置plotDiss = selectTOM^7;diag(plotDiss) = NA;# 绘制TOM图TOMplot(plotDiss, #拓扑矩阵 ， 该矩阵记录了每个节点之间的相关性selectTree, #基因的聚类结果selectColors, #基因对应的模块颜色main = "Network heatmap plot, selected genes")dev.off()