WGCNA 算法

基因组学与应用生物学，2013年，第32卷，第1期，第135-141页

GenomicsandAppliedBiology,2013,Vol.32,No.1,135-141

评述与展望

ReviewandProgress

基于WGCNA算法的基因共表达网络构建理论及其R软件实现

宋长新1*雷萍2王婷1

1青海师范大学计算机学院,青海,810008;2深圳市水质检测中心,深圳,518000

*通讯作者,scx@http://wendang.chazidian.com

摘要WGCNA(weightedgeneco-expressionnetworkanalysis)算法是一种构建基因共表达网络的典型系统生物学算法，该算法基于高通量的基因信使RNA(mRNA)表达芯片数据，被广泛应用于国际生物医学领域。本文旨在介绍WGCNA的基本数理原理，并依托R软件包WGNCA以实例的方式介绍其应用。WGCNA算法首先假定基因网络服从无尺度分布，并定义基因共表达相关矩阵、基因网络形成的邻接函数，然后计算不同节点的相异系数，并据此构建分层聚类树(hierarchicalclusteringtree)，该聚类树的不同分支代表不同的基因模块(module)，模块内基因共表达程度高，而分数不同模块的基因共表达程度低。最后，探索模块与特定表型或疾病的关联关系，最终达到鉴定疾病治疗的靶点基因、基因网络的目的。

关键词WGCNA算法,基因共表达网络,R软件

GeneCo-expressionNetworkAnalysisBasedonWGCNAAlgorithm-TheoryandImplementationinRSoftware

SongChangxin1*LeiPing2WangTing1

1DepartmentofComputerScience,QinghaiNormalUniversity,Qinghai,810008;2ShenzhenWaterQualityCentre,Shenzhen,518000

*Correspondingauthor,scx@http://wendang.chazidian.com

DOI:10.3969/gab.032.000135

AbstractWGCNA(weightedgeneco-expressionnetworkanalysis)isatypicalalgorithmwhichisusedingeneco-expressionnetworkidentification.Thisalgorithmisbasedonhigh-throughoutmRNAgeneexpressionpro-filesandbeingwidelyusedintheinternationalbiomedicalfield.Inthisarticle,wewillintroducethebasicthe-oryandit'simplementationinRsoftware.Firstly,thescale-freeofgenenetworkconditionshouldbesatisfiedbeforeconductingWGCNA,what'smore,itwasnecessarytodefinethecorrelationmatrixofgeneco-expres-sionandadjacencyfunction.Secondly,thedissimilaritymeasurementsofdifferentnodeswerecalculated,andthenhierarchicalclusteringtreewasbuiltbasedonthesedata.Differentdendrogrambranchesrepresentedvari-ousmodules.Thereismuchhigherco-expressionstrengthamonggenesinthesamemodulethanthatindiffer-entmodules.Atlast,itiscriticaltoconnectthemoduleswithinterestingphenotypesordiseaseandidentitythetargetgenesfordiseasetreatment.KeywordsWGCNA,Geneco-expressionnetwork,Rsoftware

基因共表达网络(geneco-expressionnetwork)分

析致力于寻找协同表达的基因模块(module)，并探索

基因网络与研究者关注的表型之间的关联关系。它

基于高通量的微阵列技术(microarray)，应用基因表

达芯片得到实验数据，从转录组(mRNA)水平探索基

因网与疾病或者研究者关注的性状之间的关联关

系，因此应用于复杂疾病的易感基因鉴定、新药开发

等生物医学研究领域。

基金项目：本研究由青海省135高层次人才培养基金资助加权基因共表达网络构建(weightedgeneco-expressionnetwork)算法作为一种高效、准确的生物信息学、生物数据挖掘方法，理论不断完善，应用日渐广泛。首先，在理论研究上，Zhang和Horvath(2005)学者于2005年在发表了第一篇理论分析文章，截至到2012年5月，已有12篇WGCNA相关的文章先后发表在生物信息、系统生物学领域权威期刊上，如BMCBioinformatics、PloSComputationalBiology和Bioin原

136

基因组学与应用生物学

GenomicsandAppliedBiology

formatics等。2008年，

基于该研究方法理论，Langfelder和Horvath(2008)编写并发布了WGCNAR软件包，截至到2012年4月引用次数达到108次之多。

而在实际应用层面，由于基因共表达网络分析作为一种挖掘和呈现基因在不同样本中表达形式的有效方法，可以鉴定高度共表达的基因模块，模块的特征值或模块中包含的关键基因(hupgenes:即与基因网络中的其它基因联系最紧密,在基因网中起关键作用的基因)可用于提炼模块信息，以此深入探讨基因模块或模块中的关键基因和研究人员关注的样本特征间的关联关系。

