文章目錄

• 1 介紹
• 2 名詞解釋
• • Co-expression network
  • Module
  • Connectivity
  • Intramodular connectivity $k_{IM}$
  • Module eigengene E
  • Eigengene signi?cance
  • Module Membership $k_{ME}$
  • Hub gene
  • Gene signi?cance GS
  • Module signi?cance
• 3 流程
• • 3.1 數據輸入、清洗、預處理
  • • 3.1.1 輸入表達數據
    • 3.1.2 檢查基因和樣品是否有過多缺失值，以及異常樣品
    • 3.1.3 輸入臨床特征數據
  • 3.2 一步構建共表達網絡和劃分模塊
  • • 3.2.1 篩選軟閥值 $\beta$
    • 3.2.2 構建網絡和劃分模塊
  • 3.3 模塊與性狀關聯分析
  • • 3.3.1 量化模塊特征關聯
    • 3.3.2 基因與性狀和重要模塊的關系：基因重要性和模塊成員
    • 3.3.3 模內分析：鑒定具有高GS和MM的基因
• 4 總結
• 參考資料

1 介紹

WGCNA相比于差異表達基因可以獲得跟多信息，通過考慮測得的轉錄本之間的關系可以更完整地表示微列陣數據，這可以通過基因表達譜之間的成對相關性進行評估。WGCNA從數千個基因的層次開始，確定臨床上感興趣的基因模塊，最后使用模塊內連接性，基因重要性來識別疾病途徑中的關鍵基因，以進行進一步驗證。與其將成千上萬的基因與微陣列樣品特征相關聯，不如著重于幾個（通常少于10個）模塊與樣品特征之間的關系。為此，它計算每個模塊的特征基因顯著性（樣本特征與特征基因之間的相關性）和相應的p值。 模塊定義不使用先驗定義的基因集。取而代之的是，通過使用層次聚類從表達式數據構造模塊。 盡管建議將生成的模塊與基因本體信息相關聯以評估其生物學合理性，但這不是必需的。由于模塊可能對應于生物學途徑，因此將分析重點放在模塊內Hub基因（或模塊特征基因）上就等于一種出于生物學目的數據降維。由于模塊內hub基因的表達譜高度相關，通常會產生數十種候選生物標記。盡管這些候選者在統計學上是等效的，但它們可能在生物學上的合理性或臨床效用之間存在差異。基因本體信息可用于進一步確定模塊內Hub基因的優先級。圖一顯示了WGCNA分析流程。

• 構建共表達網絡
• 劃分模塊
• 模塊與性狀關聯分析
• 模塊之間關聯分析
• 模塊中核心基因鑒定
  

2 名詞解釋

Co-expression network

將共表達網絡定義為無方向的加權基因網絡。 這種網絡的節點對應于基因表達譜，基因之間的連線由基因表達之間的成對相關性決定。通過將相關系數的絕對值提高到冪 $\beta >1$ （軟閾值）， 加權基因共表達網絡的構建強調了高相關性，卻以低相關性為代價。具體來說 $a_{ij}=|cor(x_i,x_j){|}^{\beta }$ 表示無符號（unsigned）網絡的連接，另外，也可以指定有符號的加權共表達網絡，如 $a_{ij}=|(1+cor(x_i,x_j))/2{|}^{\beta }$。

Module

模塊是高度互連的基因簇。 在無符號共表達網絡中，模塊對應于具有高度絕對相關性的基因簇。 在有符號網絡中，模塊對應于正相關的基因。

Connectivity

對于每個基因，連通性定義為與其他網絡基因的連接強度之和 $k_i={\Sigma }_{\mu =i}{a}_{\mu i}$，在共表達網絡中，連通性可衡量基因與所有其他網絡基因之間的相關性。

Intramodular connectivity $k_{IM}$

模內連通性測量給定基因相對于特定模塊的基因如何連接或共表達。 模內連通性可以解釋為模塊成員資格的量度。

Module eigengene E

模塊特征向量E被定義為給定模塊的主成分一。 可以認為是模塊中基因表達譜的代表。

Eigengene signi?cance

