Weighted Gene Co-Expression Network Analysis（加權(quán)基因共表達(dá)網(wǎng)絡(luò)分析锌历，簡稱WGCNA）可鑒定表達(dá)模式相似的基因集合（module）柿祈，解析基因集合與樣品表型之間的聯(lián)系，繪制基因集合中基因之間的調(diào)控網(wǎng)絡(luò)并鑒定關(guān)鍵調(diào)控基因。
WGCNA的流程主要分為：輸入數(shù)據(jù)——構(gòu)建共表達(dá)網(wǎng)絡(luò)——?jiǎng)澐帜K——模塊與性狀關(guān)聯(lián)分析——模塊之間的關(guān)聯(lián)分析——模塊中Hub基因的鑒定。此處輸入的數(shù)據(jù)是經(jīng)過其他方法篩出來的，比如之前TCseq分析中熟丸，對照和處理之間有差異的Cluster集群。具體原理網(wǎng)上有很多伪节，這里就不多說光羞，直接開搞。
1怀大、文件準(zhǔn)備
首先纱兑，就是準(zhǔn)備兩個(gè)文件：
1、基因表達(dá)矩陣文件化借，即橫軸是不同樣本以及生物重復(fù)ID潜慎，縱軸是基因ID，之后是基因的FPKM值蓖康，在此將該文件命名為“FPKM.txt”铐炫。
2、性狀數(shù)據(jù)文件蒜焊，橫軸是不同性狀的名稱倒信，縱軸是對應(yīng)樣品的名稱，格式必須與上個(gè)文件保持一致泳梆，在此將該文件命名為“trait.txt”鳖悠。
2榜掌、構(gòu)建相關(guān)性矩陣和鄰接矩陣
通過R語言中的WGCNA包開始分析，要注意的是WGCNA依賴的包較多竞穷，根據(jù)系統(tǒng)提示唐责，逐一安裝其他包。
#安裝WGCNA包
install.packages('BiocManager')
BiocManager::install('WGCNA')
library(WGCNA) #若加載失敗瘾带，安裝提示中依賴的包
#基因表達(dá)值矩陣
gene <- read.delim('FPKM.txt', row.names = 1, check.names = FALSE)
#只保留平均表達(dá)值在 1 以上的基因
gene <- subset(gene, rowSums(gene)/ncol(gene) >= 1)
#轉(zhuǎn)置矩陣
gene <- t(gene)
接著需要我們確定一個(gè)軟閾值powers（soft threshold）鼠哥，以建立鄰接矩陣并根據(jù)連通度使基因分布符合無尺度網(wǎng)絡(luò)。
#power 值選擇
powers <- 1:20
sft <- pickSoftThreshold(gene, powerVector = powers, verbose = 8)
#擬合指數(shù)與 power 值散點(diǎn)圖
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.80, col = 'red') #R方的值可以自己確定看政，一般>0.8最好
#平均連通性與 power 值散點(diǎn)圖
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')
左圖縱坐標(biāo)是R平方的值朴恳，R方越接近1，表明該網(wǎng)絡(luò)就越接近無尺度網(wǎng)絡(luò)允蚣，通常要求R方的值＞0.8于颖，負(fù)值是因?yàn)槌藄lope列數(shù)值的負(fù)方向。軟閾值一般選擇R方大于0.8時(shí)的power值嚷兔，即10森渐。如果不確定，可以在代碼中輸入“sft$powerEstimate”讓系統(tǒng)自身估計(jì)冒晰。
可以看到系統(tǒng)評估的是13同衣，emmm...當(dāng)然也要綜合上面右圖——縱坐標(biāo)是平均連通性，該值隨β的增加而降低壶运。那么接下來就用系統(tǒng)評估的最佳閾值13構(gòu)建無標(biāo)度網(wǎng)絡(luò)耐齐。
3、構(gòu)建共表達(dá)網(wǎng)絡(luò)
#估計(jì)的最佳 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 相異度 = 1 – TOM 相似度
tom_dis <- 1 - tom_sim
#層次聚類樹
geneTree <- hclust(as.dist(tom_dis), method = 'average')
plot(geneTree, xlab = '', sub = '', main = 'Gene Dendrogram',
labels = FALSE, hang = 0.04)
#使用動(dòng)態(tài)剪切樹挖掘模塊
minModuleSize <- 30 #模塊基因數(shù)目
dynamicMods <- cutreeDynamic(dendro = geneTree, distM = tom_dis,
deepSplit = 2, pamRespectsDendro = FALSE, minClusterSize = minModuleSize)
table(dynamicMods) #查看共劃分多少個(gè)模塊
該步驟是對基因進(jìn)行聚類蒋情，每條線代表一個(gè)基因埠况，相似的基因被聚到一個(gè)分支。輸入“table(dynamicMods)”后可以看到棵癣，將這些基因劃分為14個(gè)模塊辕翰。
接下來，為這些模塊增添顏色浙巫，并通過顏色名稱命名不同的模塊金蜀。
#模塊顏色指代
dynamicColors <- labels2colors(dynamicMods)
table(dynamicColors)
plotDendroAndColors(geneTree, dynamicColors, 'Dynamic Tree Cut',
dendroLabels = FALSE, addGuide = TRUE, hang = 0.03, guideHang = 0.05,
main = 'Gene dendrogram')
4、共表達(dá)拓?fù)錈釄D
#基因表達(dá)聚類樹和共表達(dá)拓?fù)錈釄D
plot_sim <- -(1-tom_sim)
diag(plot_sim) <- NA
TOMplot(plot_sim, geneTree, dynamicColors,
main = 'Network heatmap plot')
在該熱圖中的畴，基因間表達(dá)相似度越高渊抄，顏色就會(huì)越深。
#計(jì)算基因表達(dá)矩陣中模塊的特征基因
MEList <- moduleEigengenes(gene, colors = dynamicColors)
MEs <- MEList$eigengenes
head(MEs)[1:6]
#表征模塊間相似度
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 值確定一個(gè)合適的閾值作為剪切高度
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))
#相關(guān)程度高于 0.75 的模塊將合并到一起
merge_module <- mergeCloseModules(gene, dynamicColors, cutHeight = 0.25, verbose = 3)
mergedColors <- merge_module$colors
table(mergedColors)
#基因表達(dá)和模塊聚類樹
plotDendroAndColors(geneTree, cbind(dynamicColors, mergedColors), c('Dynamic Tree Cut', 'Merged dynamic'),
dendroLabels = FALSE, addGuide = TRUE, hang = 0.03, guideHang = 0.05)
5丧裁、模塊與性狀關(guān)聯(lián)
接下來就是性狀數(shù)據(jù)與基因模塊之間的關(guān)聯(lián)分析护桦。
#導(dǎo)入trait性狀數(shù)據(jù)
trait <- read.delim('trait.txt', row.names = 1, check.names = FALSE)
#使用上一步新組合的共表達(dá)模塊的結(jié)果
module <- merge_module$newMEs
#基因共表達(dá)模塊和表型的相關(guān)性分析
moduleTraitCor <- cor(module, trait, use = 'p')
moduleTraitCor[1:6,1:6] #相關(guān)矩陣，根據(jù)自己的性狀個(gè)數(shù)進(jìn)行修改
#相關(guān)系數(shù)的 p 值矩陣
moduleTraitPvalue <- corPvalueStudent(moduleTraitCor, nrow(module))
#相關(guān)圖繪制
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 = blueWhiteRed(50), cex.text = 0.7, zlim = c(-1,1),
textMatrix = textMatrix, setStdMargins = FALSE)
每一個(gè)模塊中煎娇，上面的數(shù)字代表皮爾遜相關(guān)系數(shù)二庵，下面（括號）里的值代表顯著性P值贪染。該圖表示哪些模塊中的基因可能與處理表型密切相關(guān)。從圖中不難發(fā)現(xiàn)催享，“green”模塊“MTA10”時(shí)期存在較為顯著的正相關(guān)杭隙，表示該模塊中的基因參與了MTA處理引起的表型差異，進(jìn)一步獲取該模塊中的基因（ps：該數(shù)據(jù)并不好因妙，相關(guān)性系數(shù)不是太高痰憎，只有0.7，一般最好在0.9左右）攀涵。
6铣耘、篩選關(guān)鍵基因集
#基因與模塊的對應(yīng)關(guān)系列表
gene_module <- data.frame(gene_name = colnames(gene), module = mergedColors, stringsAsFactors = FALSE)
head(gene_module)
#“green”模塊內(nèi)的基因名稱
gene_module_select <- subset(gene_module, module == 'green')$gene_name
#“green”模塊內(nèi)基因在各樣本中的表達(dá)值矩陣
gene_select <- t(gene[,gene_module_select])
#“green”模塊內(nèi)基因的共表達(dá)相似度
tom_select <- tom_sim[gene_module_select,gene_module_select]
#選擇 green 模塊內(nèi)的基因
gene_green <- gene[ ,gene_module_select]
#基因的模塊成員度計(jì)算
geneModuleMembership <- signedKME(gene_green, module['MEgreen'], outputColumnName = 'MM')
MMPvalue <- as.data.frame(corPvalueStudent(as.matrix(geneModuleMembership), nrow(module)))
#各基因表達(dá)值與表型的相關(guān)性分析
geneTraitSignificance <- as.data.frame(cor(gene_green, trait['MTA10'], use = 'p'))
GSPvalue <- as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance), nrow(trait)))
#選擇顯著（p<0.01）、高 green 模塊成員度（MM>=0.8）以故，與 MTA10表型高度相關(guān)（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)
nrow(select)#查看相關(guān)系高的基因個(gè)數(shù)
最終篩選到了16個(gè)候選基因（右側(cè)“Environment”中點(diǎn)擊“select”即可查看）蜗细，進(jìn)一步結(jié)合注釋信息確認(rèn)其是否是自己想要的hub基因。WGCNA分析過程較為繁瑣怒详，一環(huán)扣一環(huán)炉媒，要確保每一步正確，再進(jìn)行下一步操作昆烁。