WGCNA算法已被用于鉴定复杂疾病的候选标记或药物靶点，并用于多项人类复杂疾病的研究。如鉴定家族性混合型高脂血症(Plaisieretal.,2009)、胶质母细胞瘤(Horvathetal.,2006)、自闭症(Voineaguetal.,2011)、阿尔兹海默病(Milleretal.,2010)、骨质疏松症(Farber,2010)的关联基因、生物学通路和肿瘤治疗靶点。

Horvath等(2006)学者将WGCNA方法用于胶质母细胞瘤的研究，充分利用两组共120个患者的胶质母细胞的基因表达数据，挖掘得到一个包含的基因与已知的癌症相关模块“metasignature”中基因高度重叠的基因共表达模块，且异常纺锤型小脑畸形症相关基因(abnormalspindle-likemicrocephalyasso-ciated,ASPM)为该模块中的一个关键基因，ASPM基因已被证实为一个胶质母细胞瘤治疗的靶点基因。

Farber(2010)应用12名低骨密度和14名高骨密度绝经前妇女的单核细胞mRNA基因表达数据，使用WGCNA方法构建表达共基因模块，并发现模块9和骨密度存在显著关联关系，且该结论得到了弗莱明翰(Kieletal.,2007)基于家系和群体数据的基于单核苷酸多态性(SNP)的全基因组关联研究(GWAS)结果和dCOD基因遗传学研究(Styrkarsdottiretal.,2008)的GWAS结果的支持。

综上所述，深入了解WGCNA算法的基本原理，并掌握该方法，将其运用到实际的科学研究中具有极其重要的意义。本文首次在中文期刊上详细介绍WG-CNA算法，旨在在国内推广该系统生物学研究方法。

1原理介绍

WGCNA算法是构建基因共表达网络的常用算法。首先基因共表达网络的概念为：每个节点代表一个基因，在不同样本中存在表达共性的基因处于同一个基因网络，而基因间的共表达关系一般由它们之间的表达相关系数衡量。一般而言，同一基因共表

达网络中的基因表达形式相似、而不同基因共表达

网络中的基因表达形式差别较大。构建了基因模块后，再将基因模块与一些表型信息，例如是否患病、身高和体重等联系起来，以探索基因模块形成的根本原因，鉴定出与表型相关的基因模块。以下具体介绍WGCAN算法(陈超,2011)。1.1网络构建前提

在WGCNA算法中，定义的基因共表达矩阵中的元素为基因的相关系数的加权值，权重的选择标准为使得每个基因网中包含的基因之间的连接服从无尺

度网络分布(scale-freenetworks)(Barab佗si,2009)，

即连接数为i的概率p(i)与i的n次方成反比，即p(i)~i-r。在实际应用中研究者通过选择加权系数来逼近无尺度网络分布，使之满足如下条件：连接节点个数的对数(log(i))与此节点出现概率的对数值(log(p(i)))为负相关关系(注意:相关系数至少应达到0.8)，与此同时，不同模块中的基因的平均连接度应比较高。1.2网络构建步骤

1.2.1基因共表达相关矩阵的定义

基因共表达相关矩阵中的元素为基因之间的两两相关系数，即每对基因m和基因n的相关系数为Smn=|cor(m,n)|，据此构成了基因共表达相关矩阵S=[Smn]。1.2.2定义邻接函数

最直接的邻接函数通过指定基因之间的相关系数的阈值(如R=0.85)将基因对划分为相关的和不相关的，这种分法虽然简单易行，但将损失大量信息，如将阈值定为0.85时，即便是相关系数为0.84的基因对也将被划分到“不相关”的组中。鉴于此，WGCNA算法中应用了幂指数邻接函数。即对于任何基因对，用邻接系数a指数加权。

Smn作为衡量它们之间相关关系的指标：

amn=power(mn,茁)=|Smn|茁，

即对相关系数进行次方的幂1.2.3确定邻接函数的参数

根据无尺度网络原则确定加权系数茁，即：连接节点个数取对数(log(k))和节点出现的概率的对数值(log(p(k)))之间的相关系数至少达到0.8。1.2.4节点间的相异度衡量

在确定了邻接函数参数茁后，接下来将相关矩阵S=[Smn]转换成邻接矩阵A=[amn]。将某个基因和分析中其他所有基因之间的关系纳入考虑，邻接矩阵被转换成拓扑矩阵赘=[棕mn](Ravaszetal.,2002)，矩阵中的元素如下所示：

基于WGCNA算法的基因共表达网络构建理论及其R软件实现

GeneCo-expressionNetworkAnalysisBasedonWGCNAAlgorithm-TheoryandImplementationinRSoftware

137

棕lmnmn=

km,mn

min{kn}+1-amn

