WGCNA：加權(quán)基因共表達網(wǎng)絡(luò)分析，簡而言之，就是將基因劃分為若干個模塊，探究與表型數(shù)據(jù)與基因模塊之間的相關(guān)關(guān)系，并找到模塊中的核心基因。
適用于復(fù)雜的數(shù)據(jù)模式，推薦5組(或者15個樣品)以上的數(shù)據(jù)。一般可應(yīng)用的研究方向有：不同器官或組織類型發(fā)育調(diào)控、同一組織不同發(fā)育調(diào)控、非生物脅迫不同時間點應(yīng)答、病原菌侵染后不同時間點應(yīng)答。
更多原理還有專業(yè)術(shù)語可以看生信技能樹的帖子：https://cloud.tencent.com/developer/article/1516749
1.R包下載
BiocManager::install("WGCNA")
2.數(shù)據(jù)導(dǎo)入和清洗
這是任何網(wǎng)絡(luò)分析的第一步。在這里展示了如何加載典型的表達數(shù)據(jù)，將其預(yù)處理成適合網(wǎng)絡(luò)分析的格式，并通過去除明顯的離群樣本以及缺失條目數(shù)量過多的基因和樣本來清洗數(shù)據(jù)。
2.1 導(dǎo)入表達數(shù)據(jù)
#設(shè)置工作目錄
> setwd("D:/RNA-seq/WGCNA/fpkm")
#載入WGCNA包
> library('WGCNA')
#這個設(shè)置很重要，不要忽略，允許R語言以最大線程運行
> options(stringsAsFactors = FALSE)
> allowWGCNAThreads()
Allowing multi-threading with up to 4 threads.
#導(dǎo)入表達數(shù)據(jù)
> fpkm <- read.csv(file="fpkm.csv")
#簡單查看數(shù)據(jù)
> dim(fpkm)
[1] 138749 61
請注意，每一行對應(yīng)一個基因，而每一列對應(yīng)一個樣本或輔助信息。需要移除輔助信息，因為WGCNA針對的是基因進行聚類，而一般我們的聚類是針對樣本，所以這個時候需對表達數(shù)據(jù)進行轉(zhuǎn)置以進行進一步分析。
因為我的基因數(shù)目比較多，所以在這里選取差異最大的前30%進行模塊構(gòu)建。
> rownames(fpkm) <- fpkm[,1]
> fpkm<- fpkm[,-1]
> WGCNA_matrix = t(fpkm[order(apply(fpkm,1,mad), decreasing = T)[1:41625],]) #按mad進行排序，取前30%的基因進行后續(xù)分析
> datExpr0 = as.data.frame(WGCNA_matrix)
#全部基因
> datExpr_all = as.data.frame(t(fpkm))
2.2 檢查缺失值和離群樣本
如果是全部基因，那首先檢查基因和樣本是否有過多的缺失值。如果是Top xx%選擇后的基因，那通常 gsg$allOK會返回TRUE，因為前期已經(jīng)進行了選擇。
#全部基因
> gsg_all = goodSamplesGenes(datExpr_all, verbose = 3)
> gsg_all$allOK
[1] FALSE
#篩選后基因
> gsg = goodSamplesGenes(datExpr0, verbose = 3)
> gsg$allOK
[1] TRUE
如果返回值是TURE，證明所有的基因都符合要求，如果返回值是FALSE，需要去除離群值。
> if (!gsg_all$allOK)
+ {
+ # Optionally, print the gene and sample names that were removed:
+ if (sum(!gsg_all$goodGenes)>0)
+ printFlush(paste("Removing genes:", paste(names(datExpr_all)[!gsg_all$goodGenes], collapse = ", ")));
+ if (sum(!gsg_all$goodSamples)>0)
+ printFlush(paste("Removing samples:", paste(rownames(datExpr_all)[!gsg_all$goodSamples], collapse = ", ")));
+ # Remove the offending genes and samples from the data:
+ datExpr_all = datExpr_all[gsg_all$goodSamples, gsg_all$goodGenes]
+ }
#查看新的datExpr_all
> dim(datExpr_all)
[1] 60 122752
#可以看到基因數(shù)目變少
再次檢測
> gsg_all = goodSamplesGenes(datExpr_all, verbose = 3)
> gsg_all$allOK
[1] TRUE
#這次返回值是TRUE，可以繼續(xù)后續(xù)計算
由于基因數(shù)目太多，后續(xù)計算都是使用datExpr0數(shù)據(jù)集。
對樣本進行聚類
> 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)
這里的樣本有明顯離群值，可以通過設(shè)定閾值去掉離群值。
# 畫一條閾值線，根據(jù)實際情況人工設(shè)定，這里是15
> abline(h = 15, col = "red")
# 確定閾值線下的集群
> clust = cutreeStatic(sampleTree, cutHeight = 15, minSize = 10)
> table(clust)
# clust 1包含想要留下的樣本
> keepSamples = (clust==1)
> datExpr = datExpr0[keepSamples, ]
現(xiàn)在變量datExpr包含用于網(wǎng)絡(luò)分析的表達量數(shù)據(jù)
2.3 導(dǎo)入表型數(shù)據(jù)
> datTraits <- read.csv("phe_blue.csv",header=TRUE)
> dim(datTraits)
[1] 60 5
#去掉第一列樣本名
> rownames(datTraits) <- datTraits[,1]
> datTraits <- datTraits[,-1]
> dim(datTraits)
[1] 60 4
現(xiàn)在表型數(shù)據(jù)和基因表達量數(shù)據(jù)都已經(jīng)成功導(dǎo)入，在繼續(xù)網(wǎng)絡(luò)構(gòu)建之前，可以可視化樣本和表型之間的關(guān)系。
# 將表型值轉(zhuǎn)換為顏色表示:白色表示低，紅色表示高，灰色表示缺失
> traitColors = numbers2colors(datTraits, signed = FALSE)
#在樣本聚類圖的基礎(chǔ)上，增加表型值熱圖
> plotDendroAndColors(sampleTree, traitColors,
+ groupLabels = names(datTraits),
+ main = "Sample dendrogram and trait heatmap")
2.4 保存數(shù)據(jù)
最后一步是保存處理好的表達量數(shù)據(jù)和表型數(shù)據(jù)，以便在本教程的后續(xù)步驟中使用。
save(datExpr, datTraits, file = "WGCNA0.3-dataInput.RData")
