WGCNA.jpg

WGCNA 簡明指南|1. 基因共表達網(wǎng)絡(luò)構(gòu)建及模塊識別

[TOC]

參考

本文主要參考官方指南Tutorials for WGCNA R package (ucla.edu)混驰，詳細(xì)內(nèi)容可參閱官方文檔龙填。

其它資料：

1. WGCNA - 文集 - 簡書 (jianshu.com)
2. WGCNA分析，簡單全面的最新教程 - 簡書 (jianshu.com)

簡介

加權(quán)基因共表達網(wǎng)絡(luò)分析 (WGCNA, Weighted correlation network analysis)是用來描述不同樣品之間基因關(guān)聯(lián)模式的系統(tǒng)生物學(xué)方法锤窑，可以用來鑒定高度協(xié)同變化的基因集, 并根據(jù)基因集的內(nèi)連性和基因集與表型之間的關(guān)聯(lián)鑒定候補生物標(biāo)記基因或治療靶點涧卵。簡而言之，就是將基因劃分為若干個模塊，探究表型數(shù)據(jù)與基因模塊之間的相關(guān)關(guān)系吼句，并找到模塊中的hub基因。

22

數(shù)據(jù)導(dǎo)入事格、清洗及預(yù)處理

數(shù)據(jù)導(dǎo)入

示例數(shù)據(jù)下載：https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/Tutorials/FemaleLiver-Data.zip

# BiocManager::install("WGCNA") 
library('WGCNA')
# 在讀入數(shù)據(jù)時惕艳，遇到字符串之后搞隐，不將其轉(zhuǎn)換為factors，仍然保留為字符串格式
options(stringsAsFactors = FALSE)
# 導(dǎo)入示例數(shù)據(jù)
femData = read.csv("LiverFemale3600.csv")
# 查看數(shù)據(jù)
dim(femData)
names(femData)

> dim(femData)
[1] 3600  143

> names(femData)
  [1] "substanceBXH"   "gene_symbol"    "LocusLinkID"    "ProteomeID"    
  [5] "cytogeneticLoc" "CHROMOSOME"     "StartPosition"  "EndPosition"   
  [9] "F2_2"           "F2_3"           "F2_14"          "F2_15"  
  ...

# 提取樣本-基因表達矩陣
datExpr0 = as.data.frame(t(femData[, -c(1:8)]))
names(datExpr0) = femData$substanceBXH
rownames(datExpr0) = names(femData)[-c(1:8)]

> datExpr0[1:6,1:6]
      MMT00000044 MMT00000046 MMT00000051 MMT00000076 MMT00000080 MMT00000102
F2_2   -0.0181000     -0.0773 -0.02260000    -0.00924 -0.04870000  0.17600000
F2_3    0.0642000     -0.0297  0.06170000    -0.14500  0.05820000 -0.18900000
F2_14   0.0000644      0.1120 -0.12900000     0.02870 -0.04830000 -0.06500000
F2_15  -0.0580000     -0.0589  0.08710000    -0.04390 -0.03710000 -0.00846000
F2_19   0.0483000      0.0443 -0.11500000     0.00425  0.02510000 -0.00574000
F2_20  -0.1519741     -0.0938 -0.06502607    -0.23610  0.08504274 -0.01807182

檢查過度缺失值和離群樣本

# 檢查缺失值太多的基因和樣本
gsg = goodSamplesGenes(datExpr0, verbose = 3);
gsg$allOK

如果最后一個語句返回TRUE远搪，則所有的基因都通過了檢查劣纲。如果沒有，我們就從數(shù)據(jù)中剔除不符合要求的基因和樣本谁鳍。

if(!gsg$allOK)
{
  #(可選)打印被刪除的基因和樣本名稱:
  if(sum(!gsg$goodGenes)>0)
    printFlush(paste("Removinggenes:",paste(names(datExpr0)[!gsg$goodGenes], collapse =",")));
  if(sum(!gsg$goodSamples)>0)
    printFlush(paste("Removingsamples:",paste(rownames(datExpr0)[!gsg$goodSamples], collapse =",")));
  #刪除不滿足要求的基因和樣本:
  datExpr0 = datExpr0[gsg$goodSamples, gsg$goodGenes]
}