公式中lmn=移滋am滋a滋n代表和基因m、n都存在连接的节点的邻接系数乘积和，km=移滋am滋为基因m单独连接的节点的邻接系数加和，类似地，kn表基因n单独连接的节点的邻接系数加和。在基因=移滋an滋代m和基因n之间无连接，且无任何其它的基因将这

两个基因连接的情况下棕mn异程度用dmn棕

=1-棕mn来衡量，=0d。与此同时，

mn棕是网络构建的基础。节点的相1.2.5聚类分析鉴定基因模块

以基因之间的相异系数dij棕为基本元素构建分层聚类树(hierarchicalclusteringtree)，聚类数的不同分支代表不同的基因模块。构建聚类树有静态剪切树和动态剪切树两种算法(Langfelderetal.,2008)。简而言之，静态剪切树是通过定义一个固定的高度将群集分支，该方法识别群集的准确度不高。而动态剪切算法基于树状图的分支形状，它可挖掘得到静态剪切无法检测出的基因模块，更重要的是，用动态剪切算法鉴定出的基因网络以往的生物实验结果比较一致(DongandHorvath,2007)。两种动态剪切树算法简介如下：(1)动态树法，为一种“自上而下”的算法，由不断的分解和组合这一反复迭代过程，在迭代结果稳定后，鉴定出基因模块，该算法由根据静态剪切法得到的几个较大的基因模块开始，对每个群集加入波动特征模式进行精细的分析和分割；(2)动态混合切割法，是一种“自下而上”的算法，主要由识别和测试步骤组成，首先识别满足如下条件的基因模块：(a)模块

中基因个数满足设定的最低数目；

(b)从固定的基因模块中移除距离过远的分支；(c)基因模块间的区别明显；(d)基因模块中的关键基因(hupgene)应紧紧连接；测试阶段为测试未分配的基因，将其分配进足够接近的初步模块，并最终形成基因网络/模块结构。

基因模块与已知特征相联系的具体办法有：(1)计算得到基因网络/模块的特征值，再计算模块的特征向量与关注表型的相关系数；(2)对于分组表型数据(如疾病状态等)，可以首先定义应用T检验计算每个基因在不同组(如疾病组和正常组)间的基因mRNA差异表达的显著性检验P值，并将显著性P值的以10为底的对数值定义为基因显著性(GS:genesignificance)，再将每一个模块显著性(modulesignificance,MS)定义为模块中所包含基因的GS的平均值。然后比较MS的值，一般而言，某模块的MS值显著高于其它模块说明这一个模块可能与疾病存在关联关系；(3)利用基因网络中的包含和其它基因连接度较高的关键基因的信息来推测该基因网络/模块的成因。

2R软件实例

本小节将详细介绍在R软件应用中(R2.3.12,http://wendang.chazidian.com/)，应用WGCNA软件包(LangfelderandHorvath,2008)的实例。具体为生成模拟基因表达谱数据以及表型数据，构建基因模块，并探索其与表型之间的关联关系。有关WGCNA的最新研究进展参见其主页http://www.genetics.ucla.edu/labs/horvath/CoexpressionNetwork/。2.1模拟数据

首先在WGCNA算法主页下载WGCNAPack-age并将其安装在PC上，

设置好程序运行的工作空间，并加载程序包，然后模拟基因mRNA表达谱以及样本表型数据，得到30个样本的3600个基因的表达数据和样本的患病状态，用变量y表示(y=0代表正常组,y=1代表患病组)。具体的R程序如下所示(符号“#”后为注释)：###模拟数据

setwd("E:/WGCNA/simulate")#设定工作空间,请根据实际情况设置

library("WGCNA")#加载程序包options(stringsAsFactors=FALSE)#预设置#步骤一：预设参数no.obs<-30#设置样本个数ESturquoise<-0;ESbrown<--.6;ESgreen<-.6;ESyellow<-0#设置基因模块的显著性水平

ESvector<-c(ESturquoise,ESbrown,ESgreen,ESyellow)

nGenes1<-3000#设置基因个数为3000simulateProportions1<-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^

138

基因组学与应用生物学

GenomicsandAppliedBiology

2)*rnorm(no.obs)

MEN1<-data.frame(y,MEturquoise,MEblue,MEbrown,MEgreen,MEyellow)

dat1<-simulateDatExpr5Modules(MEturquoise=MEN1$MEturquoise,

MEblue=MEN1$MEblue,MEbrown=MEN1$ME-brown,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]]=ArrayNameDimnames(datExpr)[[2]]=GeneName

为考察数据中是否包含有离群的样本，我们提交如下程序绘制样本聚类图，并且将各样本的患病状态标注在图中：

plotClusterTreeSamples(datExpr=datExpr,y=y)#绘制样本聚类图。

