wgcna基因共表达网络构建理论及R软件实现
记录下大佬的文章，方便自己查阅
一、相关概念：
无尺度网络
矩阵乘法
邻接矩阵
邻接矩阵 n次幂 传球问题终极解法
集群系数/邻接系数
皮尔逊积矩相关系数
数学期望 标准差 协方差
基因共表达网络
拓扑矩阵
二、基本原理：
基因共表达网络中，每个节点代表一个基因，不同样本中表达共性的基因处于同一个基因网络，共性关系由表达相关系数衡量。同一个模块内的基因表达形式相似，不同模块内的基因表达形式不相似。构建好表达模块后，将表达模块与患病，身高、体重等表型关联起来
1、构建什么样的网络
构建一个无尺度网络，已有生物学研究表明某些关键基因在一些生理过程中所起的作用为不可替代的，而某些基因在生物学过程中所执行的功能被限制在了某些特定的信号通路中间，以上发现和无尺度网络分布的概念是契合的
2、网络构建步骤
关联矩阵-->确定加权系数-->幂指数加权-->邻接矩阵-->拓扑矩阵-->聚类-->数据检验
s[mn]-->β-->|Smn|β-->A=[amn]-->dωmn
(1)基因共表达矩阵的定义
基因共表达相关矩阵中的元素为基因之间的两两相关系数,即每对基因m和基因n的相关系数为Smn=|cor(m,n)|,此构成了基因共表达相关矩阵S=[Smn]。
(2)定义邻接函数
最直接的邻接函数通过指定基因之间的相关系数的阈值(如R=0.85)将基因对划分为相关的和不相关的,这种分法虽然简单易行,但将损失大量信息,如将阈值定为0.85时,即便是相关系数为0.84的基因对也将被“不相关”的组中。鉴于此,WGCNA算法中应用了幂指数邻接函数。即对于任何基因对,用邻接系数Aij作为衡量它们之间相关关系的指标:Amn=power(Smn,β)=|Smn|β,即对相关系数进行次方的幂指数加权。
(3)确定邻接函数的参数
根据无尺度网络原则确定加权系数β,即:连接节点个数取对数(log(k))和节点出现的概率的对数值(log(p(k)))之间的相关系数至少达到0.8。
(4)节点间的相异度衡量
在确定了邻接函数参数β后,接下来将相关矩阵S=[Smn]转换成邻接矩阵A=[amn]。将某个基因和分析中其他所有基因之间的关系纳入考虑,邻接矩阵被转换成拓扑矩阵Ω=[ωmn](Ravasz et al.,2002),矩阵中的元素如下所示:
公式中lmn=∑μamμaμn代表和基因m、n都存在连接的节点的邻接系数乘积和,km=∑μ Amμ为基因m单独连接的节点的邻接系数加和,类似地,kn=∑μ Anμ代表基因n单独连接的节点的邻接系数加和。在基因m和基因n之间无连接,且无任何其它的基因将这两个基因连接的情况下ωmn=0。与此同时,节点的相异程度用dωmn=1-ωmn来衡量,dωmn是网络构建的基础。
(5) 聚类分析鉴定基因模块
以基因之间的相异系数 d ij ω 为基本元素构建分层聚类树(hierarchical clustering tree),聚类数的不同分支代表不同的基因模块。构建聚类树有静态剪切树和动态剪切树两种算法(Langfelder et al., 2008)。 简而言之,静态剪切树是通过定义一个固定的高度将群集分支,该方法识别群集的准确度不高。而动态剪切算法基于树状图的分支形状,它可挖掘得到静态剪切无法检测出的基因模块,更重要的是,用动态剪切算法鉴定出的基因网络以往的生物实验结果比较一致(Dong and Horvath, 2007)。 两种动态剪切树算法简介如下: (1)动态树法,为一种 “自上而下” 的算法,由不断的分解和组合这一反复迭代过程,在迭代结果稳定后,鉴定出基因模块,该算法由根据静态剪切法得到的几个较大的基因模块开始,对每个群集加入波动特征模式进行精细的分析和分割;
(2)动态混合切割法,是一种 “自下而上” 的算法,主要由识别和测试步骤组成,首先识别满足如下条件的基因模块:
(a)模块中基因个数满足设定的最低数目;
(b)从固定的基因模块中移除距离过远的分支;
(c)基因模块间的区别明显;
(d)基因模块中的关键基因(hup gene)应紧紧连接;测试阶段为测试未分配的基因,将其分配进足够接近的初步模块,并最终形成基因网络 / 模块结构。
基因模块与已知特征相联系的具体办法有:
(1)计算得到基因网络 / 模块的特征值,再计算模块的特征向量与关注表型的相关系数;
(2)对于分组表型数据(如疾病状态等),可以首先定义应用 T 检验计算每个基因在
不同组(如疾病组和正常组)间的基因 mRNA 差异表达的显著性检验 P 值,并将显著性 P 值的以 10 为底的对再将每数值定义为基因显著性(GS: gene significance),一个模块显性(module significance, MS)定义为模块中所包含基因的 GS 的平均值。然后比较 MS 值,
一般而言,某模块的 MS 值显著高于其它模块说明这一个模块可能与疾病存在关联关系; (3)利用基因网络中的包含和其它基因连接度较高的关键基因的信息来推测该基因网络 / 模块的成因
三示例
### 模拟数据
library(fastcluster)
library("WGCNA")
# 加载程序包
options(stringsAsFactors=FALSE) # 预设置
# 步骤一:预设参数
no.obs<-30
# 设置样本个数
ESturquoise<-0; ESbrown<--.6; ESgreen<-.6; ESyellow<-0
# 设置基因模块的显著性水平
ESvector<-c (ESturquoise, ESbrown, ESgreen, ESyellow)
# 设置基因个数为 3 000
nGenes1<-3000
simulateProportions1 <-c (0.2, 0.15, 0.08, 0.06,0.04)
# 模块中基因个数比
set.seed (1)
# 步骤二:模拟基因网
MEgreen<-rnorm (no.obs)
scaledy<-MEgreen*ESgreen+sqrt (1-ESgreen^2)*rnorm (no.obs)
y<-ifelse(scaledy>median (scaledy), 1, 0)
# 模拟的表型变量
MEturquoise <-ESturquoise*scaledy+sqrt (1-ESturquoise^2)*rnorm (no.obs)
MEblue<-.6*MEturquoise+sqrt (1-.6^2)*rnorm (no.obs)
MEbrown<-ESbrown*scaledy+sqrt (1-ESbrown^2 )*rnorm(no.obs)
MEyellow<-ESyellow*scaledy+sqrt (1-ESyellow^2)*rnorm (no.obs)
MEN1<-data.frame (y, MEturquoise, MEblue, MEbrown, MEgreen, MEyellow)
dat1 <-simulateDatExpr5Modules (MEturquoise = MEN1$MEturquoise,
MEblue=MEN1$MEblue,
MEbrown=MEN1$MEbrown,
MEyellow=MEN1$MEyellow,
MEgreen=MEN1$MEgreen,
nGenes=nGenes1,
simulateProportions=simulateProportions1)
datExpr<-dat1$datExpr;
truemodule<-dat1$truemodule;
datME<-dat1$datME;
attach (MEN1)
datExpr<-data.frame (datExpr)
# 模拟的表达谱数据
ArrayName<-paste ("Sample",1:dim (datExpr)[[1]],sep="")
GeneName <-paste ("Gene",1:dim (datExpr) [[2]],sep="")
dimnames (datExpr)[[1]]=ArrayName
dimnames (datExpr)[[2]]=GeneName
plotClusterTreeSamples (datExpr=datExpr, y=y)
# 绘制样本聚类图
powers=c (c (1:10), seq (from=12, to=20, by=2))
sft=pickSoftThreshold(datExpr, powerVector=powers,verbose=5)
# 图形绘制
par(mfrow=c(1,2));
cex1=0.9;
plot(sft$fitIndices[,1],-sign(sft$fitIndices[,3])*sft$fitIndices[,2],
xlab="Soft Threshold (power)", ylab=" Scale FreeTopology Model Fit, signed R^2", type="n",
main=paste ("Scale independence"), ylim=c(0,1));
text (sft$fitIndices[,1], -sign (sft$fitIndices[,3])*sft$fitIndices[,2],
labels=powers,cex=cex1, col="black");
abline (h=0.90, col="black")
plot (sft$fitIndices[,1], sft$fitIndices[,5], type ="n", main=paste ("Mean connectivity"),
ylab =" MeanConnectivity", xlab ="SoftThreshold(power)")
text (sft$fitIndices[,1], sft$fitIndices[,5], labels = powers, cex=cex1, col="black")
softPower<-2;
# 设定无尺度参数为 2
adjacency= adjacency(datExpr, power=softPower);
TOM=TOMsimilarity (adjacency);
dissTOM=1-TOM
# 聚类分析
geneTree=hclust(as.dist(dissTOM), method="average");
minModuleSize=30;
# 设置基因网中至少包含 30 个基因
dynamicMods=cutreeDynamic (dendro=geneTree,
distM=dissTOM,
deepSplit=2, pamRespectsDendro=FALSE,
minClusterSize=minModuleSize);
dynamicColors=labels2colors (dynamicMods)
table (dynamicColors)
# 计算基因网的特征值
MEList=moduleEigengenes(datExpr, colors=dynamicColors)
MEs=MEList$eigengenes
MEDiss=1-cor (MEs);
METree=hclust(as.dist(MEDiss), method="average");
# 将特征值距离小于 0.2 的基因网合并
MEDissThres=0.2
merge=mergeCloseModules(datExpr,dynamicColors,cutHeight=MEDissThres,verbose=3)
mergedColors=merge$colors;
mergedMEs=merge$newMEs;
## 重新命名合并后的基因网
moduleColors=mergedColors
colorOrder=c ("grey", standardColors (50));
moduleLabels=match(moduleColors,colorOrder)-1;
MEs=mergedMEs;
## 绘制最终基因网络构建图
plotDendroAndColors (geneTree, mergedColors,
"Merged dynamic",
dendroLabels=FALSE, hang=0.03,
addGuide=TRUE, guideHang=0.05)
#计算每个基因模块的特征值与患病状态(变量 y)的相关系数，编写程序如下所示
signif(cor(y, datME, use="p"), 2)
p.values=corPvalueStudent(cor (y, datME, use="p"),nSamples=length(y))
# 计算基因的显著性
GS1=as.numeric (cor (y, datExpr, use="p"))
GeneSignificance=abs (GS1)
# 计算基因模块的显著性
ModuleSignificance=tapply (GeneSignificance, moduleColors, mean, na.rm=T)
# 图形绘制
plotModuleSignificance(GeneSignificance, moduleColors, ylim=c(0,0.5), main="Module Significance")