當可獲得微陣列樣品特征y（例如病例對照狀態或體重）時，可以將模塊特征基因與此結果相關聯。 相關系數稱為特征基因顯著性

Module Membership $k_{ME}$

特征基因連通性，對于每個基因，我們通過將其基因表達譜與給定模塊的模塊本征基因相關聯來定義模塊成員的“模糊”度量。如，$M{M}^{blue}(i)=K_{cor,i}^{blue}=cor(x_i,E^{blue})$ 用來測量基因 $i$ 與藍色模塊特征基因的相關程度。如果 $M{M}^{blue}(i)$ 接近0，則第 $i$ 個基因不屬于blue模塊的一部分。 另一方面，如果 $M{M}^{blue}(i)$ 接近于1或 -1，則它與藍色模塊基因高度相關。模塊成員的符號編碼該基因與藍色模塊本征基因是正相關還是負相關， 可以為所有輸入基因定義模塊成員度量（無論其原始模塊成員如何）。 事實證明，模塊成員資格度量與模塊內連接性 $k_{IM}$ 高度相關。 高度連接的模塊內Hub基因傾向于對各個模塊具有較高的模塊成員值。

Hub gene

這個松散定義的術語被用作“高度連接的基因”的縮寫。通過定義，共表達模塊內部的基因往往具有高度的連通性。

Gene signi?cance GS

為了將外部信息整合到共表達網絡中，我們利用了基因顯著性方法。抽象地說，$G{S}_i$ 的絕對值越高，第 $i$ 個基因的生物學意義就越大。例如，$G{S}_i$ 可以編碼通路成員（例如，如果該基因是已知的凋亡基因，則為1，否則為0），敲除必需性或與外部微陣列樣品性狀的相關性。基因顯著性度量也可以通過減去p值的對數來定義。 唯一的要求是，基因顯著性0表示該基因對于所關注的生物學問題不重要。 基因顯著性可以取正值或負值。

Module signi?cance

模塊顯著性被定義為給定模塊中所有基因的平均絕對基因顯著性的度量。 當將基因顯著性定義為基因表達與外部性狀y的相關性時，此度量往往與模塊特征基因與y的相關性高度相關。

3 流程

3.1 數據輸入、清洗、預處理

表達數據和性狀數據下載

3.1.1 輸入表達數據

library(WGCNA) options(stringsAsFactors = FALSE) femData = read.csv("LiverFemale3600.csv") datExpr0 = as.data.frame(t(femData[, -c(1:8)])) names(datExpr0) = femData$substanceBXH

3.1.2 檢查基因和樣品是否有過多缺失值，以及異常樣品

檢測缺失值

gsg = goodSamplesGenes(datExpr0, verbose = 3) gsg$allOK [1] TRUE

如果輸出語句為TRUE，不需要刪除基因和樣本數據；如果沒有返回TURE，需要從數據中刪除潛在的基因和樣本。

if (!gsg$allOK) {if (sum(!gsg$goodGenes)>0)printFlush(paste("Removing genes:", paste(names(datExpr0)[!gsg$goodGenes], collapse = ", ")))if (sum(!gsg$goodSamples)>0)printFlush(paste("Removing samples:", paste(names(datExpr0)[!gsg$goodSamples], collapse = ", ")))datExpr0 = datExpr0[gsg$goodSamples, gsg$goodGenes] }

對樣本進行聚類，看是否有異常樣本

sampleTree = hclust(dist(datExpr0), method = "average") sizeGrWindow(12,9) par(cex = 0.6) par(mar = c(0,4,2,0)) plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5, cex.axis = 1.5, cex.main = 2) abline(h = 15, col = 'red')

從上圖可以看出，樣品F2_221樣品疑似異樣，考慮去除。

clust = cutreeStatic(sampleTree, cutHeight = 15, minSize = 10) table(clust) clust0 1 1 134

去除異常樣品后的表達矩陣。

datExpr = datExpr0[clust == 1, ] nGenes = ncol(datExpr) nSamples = nrow(datExpr)

3.1.3 輸入臨床特征數據