观察样本聚类图(图1)发现模拟数据中并未包含有离群样本。图中下半部分中样本对应为红色的表示其属于患病组，而黑色部分代表正常组。2.2网络构建

2.2.1选择构建网络参数

为了尽量满足无尺度网络分布前提条件，编写如下R程序探索邻接矩阵权重参数茁的取值：

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

powers=c(c(1:10),seq(from=12,to=20,by=2))sft=pickSoftThreshold(datExpr,powerVector=powers,verbose=5)#图形绘制

parb(mfrow=c(1,2));cex1=0.9;

plotb(sft$fitIndices[,1],-signb(sft$fitIndices[,3])*sft$fitIndices[,2],

xlab="SoftThreshold(power)",ylab="ScaleFreeTopologyModelFit,signedR^2",type="n",

main=paste("Scaleindependence"),ylim=c(0,1));text(sft$fitIndices[,1],-sign(sft$fitIndices[,3])*sft

图1系统聚类树的结果与对应的样本信息

Figure1Thephylogenetictreeandtheinformationofthe

sample

$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("Meanconnectivity"),

ylab="MeanConnectivity",xlab="SoftThreshold(power)")

text(sft$fitIndices[,1],sft$fitIndices[,5],labels=powers,cex=cex1,col="black")

使参数茁取值从1~18，计算相应的模型选择统计量绘制图形(图2)，这两个分图的横轴均代表权重参数茁相关系数的平方。相关系数的平方取值越高，

左图的纵轴代表对应的网络中log(k)与log(最小应(p(k))达到0.8)，说明该网络越逼近无网络尺度的分布。B图的纵轴代表对应的基因模块中所有基因邻接函数的均值。本研究发现在模拟的数据中，选择log(k)与log(p(k))的相关系数的平方取值均比较高，因此我们选择其值

基于WGCNA算法的基因共表达网络构建理论及其R软件实现

GeneCo-expressionNetworkAnalysisBasedonWGCNAAlgorithm-TheoryandImplementationinRSoftware

139

图2茁参数选择图

注:A:刻度独立值;B:平均连通性

Figure2Figuresforselect茁

Note:A:Scaleindependence;B:Mean

connectivity首次接近0.90时茁的取值，即茁=2，据此构建基因网络。2.2.2构建网络

首先计算基因间的相异性系数，并得到系统聚类

树，然后按照混合动态剪切树的标准，并设置每个基因网络/模块最少的基因数目为30。在使用动态剪切法在确定基因模块后，我们依次计算每个模块的特征

向量值(eigengenes)，

然后对模块进行聚类分析，将距离较近的模块合并成新的模块(可设置height=0.2)。本

研究最终生成的模块聚类图如图3所示。图中不同颜色代表不同的基因模块，而灰色部分代表无法合并到任何其它模块中的基因。本步骤的R程序如下所示：softPower<-2;#设定无尺度参数为2

adjacency=adjacency(datExpr,power=softPower);TOM=TOMsimilarity(adjacency);dissTOM=1-TOM

#聚类分析

geneTree=flashClust(as.dist(dissTOM),method=

"average");

minModuleSize=30;

#设置基因网中至少包含30个基因

dynamicMods=cutreeDynamic(dendro=geneTree,

图3基因的系统聚类数和动态剪切法生成的基因网络/模块Figure3Clusterdendrogrambasedondynamictree

cut

distM=dissTOM,

deepSplit=2,pamRespectsDendro=FALSE,minClusterSize=minModuleSize);

dynamicColors=labels2colors(dynamicMods)table(dynamicColors)#计算基因网的特征值

MEList=moduleEigengenes(datExpr,colors=dy-namicColors)

MEs=MEList$eigengenesMEDiss=1-cor(MEs);METree=flashClust(as.dist(MEDiss),method="av

erage");

#将特征值距离小于0.2的基因网合并

MEDissThres=0.2

merge=mergeCloseModules(datExpr,dynamicCol-ors,cutHeight=MEDissThres,verbose=

3)

mergedColors=merge$colors;mergedMEs=merge$newMEs;

##重新命名合并后的基因网

moduleColors=mergedColorscolorOrder=c("grey",standardColors(50));

moduleLabels=match(moduleColors,colorOrder)-1;MEs=mergedMEs;##绘制最终基因网络构建图

plotDendroAndColors(geneTree,mergedColors,

"Mergeddynamic",dendroLabels=FALSE,hang=0.03,addGuide=TRUE,guideHang=0.05)

2.3网络与表型关联分析

2.3.1相关系数法计算每个基因模块的特征值与患病状态(变量y)

的相关系数，编写程序如下所示：

Signif(cor(y,datME,use="p"),2)p.values=corPvalueStudent(cor(y,datME,use="p"),nSamples=length(y))得到结果见表1。