找完差異基因后，現(xiàn)在開始這篇文章的WGCNA分析昼扛，這個分析最主要的是要得到與表型相關(guān)的基因贵扰；
1.png
先看看文章中怎么做的：
- 首先對數(shù)據(jù)進行過濾：過濾后為5435個基因；
2.png
構(gòu)建基因共表達網(wǎng)絡(luò)：軟閾值（20）坚冀；
3.png
4.png
模塊與臨床特征的相關(guān)性分析；
5.png
下面是代碼部分：
一 數(shù)據(jù)過濾
rm(list = ls())
options(stringsAsFactors = F)
load('./Rdata/exp_group.Rdata')
#提取文章中提到的25% of the maximum variationin基因
data=as.data.frame(t(data))
m.vars=apply(data,2,var)
expro.upper=data[,which(m.vars>quantile(m.vars,probs = seq(0,1,0.25))[4])]
data=as.data.frame(t(expro.upper))
dim(data)
6.png
文章中為5435個基因鉴逞；
二 構(gòu)建基因共表達網(wǎng)絡(luò)
library(WGCNA)
#1.判斷數(shù)據(jù)質(zhì)量
gsg <- goodSamplesGenes(data, verbose = 3)
gsg$allOK #這個數(shù)據(jù)結(jié)果是TRUE记某，如果質(zhì)量不好，可以去百度下怎么處理构捡；
#2.按文章中的設(shè)置從12開始間隔2
powers = c(c(1:10), seq(from = 12, to=30, by=2))
dat = as.data.frame(t(data))#注意倒置液南；
sft = pickSoftThreshold(dat, powerVector = powers, verbose = 5)
#3.畫文章中的βvalue圖
png("./picture/beta-value.png",width = 800,height = 600)
par(mfrow = c(1,2))
cex1 = 0.8
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.80,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")
dev.off()
7.png
這里可以發(fā)現(xiàn)軟閾值結(jié)果與文章中一致，選擇20勾徽；
#4.構(gòu)建共表達網(wǎng)絡(luò)
#設(shè)置一下power和maxBlockSize兩個數(shù)值滑凉，power就是上面得到的最先達到0.8的數(shù)值20；
net = blockwiseModules(dat, power = 20, maxBlockSize = 5439,
TOMType ='unsigned', minModuleSize = 100, #minModuleSize這個參數(shù)需要自己設(shè)置模塊最小基因數(shù)
reassignThreshold = 0, mergeCutHeight = 0.25,
numericLabels = TRUE, pamRespectsDendro = FALSE,
saveTOMs=TRUE, saveTOMFileBase = "drought",
verbose = 3)
table(net$colors)
8.png
這里有10個模塊，文章中是11個畅姊，這可能是由于minModuleSize設(shè)置的數(shù)值不一樣咒钟，文章中好像沒有標明這個參數(shù)的數(shù)值；
#5.繪制基因聚類和模塊顏色組合
moduleLabels = net$colors
moduleColors = labels2colors(moduleLabels)
png("./picture/genes-modules.png",width = 800,height = 600)
plotDendroAndColors(net$dendrograms[[1]], moduleColors[net$blockGenes[[1]]],
"Module colors",
dendroLabels = FALSE, hang = 0.03,
addGuide = TRUE, guideHang = 0.05)
dev.off()
9.png
三 模塊與臨床特征的相關(guān)性分析
這一步是WGCNA最重要的結(jié)果若未；
#1.導(dǎo)入臨床信息（需要整理好的GSE140797_clinical.csv數(shù)據(jù)可以留言給我）
data_clin=read.csv('./data/GSE140797_clinical.csv',header = T)
rownames(data_clin)=data_clin$Sample_geo_accession
data_clin=data_clin[,-1]
data_clin=as.data.frame(t(data_clin))
data_clin$normal = ifelse(grepl('normal',data_clin$Sample_characteristics_ch1),1,0)
data_clin$tumor = ifelse(grepl('adenocarcinoma',data_clin$Sample_characteristics_ch1),1,0)
data_clin=data_clin[,-1]
#2.模塊與性狀相關(guān)性
nGenes = ncol(dat)
nSamples = nrow(dat)
# 獲取eigengenes朱嘴，用顏色標簽計算ME值
MEs0 = moduleEigengenes(dat, moduleColors)$eigengenes
MEs = orderMEs(MEs0) ##不同顏色的模塊的ME值矩 (樣本vs模塊)
#計算每個模塊和每個性狀之間的相關(guān)性
moduleTraitCor = cor(MEs, data_clin , use = "p")
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)
# 可視化相關(guān)性和P值
textMatrix = paste(signif(moduleTraitCor, 2), "\n(",
signif(moduleTraitPvalue, 1), ")", sep = "");
dim(textMatrix) = dim(moduleTraitCor)
png("./picture/Module-trait-relationships.png",width = 800,height = 1200,res = 120)
par(mar = c(6, 8.5, 3, 3));
labeledHeatmap(Matrix = moduleTraitCor,
xLabels = colnames(data_clin),
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-trait relationships"))
dev.off()
10.png
我們的結(jié)果是選取最正相關(guān)的yellow模塊和最負相關(guān)的turquoise模塊，分別為0.95和-0.84粗合；
#3.對模塊里的具體基因分析
#跟tumor最相關(guān)的是yellow模塊和turquoise模塊萍嬉；
nSamples=nrow(dat)
#計算MM值和GS值
modNames = substring(names(MEs), 3) ##切割，從第三個字符開始保存
## 算出每個模塊跟基因的皮爾森相關(guān)系數(shù)矩陣
geneModuleMembership = as.data.frame(cor(dat, MEs, use = "p"))
##計算MM值對應(yīng)的P值
MMPvalue = as.data.frame(
corPvalueStudent(as.matrix(geneModuleMembership), nSamples)
)
names(geneModuleMembership) = paste("MM", modNames, sep="")
names(MMPvalue) = paste("p.MM", modNames, sep="")
##計算基因與每個性狀的顯著性（相關(guān)性）及pvalue值
geneTraitSignificance = as.data.frame(cor(dat, data_clin, use = "p"))
GSPvalue = as.data.frame(
corPvalueStudent(as.matrix(geneTraitSignificance), nSamples)
)
names(geneTraitSignificance) = paste("GS.", colnames(data_clin), sep="")
names(GSPvalue) = paste("p.GS.", colnames(data_clin), sep="")
#3.1 分析yellow模塊
module = "yellow"
column = match(module, modNames) ##在所有模塊中匹配選擇的模塊隙疚，返回所在的位置
yellow_moduleGenes <- names(net$colors)[which(moduleColors == module)]
MM <- abs(geneModuleMembership[yellow_moduleGenes, column])
GS <- abs(geneTraitSignificance[yellow_moduleGenes, 1])
#畫圖
png("./picture/Module_yellow_membership_gene_significance.png",width = 800,height = 600)
par(mfrow = c(1,1))
verboseScatterplot(MM,
GS,
xlab = paste("Module Membership in", module, "module"),
ylab = "Gene significance for Basal",
main = paste("Module membership vs. gene significance\n"),
cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
dev.off()
#提取HUB基因：MM > 0.8,GS > 0.8
yellow_hub <- yellow_moduleGenes[(GS > 0.8 & MM > 0.8)]
length(yellow_hub)
#[1] 258
#保存至表格帚湘，后面做富集分析需要
write.csv(yellow_hub,'yellow_hub_gene.csv',
row.names = F,sep = '\t',quote = F)
#3.2 分析turquoise模塊
module = "turquoise"
column = match(module, modNames) ##在所有模塊中匹配選擇的模塊，返回所在的位置
turquoise_moduleGenes <- names(net$colors)[which(moduleColors == module)]
MM <- abs(geneModuleMembership[turquoise_moduleGenes, column])
GS <- abs(geneTraitSignificance[turquoise_moduleGenes, 1])
#畫圖
png("./picture/Module_turquoise_membership_gene_significance.png",width = 800,height = 600)
par(mfrow = c(1,1))
verboseScatterplot(MM,
GS,
xlab = paste("Module Membership in", module, "module"),
ylab = "Gene significance for Basal",
main = paste("Module membership vs. gene significance\n"),
cex.main = 1.2, cex.lab = 1.2, cex.axis = 1.2, col = module)
dev.off()
#提取HUB基因：MM > 0.8,GS > 0.8
turquoise_hub <- turquoise_moduleGenes[(GS > 0.8 & MM > 0.8)]
length(turquoise_hub)
#[1] 690
#保存至表格甚淡，后面做富集分析需要
write.csv(turquoise_hub,'turquoise_hub_gene.csv',
row.names = F,sep = '\t',quote = F)
小結(jié)：
通過WGCNA分析，我們找到了與腫瘤表型密切相關(guān)的yellow模塊和turquoise模塊捅厂，并得到了這兩個模塊里的基因yellow_hub_gene.csv和turquoise_hub_gene.csv表格贯卦。
往期文章復(fù)現(xiàn)：