接下來癞季，我們對樣本進行聚類(與稍后的聚類基因相反)，看看是否有任何明顯的異常值倘潜。

sampleTree = hclust(dist(datExpr0), method ="average");
# 繪制樣本樹:打開一個尺寸為12 * 9英寸的圖形輸出窗口
# 可對窗口大小進行調(diào)整
sizeGrWindow(12,9)
# 如要保存可運行下面語句
# pdf(file="Plots/sampleClustering.pdf",width=12,height=9);
par(cex = 0.6)
par(mar =c(0,4,2,0))
plot(sampleTree, main ="Sampleclusteringtodetectoutliers",sub="", xlab="", cex.lab = 1.5,
cex.axis= 1.5, cex.main = 2)

Fig1a.樣本樹狀圖

Fig1a可見有一個異常值绷柒。可以手動或使用自動方法刪除它窍荧。選擇一個高度(height)進行切割將刪除異常樣本辉巡，比如15(Fig1b)，并在該高度使用分支切割

# 繪制閾值切割線
abline(h = 15,col="red");
# 確定閾值線下的集群
clust = cutreeStatic(sampleTree, cutHeight = 15, minSize = 10)
table(clust)
# clust1包含想要留下的樣本.
keepSamples = (clust==1)
datExpr = datExpr0[keepSamples, ]
nGenes =ncol(datExpr)
nSamples =nrow(datExpr)

Fig1b.紅線表示切割高度

datExpr包含用于網(wǎng)絡(luò)分析的表達數(shù)據(jù)蕊退。

載入臨床特征數(shù)據(jù)

將樣本信息與臨床特征進行匹配郊楣。

traitData =read.csv("ClinicalTraits.csv")
dim(traitData)
names(traitData)
# 刪除不必要的列.
allTraits = traitData[, -c(31, 16)]
allTraits = allTraits[,c(2, 11:36) ]
dim(allTraits)
names(allTraits)
# 形成一個包含臨床特征的數(shù)據(jù)框
femaleSamples =rownames(datExpr)
traitRows =match(femaleSamples, allTraits$Mice)
datTraits = allTraits[traitRows, -1]
rownames(datTraits) = allTraits[traitRows, 1]
collectGarbage() # 釋放內(nèi)存

現(xiàn)在在variabledatExpr中有了表達數(shù)據(jù)，在變量datTraits中有了相應(yīng)的臨床特征瓤荔。在進行網(wǎng)絡(luò)構(gòu)建和模塊檢測之前净蚤，將臨床特征與樣本樹狀圖的關(guān)系可視化。

# 重新聚類樣本
sampleTree2 = hclust(dist(datExpr), method ="average")
# 將臨床特征值轉(zhuǎn)換為連續(xù)顏色:白色表示低输硝，紅色表示高今瀑，灰色表示缺失
traitColors = numbers2colors(datTraits, signed = FALSE);
# 在樣本聚類圖的基礎(chǔ)上，增加臨床特征值熱圖
plotDendroAndColors(sampleTree2, traitColors,
                    groupLabels =names(datTraits),
                    main ="Sample dendrogramand trait heatmap")

Figure 2: Clustering dendrogram of samples based on their Euclidean distance.

自動構(gòu)建網(wǎng)絡(luò)及識別模塊

確定合適的軟閾值：網(wǎng)絡(luò)拓?fù)浞治?/h3>

The soft thresholding, is a value used to power the correlation of the genes to that threshold. The assumption on that by raising the correlation to a power will reduce the noise of the correlations in the adjacency matrix. To pick up one threshold use the pickSoftThreshold function, which calculates for each power if the network resembles to a scale-free graph. The power which produce a higher similarity with a scale-free network is the one you should use.

WGCNA: What is soft thresholding? (bioconductor.org)

