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的使用教程啦，大家可以在網上找教程，當然關注本公眾號，我們後續也會推出一些軟體的使用教程。
