TCGA数据挖掘(六):WGCNA（加权基因共表达网络分析）
- 2019 年 10 月 4 日
- 笔记
共表达网络的建设从概念上来讲是简单直观的，通过基因表达的相似性可分析基因产物可能的相互作用关系，从而了解基因间相互作用脉络及寻找核心基因。核心基因是重要的枢纽，在网络模块中起关键性作用。
关于共表达分析可阅读这篇文献：
https://www.ncbi.nlm.nih.gov/pubmed/28077403
加权基因共表达网络分析(WGCNA，Weighted gene co-expression network analysis)，WGCNA能够从复杂数据中（N多分组）快速地提取出与样本特征相关的基因共表达模块，以供后续分析。简单地说，它通过计算基因之间的表达相关性，将具有表达相关性的基因聚类到一个模块中，然后再分析模块与样本特征（包括临床特征、手术方式、治疗方法等等）之间的相关性，WGCNA搭建了一座样本特征与基因表达变化之间的桥梁。
理解WGCNA，需要先理解下面几个术语和它们在WGCNA中的定义。
共表达网络：定义为加权基因网络。点代表基因，边代表基因表达相关性。加权是指对相关性值进行冥次运算(冥次的值也就是软阈值 (power,pickSoftThreshold这个函数所做的就是确定合适的power))。无向网络的边属性计算方式为abs(cor(genex, geney)) ^ power；有向网络的边属性计算方式为(1+cor(genex, geney)/2) ^ power; sign hybrid的边属性计算方式为cor(genex, geney)^power if cor>0 else 0。这种处理方式强化了强相关，弱化了弱相关或负相关，使得相关性数值更符合无标度网络特征，更具有生物意义。如果没有合适的power，一般是由于部分样品与其它样品因为某种原因差别太大导致的，可根据具体问题移除部分样品或查看后面的经验值。
Module(模块)：高度内连的基因集。在无向网络中，模块内是高度相关的基因。在有向网络中，模块内是高度正相关的基因。把基因聚类成模块后，可以对每个模块进行三个层次的分析：1. 功能富集分析查看其功能特征是否与研究目的相符；2. 模块与性状进行关联分析，找出与关注性状相关度最高的模块；3. 模块与样本进行关联分析，找到样品特异高表达的模块。
Connectivity (连接度)：类似于网络中 "度"(degree)的概念。每个基因的连接度是与其相连的基因的边属性之和。
Module eigengene E:给定模型的第一主成分，代表整个模型的基因表达谱。
Intramodular connectivity:给定基因与给定模型内其他基因的关联度，判断基因所属关系。
Module membership: 给定基因表达谱与给定模型的eigengene的相关性。
Hub gene: 关键基因 (连接度最多或连接多个模块的基因)。
Adjacency matrix(邻接矩阵)：基因和基因之间的加权相关性值构成的矩阵。
TOM (Topological overlap matrix)：把邻接矩阵转换为拓扑重叠矩阵，以降低噪音和假相关，获得的新距离矩阵，这个信息可拿来构建网络或绘制TOM图。
数据融合
下面我们以之前我们分析的2组数据为例，一个是来自文章：TCGA数据挖掘（四）：表达差异分析（2）或者TCGA数据挖掘（四）：表达差异分析（3）的mRNA的表达差异基因的矩阵，也可以是2者的交集；另外一个是miRNA表达的差异矩阵数据，来自文章：TCGA数据挖掘（五）：miRNA差异分析。
我们将2个文件放在同一个文件夹下
我们首先要做的是将2个矩阵文件合并成一文件，虽然这些数据都是来自于同一数据集，一个是mRNA的数据，一个是miRNA的数据，不同的数据类型中，同一病人的编号不同，也就是TCGA条码，我们在文章：TCGA数据挖掘（一）：TCGAbiolinks包介绍中也有关于TCGA条码的介绍，同一病人的话前12位barcode应该是相同的。所以我们合并数据就是将具有mRNA数据又有miRNA数据的病人合并，只有其中一个的病人数据样本将被抛弃。
setwd("E:\BioInfoCloud\coexpressionMerge") #设置工作目录（需修改） miRNA = read.table("diffmiRNAExp.txt", row.names=1 ,header=T,sep="t",check.names=F) RNA = read.table("diffmRNAExp.txt", row.names=1 ,header=T,sep="t",check.names=F) colnames(miRNA)=gsub("(.*?)\-(.*?)\-(.*?)\-(.*?)\-.*","\1\-\2\-\3\-\4",colnames(miRNA)) colnames(RNA)=gsub("(.*?)\-(.*?)\-(.*?)\-(.*?)\-.*","\1\-\2\-\3\-\4",colnames(RNA)) sameSample=intersect(colnames(miRNA),colnames(RNA)) merge=rbind(id=sameSample,miRNA[,sameSample],RNA[,sameSample]) write.table(merge,file="merge.txt",sep="t",quote=F,col.names=F
运行代码后，会生成一个merge.txt的文件，行名就是miRNA的名称或者是基因名称，列名就是病人的部分barcode
做共表达分析，后面要绘制网络图，需要节点文件，所以这里我们也要生成一个：
miRNALabel=cbind(rownames(miRNA),"miRNA") RNALabel=cbind(rownames(RNA),"gene") nodeLabel=rbind(c("ID","Classify"),miRNALabel,RNALabel) write.table(nodeLabel,file="nodeLabel.txt",sep="t",quote=F,col.names=F,row.names=F)
WGCN分析
install.packages(c("matrixStats", "Hmisc", "splines", "foreach", "doParallel", "reshape", "fastcluster", "dynamicTreeCut", "survival")) install.packages("WGCNA") setwd("E:\BioInfoCloud\WGCNA\2_WGCNA") #设置工作目录（需修改） library(WGCNA) rt=read.table("merge.txt",sep="t",row.names=1,header=T,check.names=F,quote="!") datExpr = t(rt) ######select beta value###### powers1=c(seq(1,10,by=1),seq(12,30,by=2)) RpowerTable=pickSoftThreshold(datExpr, powerVector=powers1)[[2]] 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") # this line corresponds to using an R^2 cut-off of h 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")
运行上面代码后会出现一个图，根据这个图，我们可以选择合适的power值，我们这里选择9.
beta1=9 #计算邻近距离 ADJ= adjacency(datExpr,power=beta1) vis=exportNetworkToCytoscape(ADJ,edgeFile="edge.txt",nodeFile="node.txt",threshold = 0.8)
运行所有代码后我们就会得到一个文件edge.txt的文件。
CytoScape绘制网络图
导入文件
导入node文件
这就构建好网络图了，当然这图没有美化，自己调整网络图了，这就是cytoscape的使用教程啦，大家可以在网上找教程，当然关注本公众号，我们后续也会推出一些软件的使用教程。
