生物医学大数据解读和分析——构建生物网络实践 WGCNA使用(五)

湖医药Freescience2019-02-19 10:11:43

上期介绍了TCGA的下载,本节开始讲解如何构建共表达网络,首先我们介绍下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)。根据数据本身选择一个合适的阈值参数值,来构建共表达网络,这样更加具有优势。

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


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


2. 安装rstudio,官网 https://www.rstudio.com/



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



4.  运行rstudio。




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



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



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



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



9. 输入数据,打入以下代码到rstudio。

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


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


10. 选择阈值参数,打入以下代码到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()


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



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


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


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



这期的内容讲到如何利用WGCNA构建共表达网络,

下期将继续讲解

cytoscape在构建生物网络中的应用


FS数据挖掘主编


赵忻艺,将大数据应用于医学科研,主要包括临床医学数据的挖掘、收集、整理和利用(标准化和科学化的数据库),医学分子大数据的整理、利用及研究(基因、蛋白及代谢)。特别针对肿瘤个体化的基因测序和数据快速处理,寻找个体化的分子标志物、药物靶标和治疗方案。目前,已建立浙大大数据挖掘团队,旨在降低研究者学习大数据的门槛,推动大数据共享与研究协作,发表更高质量的研究成果,为科研决策提供精准的预测和实验证据。



科研路,不孤单!

Freescience医学科研联盟全国火热招募ing

50家高校及医院的小伙伴已经加入啦

具体点这里


FS科研软件库,集合50+医学科研必备神器,现在统统打包分享点这里




科学自由共享
投稿请扔至freescience@zju.edu.cn

未经许可 不得转载
长按二维码关注


Copyright © 古田计算器虚拟社区@2017