cdna文库与基因组文库特点（实践演示WGCNA如何构建权重基因共表达网络）

WGCNA（weighted gene co-expression network analysis）算法是一种构建基因共表达网络的典型系统生物学算法，可以系统地反馈我们各样本基因之间的相互作用模式。

大致步骤：首先假定基因网络服从无尺度分布，并定义基因共表达相关矩阵、基因网络形成的邻接函数，然后计算不同节点的相异系数，并据此构建分层聚类树（hierarchical clustering tree），该聚类树的不同分支代表不同的基因模块，模块内基因共表达程度高，而分数不同模块的基因共表达程度低。

WGCNA是用来分析基因表达数据的，能够将基因表达与表型变化关联起来，从而挖掘影响性状变化的模块（核心基因），理论上这个模块中的基因都是发挥某一个生物学功能，具有相同的表达模式。通过对模块的进一步的分析，能够实现筛选模块的核心基因，关联性状，代谢通路建模，建立基因互作网络等高级分析。

WGCNA是基于R语言编写的一款软件包，想要了解WGCNA的研究进展请参考主页（https://wenku.baidu.com/view/faea2d2dad02de80d5d8401b.html）。网络上关于WGCNA的分析有很多，但大都局限在理论上面，鲜有实践。今天，就给大家做个简单的演示，告诉大家数据是具体怎么处理的！

#首先载入WGCNA包，设置随机种子，默认数据不进行因子转换

library(WGCNA)

set.seed(1)

options(stringsAsFactors = F)

#构造性状数据（亦或是分组数据）

obs<-paste("sam",1:10,sep="")

sample<-as.data.frame(diag(x=1,nrow = length(obs)))

rownames(sample)<-obs

colnames(sample)<-1:10

#构造表达量数据

datExpr<-as.data.frame(t(matrix(runif(30000) 5,3000,10)))

rownames(datExpr)<-obs

names(datExpr)<-paste("transcript",1:3000,sep="")

#明确样本数和基因数

nGenes = ncol(datExpr)

nSamples = nrow(datExpr)

#首先针对样本做个系统聚类树

datExpr_tree<-hclust(dist(datExpr), method = "average")

par(mar = c(0,5,2,0))

plot(datExpr_tree, main = "Sample clustering", sub="", xlab="", cex.lab = 2,

cex.axis = 1, cex.main = 1,cex.lab=1)

#针对前面构造的样品矩阵添加对应颜色

sample_colors <- numbers2colors(sample, colors = c("white","blue"),signed = FALSE)

#构造10个样品的系统聚类树及性状热图

par(mar = c(1,4,3,1),cex=0.8)

plotDendroAndColors(datExpr_tree, sample_colors,

groupLabels = colnames(sample),

cex.dendroLabels = 0.8,

marAll = c(1, 4, 3, 1),

cex.rowText = 0.01,

main = "Sample dendrogram and trait heatmap")

#选择合适“软阀值（soft thresholding power）”beta

powers = c(1:30,mar=c(4,4,4,4))

pow<-pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)

#设置网络构建参数选择范围，计算无尺度分布拓扑矩阵

par(mfrow = c(1,2))

plot(pow$fitIndices[,1], -sign(pow$fitIndices[,3])*pow$fitIndices[,2],xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n", main = paste("Scale independence"))

text(pow$fitIndices[,1], -sign(pow$fitIndices[,3])*pow$fitIndices[,2],labels=powers,cex=0.5,col="red")

plot(pow$fitIndices[,1], pow$fitIndices[,5],xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",main = paste("Mean connectivity"))

text(pow$fitIndices[,1], pow$fitIndices[,5], labels=powers, cex=0.6,col="red")

参数beta取值默认是1：30，上述图形的横轴均代表权重参数β，左图纵轴代表对应的网络中log(k)与log(p(k))相关系数的平方。相关系数的平方越高，说明该网络越逼近无网路尺度的分布。右图的纵轴代表对应的基因模块中所有基因邻接函数的均值。

接下来是非常重要的一块内容就是构架基因网络。

大体思路：计算基因间的邻接性，根据邻接性计算基因间的相似性，然后推出基因间的相异性系数，并据此得到基因间的系统聚类树。然后按照混合动态剪切树的标准，设置每个基因模块最少的基因数目为30。根据动态剪切法确定基因模块后，再次分析，依次计算每个模块的特征向量值，然后对模块进行聚类分析，将距离较近的模块合并为新的模块。

softPower = 6

#1、计算树之间的邻接性

adjacency = adjacency(datExpr, power = softPower)

#2、计算树之间的相异性

TOM = TOMsimilarity(adjacency)

dissTOM = 1-TOM

#3、聚类分析

geneTree = hclust(as.dist(dissTOM), method = "average")

#4、设置基因模块中最少基因包含30个

minModuleSize = 30

dynamicMods = cutreeDynamic(dendro = geneTree, distM = dissTOM, deepSplit = 2,pamRespectsDendro = FALSE, minClusterSize = minModuleSize)

#5、基因分组上色

dynamicColors = labels2colors(dynamicMods)

table(dynamicColors)

#6、计算基因模块的特征值

MEList = moduleEigengenes(datExpr, colors = dynamicColors)

MEs = MEList$eigengenes

MEDiss = 1-cor(MEs)

METree = hclust(as.dist(MEDiss), method = "average")

#7、建立系统聚类树

MEDissThres = 0.4

#8、基因模块合并

merge = mergeCloseModules(datExpr, dynamicColors, cutHeight = MEDissThres, verbose = 3)

mergedColors = merge$colors

mergedMEs = merge$newMEs

#9、绘制基因网络图

plotDendroAndColors(geneTree, cbind(dynamicColors, mergedColors), c("Dynamic Tree Cut", "Merged dynamic"), dendroLabels = FALSE, hang = 0.03, addGuide = TRUE, guideHang = 0.05)

plot(METree, main = "Clustering of module eigengenes", xlab = "", sub = "")

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

好了，今天的演示就到这里。想要详细了解WGCNA的分析，请参考https://labs.genetics.ucla.edu/horvath/CoexpressionNetwork/Rpackages/WGCNA/

,