這篇文章更多的是對于混亂的中文資源的梳理，并補充了一些沒有提到的重要參數(shù)，希望大家不會踩坑。
1. 簡介
1.1 背景
WGCNA（weighted gene co-expression network analysis，權(quán)重基因共表達網(wǎng)絡(luò)分析）是一種分析多個樣本基因表達模式的分析方法，可將表達模式相似的基因進行聚類，并分析模塊與特定性狀或表型之間的關(guān)聯(lián)關(guān)系，因此在基因組研究中被廣泛應(yīng)用。
相比于只關(guān)注差異表達的基因，WGCNA利用數(shù)千或近萬個變化最大的基因或全部基因的信息識別感興趣的基因集，并與表型進行顯著性關(guān)聯(lián)分析。既充分利用了信息，也把數(shù)千個基因與表型的關(guān)聯(lián)轉(zhuǎn)換為數(shù)個基因集與表型的關(guān)聯(lián)，免去了多重假設(shè)檢驗校正的問題。
WGCNA算法首先假定基因網(wǎng)絡(luò)服從無尺度分布(scale free network)，并定義基因共表達相關(guān)矩陣、基因網(wǎng)絡(luò)形成的鄰接函數(shù)，然后計算不同節(jié)點的相異系數(shù)，并據(jù)此構(gòu)建分層聚類樹(hierarchical clustering tree)，該聚類樹的不同分支代表不同的基因模塊(module)，模塊內(nèi)基因共表達程度高，而分屬不同模塊的基因共表達程度低。
1.2 無尺度網(wǎng)絡(luò)
網(wǎng)絡(luò)的數(shù)學(xué)名稱是圖，在圖論中對于每一個節(jié)點有一個重要概念，即：度(degree)。一個點的度是指圖中該點所關(guān)聯(lián)的邊數(shù)。如下圖，如果不加以思考，人們很容易認為生活中常見的網(wǎng)絡(luò)會是一種random network，即每一個節(jié)點的度相對平均。然而第二種圖，即scale-free network才是一種更穩(wěn)定的選擇。Scale-free network具有這樣的特點，即存在少數(shù)節(jié)點具有明顯高于一般點的度，這些點被稱為hub。由少數(shù)hub與其它節(jié)點關(guān)聯(lián)，最終構(gòu)成整個網(wǎng)絡(luò)。這樣的網(wǎng)絡(luò)的節(jié)點度數(shù)與具有該度數(shù)的節(jié)點個數(shù)間服從power distribution。生物體選擇scale-free network而不是random network尤其進化上的原因，對于scale-free network，少數(shù)關(guān)鍵基因執(zhí)行主要功能，這種網(wǎng)絡(luò)具有非常好的魯棒性(Robust)，即只要保證hub的完整性，整個生命體的基本活動在一定刺激影響下將不會受到太大影響，而random network若受到外界刺激，其受到的傷害程度將直接與刺激強度成正比。
1.3 相關(guān)術(shù)語
共表達網(wǎng)絡(luò)：點代表基因，邊代表基因表達相關(guān)性。加權(quán)是指對相關(guān)性值進行冥次運算 (冥次的值也就是軟閾值 (power, pickSoftThreshold這個函數(shù)所做的就是確定合適的power))。無向網(wǎng)絡(luò)(unsigned network)的邊屬性計算方式為
abs(cor(genex, geney)) ^ power；有向網(wǎng)絡(luò)(signed network)的邊屬性計算方式為(1+cor(genex, geney)/2) ^ power; sign hybrid的邊屬性計算方式為cor(genex, geney)^power if cor>0 else 0， sign hybrid意味著它既包含加權(quán)網(wǎng)絡(luò)也包含非加權(quán)網(wǎng)絡(luò)。這種處理方式強化了強相關(guān)，弱化了弱相關(guān)或負相關(guān)，使得相關(guān)性數(shù)值更符合無標(biāo)度網(wǎng)絡(luò)特征，更具有生物意義。除了軟閾值還有硬閾值一說，計算方式是 a_ij = 1 if s_ij > β otherwise a_ij = 0。這里的β就是硬閾值(hard threshold)。Module(模塊)：高度內(nèi)連的基因集。在無向網(wǎng)絡(luò)中，模塊內(nèi)是高度相關(guān)的基因。在有向網(wǎng)絡(luò)中，模塊內(nèi)是高度正相關(guān)的基因。
Connectivity (連接度)：類似于網(wǎng)絡(luò)中 “度” (degree)的概念。每個基因的連接度是與其相連的基因的
邊屬性之和。Module eigengene E: 給定模型的第一主成分，代表整個模型的基因表達譜。即用一個向量代替了一個矩陣，方便后期計算。
Intramodular connectivity: 給定基因與給定模型內(nèi)其他基因的關(guān)聯(lián)度，判斷基因所屬關(guān)系。
Adjacency matrix (鄰接矩陣)：基因和基因之間的加權(quán)相關(guān)性值構(gòu)成的矩陣。
TOM (Topological overlap matrix)：把鄰接矩陣轉(zhuǎn)換為拓撲重疊矩陣，以降低噪音和假相關(guān)，獲得的新距離矩陣，這個信息可拿來構(gòu)建網(wǎng)絡(luò)或繪制TOM圖。
2. 一般步驟
3. 代碼
利用WGCNA有一步法建網(wǎng)絡(luò)的，也有step by step的方法，為了保證理解，建議至少過一遍step by step。
安裝WGCNA根據(jù)不同的操作系統(tǒng)可能略有不同，我在macOS下安裝著實廢了一番功夫。這一部分不提。
# 加載必須的包并做參數(shù)設(shè)置
library(MASS)
library(class)
library(cluster)
library(impute)
library(Hmisc)
library(WGCNA)
options(stringsAsFactors = F)
讀取基因表達數(shù)據(jù)，注意要做一個轉(zhuǎn)換，保證基因在列，樣品在行，官方推薦使用Deseq2的varianceStabilizingTransformation或log2(x+1)對標(biāo)準(zhǔn)化后的數(shù)據(jù)做個轉(zhuǎn)換。如果數(shù)據(jù)來自不同的批次，需要先移除批次效應(yīng)。如果數(shù)據(jù)存在系統(tǒng)偏移，需要做下quantile normalization。一般要求樣本數(shù)多于15個。樣本數(shù)多于20時效果更好，樣本越多，結(jié)果越穩(wěn)定。
dat0 <- read.csv("./01raw_data/GBM55and65and8000.csv")
datExprdataOne <- t(dat0[,15:69])
datExprdataTwo <- t(dat0[, 70:134])
datSummary <- dat0[, c(1:14)]
dim(datExprdataOne); dim(datExprdataTwo); dim(datSummary)
rm(dat0)
#[1] 55 8000
#[1] 65 8000
#[1] 8000 14
檢驗數(shù)據(jù)質(zhì)量
gsg = goodSamplesGenes(datExprdataOne, verbose = 3)
if (!gsg$allOK){
# Optionally, print the gene and sample names that were removed:
if (sum(!gsg$goodGenes)>0)
printFlush(paste("Removing genes:",
paste(names(datExprdataOne)[!gsg$goodGenes], collapse = ",")));
if (sum(!gsg$goodSamples)>0)
printFlush(paste("Removing samples:",
paste(rownames(datExprdataOne)[!gsg$goodSamples], collapse = ",")));
# Remove the offending genes and samples from the data:
datExprdataOne = datExprdataOne[gsg$goodSamples, gsg$goodGenes]
}
#Flagging genes and samples with too many missing values...
# ..step 1
檢查是否有離群值，結(jié)果顯示無
sampleTree = hclust(dist(datExprdataOne), method = "average")
plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="")
篩選軟閾值， 無向網(wǎng)絡(luò)在power小于15或有向網(wǎng)絡(luò)power小于30內(nèi)，沒有一個power值可以使無標(biāo)度網(wǎng)絡(luò)圖譜結(jié)構(gòu)R^2達到0.8或平均連接度降到100以下，則可能是由于部分樣品與其他樣品差別太大造成的。這可能由批次效應(yīng)、樣品異質(zhì)性或?qū)嶒灄l件對表達影響太大等造成，需要移除。
powers1 <- c(seq(1, 10, by=1), seq(12, 20, by=2))
sft <- pickSoftThreshold(datExprdataOne, powerVector = powers1)
RpowerTable <- pickSoftThreshold(datExprdataOne, powerVector = powers1)[[2]]
cex1 = 0.7
par(mfrow = c(1,2))
plot(RpowerTable[,1], -sign(RpowerTable[,3])*RpowerTable[,2], xlab = "soft threshold (power)", ylab = "scale free topology model fit, signes R^2", type = "n")
text(RpowerTable[,1], -sign(RpowerTable[,3])*RpowerTable[,2], labels = powers1, cex = cex1, col = "red")
abline(h = 0.95, 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")
橫軸是Soft threshold (power)，縱軸是無標(biāo)度網(wǎng)絡(luò)的評估參數(shù)，數(shù)值越高，網(wǎng)絡(luò)越符合無標(biāo)度特征 (non-scale)。
我們可以使用系統(tǒng)給的參數(shù)幫助我們得到soft threshold，可以是
sft$powerEstimate
#4
給出的是4，因為這個篩選的標(biāo)準(zhǔn)是R-square=0.85，但是我們希望R-square的值達到0.9，所以選擇power值為6.
利用power=6計算connectivity，并且可視化無尺度網(wǎng)絡(luò)圖的拓撲結(jié)構(gòu)
betal = 6
k.dataOne <- softConnectivity(datExprdataOne, power = betal) -1
k.dataTwo <- softConnectivity(datExprdataTwo, power = betal) - 1
par(mfrow=c(2,2))
scaleFreePlot(k.dataOne, main = paste("data set I, power=", betal), truncated = F)
scaleFreePlot(k.dataTwo, main = paste("data set II, power=", betal), truncated = F)
篩選連通性最高的3600個基因做為后續(xù)分析，不過一般不在這一步進行篩選，因為生物體內(nèi)的基因網(wǎng)絡(luò)圖更多的是無尺度的。
kCut <- 3601
kRank <- rank(-k.dataOne)
vardataOne <- apply(datExprdataOne, 2, var)
vardataTwo <- apply(datExprdataTwo, 2, var)
restK <- kRank <= kCut & vardataOne >0 & vardataTwo > 0
ADJdataOne <- adjacency(datExpr = datExprdataOne[,restK], power = betal)
dissTOMdataOne <- TOMdist(ADJdataOne)
hierTOMdataOne <- hclust(as.dist(dissTOMdataOne), method = "average")
par(mfrow = c(1,1))
plot(hierTOMdataOne, labels = F, main = "dendrogram, 3600 mast connected in data set I")
層級聚類樹展示各個模塊, 灰色的為未分類到模塊的基因，這里使用的cutreeStaticColor檢測module，但是對于復(fù)雜的基因結(jié)構(gòu)建議使用 cutreeDynamic。其中也有一些具體的參數(shù)可以選擇得到合適的module。
colordataOne <- cutreeStaticColor(hierTOMdataOne, cutHeight = 0.94, minSize = 125)
par(mfrow=c(2,1), mar = c(2,4,1,1))
plot(hierTOMdataOne, main = "Glioblastoma data set 1, n = 55", labels = F, xlab = "", sub = "")
plotColorUnderTree(hierTOMdataOne, colors = data.frame(module = colordataOne))
title("module membership data set I")
后續(xù)換了一種方法希望得到更多module以期得到更多的eigegene。
dataOne_net = blockwiseModules(datExprdataOne, power = 6, maxBlockSize = 200,
TOMType = "signed", minModuleSize = 30,
reassignThreshold = 0, mergeCutHeight = 0.25,
numericLabels = TRUE, pamRespectsDendro = FALSE,
saveTOMs=TRUE, corType = "pearson",
loadTOMs=TRUE,
saveTOMFileBase = "DataI.tom",
verbose = 3)
# Calculating module eigengenes block-wise from all genes
# Flagging genes and samples with too many missing values...
# ..step 1
# ....pre-clustering genes to determine blocks..
# Projective K-means:
# ..k-means clustering..
# ..merging smaller clusters...
# Block sizes:
繪制模塊之間相關(guān)性圖
dataOne_MEs <- dataOne_net$MEs
plotEigengeneNetworks(dataOne_MEs, "Eigengene adjacency heatmap",
marDendro = c(3,3,2,4),
marHeatmap = c(3,4,2,2), plotDendrograms = T,
xLabelsAngle = 90)
可視化基因網(wǎng)絡(luò) (TOM plot)
load(dataOne_net$TOMFiles[1], verbose=T)
## Loading objects:
## TOM
TOM <- as.matrix(TOM)
dissTOM = 1-TOM
# Transform dissTOM with a power to make moderately strong
# connections more visible in the heatmap
plotTOM = dissTOM^7
# Set diagonal to NA for a nicer plot
diag(plotTOM) = NA
# Call the plot function
TOMplot(plotTOM, dataOne_net$dendrograms,
main = "Network heatmap plot, all genes")
導(dǎo)出網(wǎng)絡(luò)圖
probes = colnames(dat0[,15:69])
dimnames(TOM) <- list(probes, probes)
# Export the network into edge and node list files Cytoscape can read
cyt = exportNetworkToCytoscape(TOM,
edgeFile = "edges.txt",
nodeFile = "nodes.txt",
weighted = TRUE, threshold = 0,
nodeNames = probes, nodeAttr = dataOne_net$colors)
部分參考：