traitData = read.csv("ClinicalTraits.csv") dim(traitData) [1] 361 38 names(traitData) [1] "X" "Mice" "Number" "Mouse_ID" [5] "Strain" "sex" "DOB" "parents" [9] "Western_Diet" "Sac_Date" "weight_g" "length_cm" [13] "ab_fat" "other_fat" "total_fat" "comments" [17] "X100xfat_weight" "Trigly" "Total_Chol" "HDL_Chol" [21] "UC" "FFA" "Glucose" "LDL_plus_VLDL" [25] "MCP_1_phys" "Insulin_ug_l" "Glucose_Insulin" "Leptin_pg_ml" [29] "Adiponectin" "Aortic.lesions" "Note" "Aneurysm" [33] "Aortic_cal_M" "Aortic_cal_L" "CoronaryArtery_Cal" "Myocardial_cal" [37] "BMD_all_limbs" "BMD_femurs_only"

去除不需要的臨床特征和樣品

allTraits = traitData[, -c(31, 16)] allTraits = allTraits[, c(2, 11:36)] datTraits = allTraits[match(rownames(datExpr),allTraits$Mice),] rownames(datTraits) = datTraits[,1] datTraits = datTraits[,-1]

對表達矩陣再次聚類

sampleTree2 = hclust(dist(datExpr), method = "average")

將臨床特征大小用顏色深淺表示：白色表示數值低，紅色表示數值高，灰色表示缺少輸入

traitColors = numbers2colors(datTraits, signed = FALSE)

樣本聚類熱圖

plotDendroAndColors(sampleTree2, traitColors, groupLabels = names(datTraits), main = "Sample dendrogram and trait heatmap")

3.2 一步構建共表達網絡和劃分模塊

3.2.1 篩選軟閥值 $\beta$

設置一個軟閥值向量

powers = c(c(1:10), seq(from = 12, to=20, by=2))

使用網絡拓撲分析函數

sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)

軟閾值 $\beta$ 與無尺度網絡評價系數 $R^2$ 的關系，以及軟閾值 $\beta$ 與平均連通性的關系