# 設(shè)置軟閾值調(diào)參范圍
powers =c(c(1:10),seq(from = 12, to=20,by=2))
# 網(wǎng)絡(luò)拓?fù)浞治?sft = pickSoftThreshold(datExpr, powerVector = powers, verbose = 5)
# 繪圖
sizeGrWindow(9, 5)
# 1行2列排列
par(mfrow =c(1,2));
cex1 = 0.9;
# 無標(biāo)度拓?fù)鋽M合指數(shù)與軟閾值的函數(shù)(左圖)
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     xlab="SoftThreshold(power)",ylab="ScaleFreeTopologyModelFit,signedR^2",type="n",
     main =paste("Scaleindependence"));
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     labels=powers,cex=cex1,col="red");
# 這條線對應(yīng)于h的R^2截止點
abline(h=0.90,col="red")
# Mean Connectivity與軟閾值的函數(shù)(右圖)
plot(sft$fitIndices[,1], sft$fitIndices[,5],
     xlab="SoftThreshold(power)",ylab="MeanConnectivity", type="n",
     main =paste("Meanconnectivity"))
text(sft$fitIndices[,1], sft$fitIndices[,5],labels=powers, cex=cex1,col="red")

Figure 3:各種軟閾值的網(wǎng)絡(luò)拓?fù)浞治?/div>

我們選擇power 6点把，這是無標(biāo)度拓?fù)鋽M合指數(shù)曲線在達到較高值時趨于平緩的最低 power橘荠。

一步構(gòu)建網(wǎng)絡(luò)和識別模塊

net = blockwiseModules(datExpr,power= 6,
                       TOMType ="unsigned", minModuleSize = 30,
                       reassignThreshold = 0, mergeCutHeight = 0.25,
                       numericLabels = TRUE, pamRespectsDendro = FALSE,
                       saveTOMs = TRUE,
                       saveTOMFileBase ="femaleMouseTOM",
                       verbose = 3)

• deepSplit 參數(shù)調(diào)整劃分模塊的敏感度，值越大郎逃，越敏感哥童，得到的模塊就越多，默認(rèn)是2褒翰；
• minModuleSize 參數(shù)設(shè)置最小模塊的基因數(shù)贮懈，值越小，小的模塊就會被保留下來优训；
• mergeCutHeight 設(shè)置合并相似性模塊的距離朵你，值越小，就越不容易被合并揣非，保留下來的模塊就越多

# 查看識別了多少模塊以及模塊大小
table(net$colors)

> table(net$colors)

  0   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18 
 99 609 460 409 316 312 221 211 157 123 106 100  94  91  77  76  58  47  34 

表明有18個模塊抡医，按大小遞減順序標(biāo)記為1到18，大小從609到34個基因早敬。標(biāo)簽0為所有模塊之外的基因魂拦。

# 可視化模塊
sizeGrWindow(12, 9)
# 將標(biāo)簽轉(zhuǎn)換為顏色
mergedColors = labels2colors(net$colors)
# 繪制樹狀圖和模塊顏色圖
plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
                    "Modulecolors",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)

Figure 4:基因的聚類樹狀圖毛仪，具有不同的基于拓?fù)渲丿B，以及指定的模塊顏色

保存模塊賦值和模塊特征基因信息芯勘，以供后續(xù)分析箱靴。

moduleLabels = net$colors
moduleColors = labels2colors(net$colors)
MEs = net$MEs;
geneTree = net$dendrograms[[1]];
save(MEs, moduleLabels, moduleColors, geneTree,
     file="FemaleLiver-02-networkConstruction-auto.RData")

往期

1. 跟著Nature學(xué)作圖 | 配對啞鈴圖+分組擬合曲線+分類變量熱圖
2. (免費教程+代碼領(lǐng)取)|跟著Cell學(xué)作圖系列合集
3. 跟著Nat Commun學(xué)作圖 | 1.批量箱線圖+散點+差異分析
4. 跟著Nat Commun學(xué)作圖 | 2.時間線圖
5. 跟著Nat Commun學(xué)作圖 | 3.物種豐度堆積柱狀圖
6. 跟著Nat Commun學(xué)作圖 | 4.配對箱線圖+差異分析

---