示例數據、R代碼的百度盤連接（提取碼，zb33）：node
https://pan.baidu.com/s/1UE8c5ev-gCA-48_6vCcUKQweb
WGCNA的加權原理算法
在下文介紹具體過程以前，首先來理解一下WGCNA的加權原理。數據庫
傳統方法上，描述兩個基因間的關聯程度，可經過計算表達值間的Pearson、Spearman等相關係數得到。爲了構建關聯網絡，一般指定一個篩選閾值，如相關係數大於0.8以上，做爲兩個基因間具備強關聯程度的依據。隨後在網絡中，節點表示基因，若兩個基因間的相關係數大於指定的閾值，則會以邊鏈接以表示二者間可能存在相互做用。express
可是基於固定閾值法的缺點在於，閾值是人爲定義的，將會忽略不少潛在關聯。例如，0.79就是不相關嗎？同時，這種一刀切的方法也會丟失基因的變化趨勢信息，將難以在網絡中描述相關性的強弱關係。微信
爲了解決這些問題，提出了「加權」的思想。WGCNA的作法是對基因表達值之間的相關係數取β次冪，直接結果是把基因間相關性的強弱的差別放大。這樣作的好處是使強弱關係更爲分明，有利於後續聚類（模塊）識別。網絡
以下所示，對於基因i和j，相關係數爲rij，取β次冪後得到aij，最終表徵基因的相關強度。app
所以，最終網絡中基因間的相關強度，取決於β次冪的選擇，其取值的高低直接影響模塊的構建和模塊內基因的劃分。β值根據接近無標度網絡（scale-free network）的最低值肯定，下文的計算過程當中，power值的選擇即爲肯定β值的過程。編輯器
根據這個思想，計算全部基因間表達水平的關聯強度，得到鄰接矩陣。並進一步，計算拓撲重疊矩陣（TOM），簡單來講若是基因i和j有不少相同的鄰接基因，那麼TOMi,j就高，意味着兩基因間有類似的表達模式。所以，TOM矩陣用於反映基因間表達的類似性。ide
在這些過程當中，再也不根據某關聯強度閾值或者類似度閾值作篩選，所以不會丟失基因的變化趨勢信息。過程當中綜合考慮全部給定基因間的表達關係，根據基因間表達的類似性，進行層次聚類劃分功能模塊，進而識別模塊特徵基因、肯定功能模塊中基因表達和性狀的關聯等，詳見下文過程。
輸入數據示例
網盤附件「gene_express.txt」爲某項關於腫瘤的研究中經過轉錄組測序得到的患者腫瘤組織的基因表達值矩陣，包含了24個樣本中基因表達的FPKM值信息。
「trait.txt」則記錄了24個腫瘤組織樣本所來源的患者的臨牀資料信息，包括性別、年齡、是否有吸菸史、飲酒史以及腫瘤分期、生存狀態等信息。
接下來，經過WGCNA識別腫瘤組織中基因的共表達概況，並指望尋找與疾病進展相關的基因集。
加權網絡構建及共表達模塊劃分
首先加載R包，並讀取基因表達值矩陣。
能夠進行一些前處理，比方說處理缺失值，按表達水平過濾低表達的基因等，避免它們的影響，同時減小後續運算資源消耗。
#WGCNA 包經過 bioconductor 安裝
#install.packages('BiocManager')
#BiocManager::install('WGCNA')
library(WGCNA)
#示例數據，基因表達值矩陣
gene <- read.delim('gene_express.txt', row.names = 1, check.names = FALSE)
#例如過濾低表達值的基因，只保留平均表達值在 1 以上的
gene <- subset(gene, rowSums(gene)/ncol(gene) >= 1)
#轉置矩陣，行是基因，列是樣本
gene <- t(gene)
無標度拓撲分析肯定最佳β值
爲了構建加權共表達網絡，首先須要肯定一個軟閾powers值（soft thresholding powers），做爲上文提到的相關係數的β次冪，以創建鄰接矩陣並根據連通度使基因分佈符合無尺度網絡。
pickSoftThreshold()可經過將給定的基因表達矩陣按表達類似度計算加權網絡，並分析指定powers值（構建一組向量，傳遞給powerVector）下網絡的無標度拓撲擬合指數以及連通度等信息，目的是幫助爲網絡構建選擇合適的powers值。
##power 值選擇
powers <- 1:20
sft <- pickSoftThreshold(gene, powerVector = powers, verbose = 5)
#擬合指數與 power 值散點圖
par(mfrow = c(1, 2))
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2], type = 'n',
xlab = 'Soft Threshold (power)', ylab = 'Scale Free Topology Model Fit,signed R^2',
main = paste('Scale independence'))
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
labels = powers, col = 'red');
abline(h = 0.90, col = 'red')
#平均連通性與 power 值散點圖
plot(sft$fitIndices[,1], sft$fitIndices[,5],
xlab = 'Soft Threshold (power)', ylab = 'Mean Connectivity', type = 'n',
main = paste('Mean connectivity'))
text(sft$fitIndices[,1], sft$fitIndices[,5], labels = powers, col = 'red')
左圖中縱軸展現了無標度拓撲擬合指數R2（scale free topology fitting index R2，即統計結果中的SFT.R.sq列數值），之因此有負值是由於又乘以了slope列數值的負方向，僅關注正值便可；右圖縱軸展現了網絡的平均連通度。
綜合考慮這些指標選擇合適的軟閾powers值，使R2可能大但也要保證連通度不能過低。看到網上有教程中提到，通常選擇R2大於0.85時的power值，那麼在該示例中則需選擇powers=8。
pickSoftThreshold()自身也估計了一個最佳的參考值，能夠看到一樣爲powers=8。
得到基因共表達類似度矩陣
接下來，就使用上述估計的最合適的power值，構建無標度網絡。
首先基於基因表達值矩陣計算基因間表達的類似度，得到鄰接度矩陣，並進一步計算拓撲重疊矩陣（topological overlap matrix，TOM）。
#上一步估計的最佳 power 值
powers <- sft$powerEstimate
#得到 TOM 矩陣
adjacency <- adjacency(gene, power = powers)
tom_sim <- TOMsimilarity(adjacency)
rownames(tom_sim) <- rownames(adjacency)
colnames(tom_sim) <- colnames(adjacency)
tom_sim[1:6,1:6]
#輸出 TOM 矩陣
#write.table(tom_sim, 'TOMsimilarity.txt', sep = '\t', col.names = NA, quote = FALSE)
TOM矩陣中的值就反映了兩兩基因間共表達的類似度，取值0-1，基因間共表達類似度越高，值越趨於1。
共表達模塊識別
上述得到了基因間的共表達類似度矩陣，記錄了基因間表達的類似性。能夠進一步據此計算距離並繪製聚類樹，查看類似度的總體分佈特徵。
#TOM 相異度 = 1 – TOM 類似度
tom_dis <- 1 - tom_sim
#層次聚類樹，使用中值的非權重成對組法的平均聚合聚類
geneTree <- hclust(as.dist(tom_dis), method = 'average')
plot(geneTree, xlab = '', sub = '', main = 'Gene clustering on TOM-based dissimilarity',
labels = FALSE, hang = 0.04)
可以明顯地看到，這些基因根據共表達類似性模式，能夠劃分爲不一樣的聚類簇。每一簇中的基因間具備較高的共表達類似性，簇間的基因共表達類似性較低。
所以，就能夠繼續根據這種共表達模式（聚類簇特徵），劃分出幾個更爲明顯的共表達模塊。不難想到，同一模塊內的基因具備較高的共表達類似性，即存在較高的協同，暗示它們參與了類似的調節途徑，或者在類似的細胞區域中發揮功能。
#使用動態剪切樹挖掘模塊
minModuleSize <- 30 #模塊基因數目
dynamicMods <- cutreeDynamic(dendro = geneTree, distM = tom_dis,
deepSplit = 2, pamRespectsDendro = FALSE, minClusterSize = minModuleSize)
table(dynamicMods)
這裏共劃分了19個模塊，計數值爲這些模塊中對應的基因數量。
在文獻中，這些模塊一般經過「顏色」指代，例如咱們常常會在文獻中看到以「red module」、「lightgreen module」等做爲模塊的名稱進行描述。
相似地，接下來爲這些模塊賦值顏色，並經過顏色名稱命名不一樣的模塊。
#模塊顏色指代
dynamicColors <- labels2colors(dynamicMods)
table(dynamicColors)
plotDendroAndColors(geneTree, dynamicColors, 'Dynamic Tree Cut',
dendroLabels = FALSE, addGuide = TRUE, hang = 0.03, guideHang = 0.05,
main = 'Gene dendrogram and module colors')
文獻中也常見到這種表現形式，將基因共表達類似度矩陣繪製以熱圖，和表徵模塊的聚類樹組合在一塊兒展現。
#基因表達聚類樹和共表達拓撲熱圖
plot_sim <- -(1-tom_sim)
#plot_sim <- log(tom_sim)
diag(plot_sim) <- NA
TOMplot(plot_sim, geneTree, dynamicColors,
main = 'Network heatmap plot, selected genes')
熱圖中，若基因間表達類似度越高，則顏色越深。能夠看到，同一模塊內的基因具備較高的共表達模式，而不一樣模塊間則相差明顯。
模塊特徵基因及其應用舉例
所謂模塊特徵基因，定義爲相應模塊的表達矩陣的第一主成分，表明了該模塊內基因表達的總體水平。
所以，模塊特徵基因並不是某個具體的基因，而是一種「特徵組合」。
模塊特徵基因計算
使用moduleEigengenes()可對基因表達值矩陣中的基因按模塊劃分後，執行特徵分解，獲取模塊特徵基因。
#計算基因表達矩陣中模塊的特徵基因（第一主成分）
MEList <- moduleEigengenes(gene, colors = dynamicColors)
MEs <- MEList$eigengenes
head(MEs)[1:6]
#輸出模塊特徵基因矩陣
#write.table(MEs, 'moduleEigengenes.txt', sep = '\t', col.names = NA, quote = FALSE)
模塊特徵基因可用於反映各樣本中，各基因共表達模塊（具備類似的共表達模式）中基因的總體表達水平，或者使用它們替代整個模塊中的基因執行特定的分析。
共表達模塊的進一步聚類
例如對於上述劃分的共表達模塊，若是以爲模塊過多，能夠進一步合併，將類似度較高的模塊聚合到一塊兒。此時能夠經過模塊特徵基因來實現，參考如下示例。
#經過模塊特徵基因計算模塊間相關性，表徵模塊間類似度
ME_cor <- cor(MEs)
ME_cor[1:6,1:6]
#繪製聚類樹觀察
METree <- hclust(as.dist(1-ME_cor), method = 'average')
plot(METree, main = 'Clustering of module eigengenes', xlab = '', sub = '')
#探索性分析，觀察模塊間的類似性
#height 值可表明模塊間的相異度，並肯定一個合適的閾值做爲剪切高度
#以便爲低相異度（高類似度）的模塊合併提供依據
abline(h = 0.2, col = 'blue')
abline(h = 0.25, col = 'red')
#模塊特徵基因聚類樹熱圖
plotEigengeneNetworks(MEs, '', cex.lab = 0.8, xLabelsAngle= 90,
marDendro = c(0, 4, 1, 2), marHeatmap = c(3, 4, 1, 2))
根據聚類樹中所示模塊間類似程度，定義一個合適的高度閾值，聚類樹中低於該值的模塊能夠視爲具備較高的類似程度，容許它們合併至一個簇中。
#類似模塊合併，以 0.25 做爲合併閾值（剪切高度），在此高度下的模塊將合併
#近似理解爲相關程度高於 0.75 的模塊將合併到一塊兒
merge_module <- mergeCloseModules(gene, dynamicColors, cutHeight = 0.25, verbose = 3)
mergedColors <- merge_module$colors
table(mergedColors)
#基因表達和模塊聚類樹
plotDendroAndColors(geneTree, cbind(dynamicColors, mergedColors), c('Dynamic Tree Cut', 'Merged dynamic'),
dendroLabels = FALSE, addGuide = TRUE, hang = 0.03, guideHang = 0.05)
從新裁剪聚類樹的結果中，保留了13個模塊（以前是19個）。
共表達模塊與與臨牀表型的關聯分析
上文已經提到，同一模塊內的基因具備較高的共表達類似性，即存在較高的協同，暗示它們參與了類似的調節途徑，或者在類似的細胞區域中發揮功能。
怎樣尋找具備生物學意義的共表達模塊（基因集）呢？來看一個在文獻中常常見到的例子。在腫瘤轉錄組分析中，在構建了WGCNA共表達模塊後，一般將這些模塊與患者的臨牀指標聯繫起來，以探索哪些基因間的協同做用可能與生理特徵、疾病進展等密切相關。
因爲共表達模塊是一組基因集，而非單一的變量，所以沒法直接計算其與臨牀表型資料的相關性。此時模塊特徵基因就能夠派上用場。
#患者的臨牀表型數據
trait <- read.delim('trait.txt', row.names = 1, check.names = FALSE)
#使用上一步新組合的共表達模塊的結果
module <- merge_module$newMEs
#患者基因共表達模塊和臨牀表型的相關性分析
moduleTraitCor <- cor(module, trait, use = 'p')
moduleTraitCor[1:6,1:6] #相關矩陣
#相關係數的 p 值矩陣
moduleTraitPvalue <- corPvalueStudent(moduleTraitCor, nrow(module))
#輸出相關係數矩陣或 p 值矩陣
#write.table(moduleTraitCor, 'moduleTraitCor.txt', sep = '\t', col.names = NA, quote = FALSE)
#write.table(moduleTraitPvalue, 'moduleTraitPvalue.txt', sep = '\t', col.names = NA, quote = FALSE)
#相關圖繪製
textMatrix <- paste(signif(moduleTraitCor, 2), '\n(', signif(moduleTraitPvalue, 1), ')', sep = '')
dim(textMatrix) <- dim(moduleTraitCor)
par(mar = c(5, 10, 3, 3))
labeledHeatmap(Matrix = moduleTraitCor, main = paste('Module-trait relationships'),
xLabels = names(trait), yLabels = names(module), ySymbols = names(module),
colorLabels = FALSE, colors = greenWhiteRed(50), cex.text = 0.7, zlim = c(-1,1),
textMatrix = textMatrix, setStdMargins = FALSE)
圖中，每一行表明不一樣的基因共表達模塊，每一列表明不一樣的臨牀表型，數值表明了相關係數，並分別經過紅色和綠色區分正負相關，括號中的值爲顯著性p值。
由此可幫助評估，哪些模塊中的基因可能與特定的生理特徵、疾病進展等密切相關。
模塊內基因及表達矩陣的提取
若是找到了感興趣的共表達模塊，最後就是將該模塊所涉及的基因名稱、基因的表達矩陣、或者基因間共表達類似度矩陣等提取出來，以便進行下游的分析，繼續探索它們的功能。
例如，上述相關圖中觀察到「black」模塊與腫瘤TNM分期存在很是顯著的正相關，暗示該模塊內的基因集參與了腫瘤進程。所以，指望得到與「black」模塊有關的基因信息。
#基因與模塊的對應關係列表
gene_module <- data.frame(gene_name = colnames(gene), module = mergedColors, stringsAsFactors = FALSE)
head(gene_module)
#「black」模塊內的基因名稱
gene_module_select <- subset(gene_module, module == 'black')$gene_name
#「black」模塊內基因在各樣本中的表達值矩陣（基因表達值矩陣的一個子集）
gene_select <- t(gene[,gene_module_select])
#「black」模塊內基因的共表達類似度（在 TOM 矩陣中提取子集）
tom_select <- tom_sim[gene_module_select,gene_module_select]
#輸出
#write.table(gene_select, 'gene_select.txt', sep = '\t', col.names = NA, quote = FALSE)
#write.table(tom_select, 'tom_select.txt', sep = '\t', col.names = NA, quote = FALSE)
例如後續能夠根據基因名稱，在STRING數據庫中構建這些基因的蛋白互做網絡，或者執行GO、KEGG功能富集分析或GSEA查看它們主要涉及到哪些通路的失調可能致使疾病發生和進展等，再也不多說。
尋找模塊中關鍵基因子集
已知「black」模塊中的基因表達與腫瘤TNM分期存在很是顯著的正相關，因爲該模塊中的基因數量很是多，所以指望進一步找出可能的驅動基因子集。
可經過計算各基因在「black」模塊中的模塊成員度（module membership），識別該模塊的表明基因；同時將各基因的表達值與TNM分期表型做相關性分析，看哪些基因的高表達與TNM分期的進展顯著相關。
#選擇 black 模塊內的基因
gene_black <- gene[ ,gene_module_select]
#基因的模塊成員度（module membership）計算
#即各基因表達值與相應模塊特徵基因的相關性，其衡量了基因在全局網絡中的位置
geneModuleMembership <- signedKME(gene_black, module['MEblack'], outputColumnName = 'MM')
MMPvalue <- as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nrow(module)))
#各基因表達值與臨牀表型的相關性分析
geneTraitSignificance <- as.data.frame(cor(gene_black, trait['TNM'], use = 'p'))
GSPvalue <- as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nrow(trait)))
#選擇顯著（p<0.01）、高 black 模塊成員度（MM>=0.8），與 TNM 表型高度相關（r>=0.8）的基因
geneModuleMembership[geneModuleMembership<0.8 | MMPvalue>0.01] <- 0
geneTraitSignificance[geneTraitSignificance<0.8 | MMPvalue>0.01] <- 0
select <- cbind(geneModuleMembership, geneTraitSignificance)
select <- subset(select, geneModuleMembership>=0.8 & geneTraitSignificance>=0.8)
head(select)
考慮多個條件後，識別了36個可能的候選標誌物。
此外，還可結合差別表達分析的結果，例如腫瘤組織和正常組織相比哪些基因發生顯著失調，並根據顯著性、差別倍數變化程度等繼續在這36個候選基因中進一步選擇可能推進腫瘤進程的基因。
順便做圖觀察這些基因的概況。
#候選基因與 TNM 分期的相關性散點圖
dir.create('black_TNM_cor', recursive = TRUE)
for (i in rownames(select)) {
png(paste('black_TNM_cor/', i, '-TNM.png', sep = ''),
width = 4, height = 4, res = 400, units = 'in', type = 'cairo')
plot(trait[ ,'TNM'], gene[ ,i], xlab = 'TNM', ylab = i,
main = paste('MM_black = ', select[i,'MMblack'], '\ncor_TNM = ', select[i,'TNM']), cex = 0.8, pch = 20)
fit <- lm(gene[ ,i]~trait[ ,'TNM'])
lines(trait[ ,'TNM'], fitted(fit))
dev.off()
}
#候選基因間的 TOM 矩陣
plotNetworkHeatmap(gene,
plotGenes = rownames(select),
networkType = 'unsigned',
useTOM = TRUE,
power = powers,
main = 'unsigned correlations')
關於網絡可視化
執行上述一系列分析後，也能大體理解了，WGCNA的主要側重點實際上是識別感興趣的基因集（共表達模塊）。
對於網絡的可視化是次要內容，但若是有須要，能夠將相應的結果導出爲邊列表或節點列表的形式，後續再讀入至Cytoscape等網絡可視化軟件中進行共表達網絡的可視化。
#輸出各模塊子網絡的邊列表和節點列表，後續可導入 cytoscape 繪製網絡圖
#其中，邊的權重爲 TOM 矩陣中的值，記錄了基因間共表達類似性
dir.create('cytoscape', recursive = TRUE)
for (mod in 1:nrow(table(mergedColors))) {
modules <- names(table(mergedColors))[mod]
probes <- colnames(gene)
inModule <- (mergedColors == modules)
modProbes <- probes[inModule]
modGenes <- modProbes
modtom_sim <- tom_sim[inModule, inModule]
dimnames(modtom_sim) <- list(modProbes, modProbes)
outEdge <- paste('cytoscape/', modules, '.edge_list.txt',sep = '')
outNode <- paste('cytoscape/', modules, '.node_list.txt', sep = '')
exportNetworkToCytoscape(modtom_sim,
edgeFile = outEdge,
nodeFile = outNode,
weighted = TRUE,
threshold = 0.3, #該參數可控制輸出的邊數量，過濾低權重的邊
nodeNames = modProbes,
altNodeNames = modGenes,
nodeAttr = mergedColors[inModule])
}
Cytoscape可視化不在本篇範疇內，還請自行了解了。
參考文獻
基礎統計圖（以R做圖爲主）
illustrator做圖
illustrator做圖教程2-乙烯信號轉導通路圖
網絡分析
網絡拓撲結構：節點和邊特徵
本文分享自微信公衆號 - 小明的數據分析筆記本（）。
若有侵權，請聯繫 support@oschina.cn 刪除。
本文參與「OSC源創計劃」，歡迎正在閱讀的你也加入，一塊兒分享。