WGCNA 系列
[TOC]
22
參考
本文主要參考官方指南Tutorials for WGCNA R package (ucla.edu)，詳細(xì)內(nèi)容可參閱官方文檔。
其它資料：
- WGCNA - 文集 - 簡(jiǎn)書(shū) (jianshu.com)
- WGCNA分析，簡(jiǎn)單全面的最新教程 - 簡(jiǎn)書(shū) (jianshu.com)
- WGCNA：（加權(quán)共表達(dá)網(wǎng)絡(luò)分析）_bioprogrammer-CSDN博客
- WGCNA如何從module中挖掘關(guān)鍵基因_廬州月光的博客-CSDN博客
關(guān)聯(lián)模塊與臨床特征
如果前面沒(méi)保存再跑一遍WGCNA 簡(jiǎn)明指南|1. 基因共表達(dá)網(wǎng)絡(luò)構(gòu)建及模塊識(shí)別以得到本次分析需要的數(shù)據(jù)。
量化module-trait(模塊-特征)關(guān)系
在這個(gè)分析中咧栗，我們希望識(shí)別與臨床特征顯著相關(guān)的模塊交煞。由于我們已經(jīng)有了每個(gè)模塊的特征基因，我們簡(jiǎn)單地將特征基因與臨床特征關(guān)聯(lián)起來(lái)，并尋找最顯著的關(guān)聯(lián)凤粗。
**Module eigengene **
Module eigengene is defined as the first principal component of the expression matrix of the corresponding module.
# 導(dǎo)入之前的數(shù)據(jù)(也可以重新跑一遍上一期的內(nèi)容)
lnames =load(file="FemaleLiver-02-networkConstruction-auto.RData");
lnames
# 明確基因和樣本數(shù)
nGenes =ncol(datExpr);
nSamples =nrow(datExpr);
# 用顏色標(biāo)簽重新計(jì)算MEs
# 按照模塊計(jì)算每個(gè)module的MEs（也就是該模塊的第一主成分）
MEs0 = moduleEigengenes(datExpr, moduleColors)$eigengenes
# 對(duì)給定的(特征)向量重新排序异逐，使相似的向量(通過(guò)相關(guān)性測(cè)量)彼此相鄰。
MEs = orderMEs(MEs0)
# 計(jì)算基因模塊MEs 與 臨床特征的相關(guān)性以及p值
# use 給出在缺少值時(shí)計(jì)算協(xié)方差的方法
moduleTraitCor =cor(MEs, datTraits, use ="p");
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)
可視化
# 如需保存
# pdf("Module-trait associations.pdf",width = 8, height=10)
par(mar =c(6, 8.5, 3, 3));
# 通過(guò)熱圖顯示相關(guān)性
labeledHeatmap(Matrix = moduleTraitCor,
xLabels =names(datTraits),
yLabels =names(MEs),
ySymbols =names(MEs),
colorLabels = FALSE,
colors= greenWhiteRed(50),
textMatrix = textMatrix,
setStdMargins = FALSE,
cex.text= 0.5,
zlim =c(-1,1),
main =paste("Module-traitrelationships"))
# dev.off()
基因與性狀和重要模塊的關(guān)系：基因重要性和模塊成員
Module Membership
簡(jiǎn)稱(chēng)MM, 將該基因的表達(dá)量與module的第一主成分，即module eigengene進(jìn)行相關(guān)性分析就可以得到MM值，所以MM值本質(zhì)上是一個(gè)相關(guān)系數(shù)晃琳，如果基因和某個(gè)module的MM值為0顾翼，說(shuō)明二者根本不相關(guān)，該基因不屬于這個(gè)module; 如果MM的絕對(duì)值接近1，說(shuō)明基因與該module相關(guān)性很高蕊肥。Gene Significance
簡(jiǎn)稱(chēng)GS, 將該基因的表達(dá)量與對(duì)應(yīng)的表型數(shù)值進(jìn)行相關(guān)性分析展东，最終的相關(guān)系數(shù)的值就是GS, GS反映出基因表達(dá)量與表型數(shù)據(jù)的相關(guān)性爪膊。
# 以特征變量weight 為例
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="");
模塊內(nèi)分析:鑒定高GS和MM基因
在感興趣的基因模塊中小槐，篩選出高GS和MM的基因疮方。在圖1中可以看到曾掂，與性狀weight_g相關(guān)性最高的基因模塊是MEbrown，因此選擇MEbrown進(jìn)行后續(xù)分析蝴猪，繪制brown模塊中GS和MM的散點(diǎn)圖米酬。
module ="brown"
column =match(module, modNames);
moduleGenes = moduleColors==module;
sizeGrWindow(7, 7);
par(mfrow =c(1,1));
verboseScatterplot(abs(geneModuleMembership[moduleGenes, column]),
abs(geneTraitSignificance[moduleGenes, 1]),
xlab =paste("ModuleMembershipin", module,"module"),
ylab ="Genesignificanceforbodyweight",
main =paste("Modulemembershipvs.genesignificance\n"),
cex.main = 1.2, cex.lab = 1.2, cex.axis= 1.2,col= module)
如圖2所示，GS和MM是高度相關(guān)的丧慈，這說(shuō)明與一個(gè)性狀高度相關(guān)的基因往往也是與該性狀相關(guān)的模塊中最重要的(中心)元素完域。
網(wǎng)絡(luò)分析結(jié)果匯總輸出
我們已經(jīng)找到了與我們的興趣特征高度相關(guān)的模塊姿现，并通過(guò)MM確定了它們的核心參與者。現(xiàn)在，我們將這些統(tǒng)計(jì)信息與基因注釋合并并導(dǎo)出匹涮。
names(datExpr)
# 返回brown模塊中所有ID
names(datExpr)[moduleColors=="brown"]
# 導(dǎo)入注釋文件
annot = read.csv(file = "GeneAnnotation.csv");
dim(annot)
names(annot)
probes = names(datExpr)
probes2annot = match(probes, annot$substanceBXH)
# 統(tǒng)計(jì)沒(méi)有注釋到的基因
sum(is.na(probes2annot))
# 創(chuàng)建數(shù)據(jù)集雳攘，包含探測(cè)ID 刑巧，
geneInfo0 = data.frame(substanceBXH = probes,
geneSymbol = annot$gene_symbol[probes2annot],
LocusLinkID = annot$LocusLinkID[probes2annot],
moduleColor = moduleColors,
geneTraitSignificance, GSPvalue)
# 通過(guò)顯著性對(duì)模塊進(jìn)行排序
modOrder = order(-abs(cor(MEs, weight, use = "p")))
# 添加模塊成員
for (mod in 1:ncol(geneModuleMembership))
{
oldNames = names(geneInfo0)
geneInfo0 = data.frame(geneInfo0, geneModuleMembership[,
modOrder[mod]],
MMPvalue[, modOrder[mod]]);
names(geneInfo0) = c(oldNames, paste("MM.",
modNames[modOrder[mod]], sep=""),
paste("p.MM.", modNames[modOrder[mod]], sep=""))
}
# 對(duì)基因進(jìn)行排序
geneOrder = order(geneInfo0$moduleColor, abs(geneInfo0$GS.weight))
geneInfo = geneInfo0[geneOrder, ]
# 導(dǎo)出
write.csv(geneInfo, file = "geneInfo.csv")
模塊基因富集分析
批量導(dǎo)出模塊基因后拯辙，使用在線網(wǎng)站或者R包clusterProfiler進(jìn)行富集分析。
# Get the corresponding Locuis Link IDs
allLLIDs = annot$LocusLinkID[probes2annot];
# 選擇感興趣的模塊
intModules = c("brown", "red", "salmon")
for (module in intModules)
{
# 模塊探針
modGenes = (moduleColors==module)
# 得到 entrez ID
modLLIDs = allLLIDs[modGenes];
# 導(dǎo)出
fileName = paste("LocusLinkIDs-", module, ".txt", sep="");
write.table(as.data.frame(modLLIDs), file = fileName,
row.names = FALSE, col.names = FALSE)
}
往期
- 跟著Nature學(xué)作圖 | 配對(duì)啞鈴圖+分組擬合曲線+分類(lèi)變量熱圖
- (免費(fèi)教程+代碼領(lǐng)取)|跟著Cell學(xué)作圖系列合集
- 跟著Nat Commun學(xué)作圖 | 1.批量箱線圖+散點(diǎn)+差異分析
- 跟著Nat Commun學(xué)作圖 | 2.時(shí)間線圖
- 跟著Nat Commun學(xué)作圖 | 3.物種豐度堆積柱狀圖
- 跟著Nat Commun學(xué)作圖 | 4.配對(duì)箱線圖+差異分析
