wgcna_1

knitr::opts_chunk$set(echo = TRUE)

背景

1.无尺度网络

相比于只关注差异表达的基因，WGCNA利用数千或近万个变化最大的基因或全部基因的信息识别感兴趣的基因集，并与表型进行显著性关联分析。既充分利用了信息，也把数千个基因与表型的关联转换为数个基因集与表型的关联，免去了多重假设检验校正的问题。

存在少数节点具有明显高于一般点的度，这些点被称为hub。由少数hub与其它节点关联，最终构成整个网络。这样的网络的节点度数与具有该度数的节点个数间服从power distribution。

2.相关术语

• 共表达网络：点代表基因，边代表基因表达相关性。加权是指对相关性值进行冥次运算 (冥次的值也就是软阈值 (power, pickSoftThreshold这个函数所做的就是确定合适的power))。无向网络(unsigned network)的边属性计算方式为 abs(cor(genex, geney)) ^ power；有向网络(signed network)的边属性计算方式为 (1+cor(genex, geney)/2) ^ power; sign hybrid的边属性计算方式为cor(genex, geney)^power if cor>0 else 0， sign hybrid意味着它既包含加权网络也包含非加权网络。这种处理方式强化了强相关，弱化了弱相关或负相关，使得相关性数值更符合无标度网络特征，更具有生物意义。除了软阈值还有硬阈值一说，计算方式是 a_ij = 1 if s_ij > β otherwise a_ij = 0。这里的β就是硬阈值(hard threshold)。

• Module(模块)：高度內连的基因集。在无向网络中，模块内是高度相关的基因。在有向网络中，模块内是高度正相关的基因。

• Connectivity (连接度)：类似于网络中 “度” (degree)的概念。每个基因的连接度是与其相连的基因的边属性之和。节点的重要程度

• Module eigengene E: 给定模型的第一主成分，代表整个模型的基因表达谱。即用一个向量代替了一个矩阵，方便后期计算。

• Intramodular connectivity: 给定基因与给定模型内其他基因的关联度，判断基因所属关系。

• Adjacency matrix (邻接矩阵)：基因和基因之间的加权相关性值构成的矩阵。 两个节点间的关系强弱 邻接矩阵

• TOM (Topological overlap matrix)：把邻接矩阵转换为拓扑重叠矩阵，以降低噪音和假相关，获得的新距离矩阵，这个信息可拿来构建网络或绘制TOM图。

3.主要步骤

rm(list = ls())
options(stringsAsFactors = F)

1.读取表达数据

1.1读取数据

exp0 <- read.table("sample_data/LiverFemale_Expression.txt",header = T,row.names = 1)
#过滤掉在所有样品中表达之和小于10 的基因
#exp <- exp0[rowSums(exp0) > 10, ]
exp0[0:4,0:4]

WGCNA需要的格式为每一列是一个基因，每一行是一个样品，因此需要对exp0转置。

1.2转置

datExpr <- t(exp0)
dim(datExpr)
rownames(datExpr)[1:5]

1.3保存数据

write.table(datExpr,file = "datExpr.txt",sep = "\t",row.names = T,quote = F)
#quote表示是否用引号括起来

2.寻找最佳β值

library(WGCNA)
enableWGCNAThreads(nThreads = 20)
#开启多线程模式

#设置power的迭代次数
powers <- c(1:20)
sft <- pickSoftThreshold(datExpr,
                         powerVector = powers,
                         verbose=5#区分正负号，unsigned不区分正负号。
)

sft$powerEstimate

没有返回β值，进行可视化后自己添加。

sizeGrWindow(9,5)
par(mfrow = c(1,2))#子图
cex1 = 0.6

#Scale-free topology fit index as a function of the soft-thresholding power
plot(sft$fitIndices[,1],
     -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=cex1,col="red")

abline(h=0.8,col="red")

#Mean connectivity as a function of the soft-thresholding power
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")

由图可知，R^2值第一次达到0.8时，为最佳β值13。

sft$powerEstimate <- 13

3.构建网络

net = blockwiseModules(
                 datExpr,
                 
                 #计算相关系数
                 corType = "pearson",
                 
                 #计算邻接矩阵
                 power = sft$powerEstimate,
                 networkType = "unsigned",
                 
                 #计算TOM矩阵
                 TOMType = "unsigned",
                 saveTOMs = TRUE,
                 saveTOMFileBase = "blockwiseTOM", # 保存文件前缀
                 
                 #聚类并划分模块
                 deepSplit = 2, # 0|1|2|3|4, 值越大得到的模块就越多越小
                 minModuleSize = 30,#模块内最少包含多少基因
                 
                 #合并相似模块
                 #计算模块特征向量（module eigengenes， MEs），即第一个主成分（1st PC）
                 #计算 MEs 与 datTrait 之间的相关性
                 #对距离小于 mergeCutHeight （1-cor）的模块进行合并
                 mergeCutHeight = 0.15,
                 
                 numericLabels = TRUE, #以数字命名模块
                 nThreads = 0 # 0 则使用所有可用线程

 )

#查看每个模块包含基因数目
table(net$colors) 

4.可视化

#模块名称修改为颜色
moduleColors = labels2colors(net$colors)

plotDendroAndColors(net$dendrograms[[1]], 
                    moduleColors[net$blockGenes[[1]]],
                    "Module colors",
                    dendroLabels = FALSE, 
                    addGuide = TRUE
)

save(datExpr, sft, net, moduleColors, file = "wgcna-network.Rdata")