sizeGrWindow(9, 5) par(mfrow = c(1,2)) cex1 = 0.9 plot(sft$fitIndices[,1],-sign(sft$fitIndices[,3])*sft$fitIndices[,2], xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n", main = paste("Scale independence")) text(sft$fitIndices[,1],-sign(sft$fitIndices[,3])*sft$fitIndices[,2], labels=powers,cex=cex1,col="red") abline(h=0.90,col="red") 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, cex=cex1,col="red")

$R^2$ 為無尺度網絡評價系數，一般設置為0.9，$\beta$ 取值標準：$R^2$ 第一次到達0.9時對應的 $\beta$ 值，由上圖可知，$\beta$ 選為6

3.2.2 構建網絡和劃分模塊

構建共表達網絡、劃分模塊、合并相似模塊。

net = blockwiseModules(datExpr,power = 6, # 軟閥值選為6TOMType = "unsigned", # 構建無尺度網絡minModuleSize = 30, # 最小模塊基因數為30reassignThreshold = 0, mergeCutHeight = 0.25, # 模塊合并閥值numericLabels = F, # 模塊顏色標簽pamRespectsDendro = FALSE, saveTOMs = TRUE, # 保存TOM矩陣saveTOMFileBase = "femaleMouseTOM",verbose = 3,maxBlockSize = 5000) # 可以處理的數據集的最大基因數，默認5000 # 所有模塊個數和各個模塊中基因數量 table(net$colors) black blue brown cyan green greenyellow grey 211 460 409 77 312 100 99 grey60 lightcyan lightgreen magenta midnightblue pink purple 47 58 34 123 76 157 106 red salmon tan turquoise yellow 221 91 94 609 316

結果顯示有18個模塊，模塊大小為34至609個基因，模塊grey表示所有模塊外部的基因。

moduleColors = net$colors # 模塊顏色標簽 MEs = net$MEs # 模塊特征向量

TOM矩陣層次聚類結果，其中dissTom = 1-Tom，后面進一步介紹TOM矩陣。

geneTree = net$dendrograms[[1]] geneTree Call: fastcluster::hclust(d = as.dist(dissTom), method = "average") Cluster method : average Number of objects: 3600

TOM矩陣的層次聚類以及基因所屬模塊可視化

sizeGrWindow(12, 9) mergedColors = labels2colors(net$colors) plotDendroAndColors(net$dendrograms[[1]],moduleColors,"Module colors",dendroLabels = FALSE,hang = 0.03,addGuide = TRUE,guideHang = 0.05)

TOM矩陣介紹：
首先由表達矩陣計算鄰接矩陣 $a_{ij}=|cor(x_i,x_j){|}^{\beta }$ ，$x_i$ 和 $x_j$表示任意兩個基因 ，為了最大程度地減少噪聲和虛假關聯的影響，需要將鄰接矩陣轉換為拓撲重疊矩陣（TOM矩陣），常用TOMsimilarity()函數將鄰接矩陣轉為TOM矩陣。

3.3 模塊與性狀關聯分析

3.3.1 量化模塊特征關聯

計算模塊特征向量和臨床性狀之間相關系數矩陣，并對相關系數進行檢驗。

MEs = orderMEs(MEs) moduleTraitCor = cor(MEs, datTraits, use = "p") moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)

可視化

pdf(file = 'labeledHeatmap.pdf',width = 12,height = 8) textMatrix = paste(signif(moduleTraitCor, 2), "\n(", signif(moduleTraitPvalue, 1), ")", sep = "") dim(textMatrix) = dim(moduleTraitCor) labeledHeatmap(Matrix = moduleTraitCor, xLabels = names(datTraits),yLabels = names(MEs), ySymbols = names(MEs),colorLabels = FALSE, colors = greenWhiteRed(50), textMatrix = textMatrix, setStdMargins = FALSE, cex.text = 0.7, zlim = c(-1,1), cex.lab.y = 0.6,cex.lab.x = 0.7,yColorWidth = strheight("M")/2,main = paste("Module-trait relationships")) dev.off()

3.3.2 基因與性狀和重要模塊的關系：基因重要性和模塊成員

通過將基因顯著性GS定義為基因和性狀之間的相關性的絕對值，我們可以量化單個基因與我們感興趣的性狀（體重）的關聯。 對于每個模塊，我們還定義了模塊成員MM的定量度量，作為模塊特征基因與基因表達譜的相關性。 這樣能夠量化微陣列上所有基因與每個模塊的相似性。

weight = as.data.frame(datTraits$weight_g) names(weight) = "weight" modNames = substring(names(MEs), 3) geneModuleMembership = as.data.frame(cor(datExpr, MEs, use = "p")) MMPvalue = as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nSamples)) names(geneModuleMembership) = paste("MM", modNames, sep="") names(MMPvalue) = paste("p.MM", modNames, sep="") geneTraitSignificance = as.data.frame(cor(datExpr, weight, use = "p")) GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nSamples)) names(geneTraitSignificance) = paste("GS.", names(weight), sep="") names(GSPvalue) = paste("p.GS.", names(weight), sep="")

3.3.3 模內分析：鑒定具有高GS和MM的基因

module = "brown" column = match(module, modNames) moduleGenes = moduleColors==module table(moduleGenes) sizeGrWindow(7, 7) par(mfrow = c(1,1)) verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),abs(geneTraitSignificance[moduleGenes, 1]),xlab = paste("Module Membership in", module, "module"),ylab = "Gene significance for body weight",main = paste("Module membership vs. gene significance\n"),cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)

4 總結

加權共表達網絡最主要的兩個步驟：

• 共表達網絡構建和模塊劃分
• 模塊和性狀關聯分析

通過上面兩個步驟找到與目標性狀顯著相關的模塊，對模塊里面的基因進行后續分析，如結合差異表達分析。

參考資料

WGCNA參考手冊
本博客內容將同步更新到個人微信公眾號：生信玩家。歡迎大家關注~~~

總結

以上是生活随笔為你收集整理的WGCNA：（加权共表达网络分析）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。