构建生物网络实践 WGCNA使用

“

生信草堂

将会与更多的优秀微信公众号合作，把最优秀的微信推文呈现给大家，希望可以帮助读者更多的了解生信技术，培养和提高读者的生信分析能力！

号外，号外，号外

你想和生信分析大神做好朋友么？

你想认识更多爱好生信分析的小伙伴么？

你想让自己的生信分析走上快车道么？

那就赶快加入我们的生信交流微信群吧！

正确加入我们的模式是：

添加我们的微信Edison686868或者mly-1800为好友

标注“加入生信草堂交流群”

在群里请大家注明自己本名，单位，研究领域

便于小编管理

Freescience联盟 

由高校、医院FS公众号和科研技能公众号等百家单位联合创建的科研交流分享平台；联盟的宗旨：“公正至上，自由分享，平等共赢”。欢迎您的关注，让我们共同学习进步。戳这里Freescience联盟公众号原文，请多关注哦~

本期开始讲解如何构建共表达网络，首先我们介绍下WGCNA（Weighted Correlation Network analysis），它是一种从芯片数据中挖掘模块(module)信息的算法。

如果在一个生理过程或不同组织中某些基因总是具有相类似的表达变化，那么我们有理由认为这些表达类似的基因在功能上是相关的，可以把它们定义为一个模块(module)。

当然除了寻找有生物功能的模块外，我们还可以通过WGCNA实现性状关联，代谢通路建模，基因共表达网络构建等分析。

这里主要讲解如何构建共表达网络（co-expression network）及其生物学意义。首先，我们了解下什么是共表达网络。

在共表达网络中，每一个基因在一个特定时间或空间的表达情况被视做一个点(node)。如果我们做了100张芯片，每张芯片上有2万个基因，那么我们可以用一个100*20000的矩阵文本来表示实验结果。

为了得到基因间的关联情况，我们需要计算任何两个基因间的相关系数(一般采用Person Coefficient)，经过运算后，我们得到一个20000*20000的关系矩阵S，sij表示第i个基因和第j个基因的Person Coefficient，即两个基因的表达相关性。 

为了知道两个基因的表达是否具有相关性，需要人为规定一个阈值，只有当基因间的Person Coefficient达到这一阈值后(如+0.7，正为正相关，负为负相关)。

我们才认为这两个基因是相关，否则不相关。对所有基因进行两两比较运算，这样就又得到一个20000*20000的邻接矩阵（由0或1构成，1表示相关，0表示不相关）。通过这个矩阵就可以构成网络。

这是一般构建共表达网络的方法，实现以上流程，只需调用R语言内置的cor()函数即可。

但是以上分析方法存在一个很明显的局限性，即我们没有理由认为Person Coefficient为0.7的两个基因与Person Coefficient为0.69的两个基因是有显著差别的。

而WGCNA采用了一种基于软阈值的判定方法，这很好地避免了这一问题。

软阈值的思想是通过权函数将Adjacency Matrix中的元素连续化(所以才称之为Weighted Network)。根据数据本身选择一个合适的阈值参数值，来构建共表达网络，这样更加具有优势。

接下来将对这个流程进行还原。以下软件工具都是开源免费。

01

安装R语言。官网  https://www.r-project.org/。

02

安装rstudio，官网 https://www.rstudio.com/。

03

进入WGCNA官网。
https://labs.genetics.ucla.edu/horvath/CoexpressionNetwork/Rpackages/WGCNA/

04

运行rstudio。

05

根据官网WGCNA的安装代码复制黏贴到rstudio中。注意提示的版本要求等。

06

意提示的版本要求等。

07

运行WGCNA程序。打入“library(WGCNA)”代码到rstudio。如图即安装成功。

08

准备数据，这里以还原文章中的TCGA的乳腺癌的三级数据为例，行是1,218 个样本，列是20,531个 基因的矩阵格式。(表达值由log2(rsem+1)转化过)。

09

输入数据，打入以下代码到rstudio。

红色指根据个人情况自己设置。

setwd ("数据路径") 
enableWGCNAThreads(nThreads = NULL)
rt=read.table("数据名称",sep="\t",row.names=1,header=T,check.names=F,quote="!")
datExpr = t(rt)

再输入以下代码

powers1=c(seq(1,10,by=1),seq(12,30,by=2))

RpowerTable=pickSoftThreshold(datExpr, powerVector=powers1)[[2]]

09

选择阈值参数，打入以下代码到rstudio。最终得到softThresholding.pdf文件。

cex1=0.7
pdf(file="softThresholding.pdf")
par(mfrow=c(1,2))
plot(RpowerTable[,1], -sign(RpowerTable[,3])*RpowerTable[,2],xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n")
text(RpowerTable[,1], -sign(RpowerTable[,3])*RpowerTable[,2], labels=powers1,cex=cex1,col="red")
abline(h=0.85,col="red")
plot(RpowerTable[,1], RpowerTable[,5],xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n")
text(RpowerTable[,1], RpowerTable[,5], labels=powers1, cex=cex1,col="red")
dev.off()

011

打开softThresholding.pdf。第一张图红色线相交的值为选择阈值，例如选择beta1=10。

012

根据选择的阈值生成邻接矩阵。打入以下代码到rstudio。注意这步非常耗内存，以上数据运行已超过20g内存的使用量，所以一般的台式机是无法运行，需在服务器上运行。当然可以通过减少基因数来运行。

beta1=10
ADJ= adjacency(datExpr,power=beta1)
vis=exportNetworkToCytoscape(ADJ,edgeFile="edge.txt",nodeFile="node.txt",threshold = 0.7)

013

最终我们发现目录下生成了两个文件，一个是edge.txt和node.txt的信息。将edge.txt导入到cytoscape即可构建成网络图了。

学术手拉手

长按关注生信草堂

长按关注Freescience联盟