找完差异基因后，现在开始这篇文章的WGCNA分析，这个分析最主要的是要得到与表型相关的基因；
1.png
先看看文章中怎么做的：
- 首先对数据进行过滤：过滤后为5435个基因；
2.png
构建基因共表达网络：软阈值（20）；
3.png
4.png
模块与临床特征的相关性分析；
5.png
下面是代码部分：
一 数据过滤
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个基因；
二 构建基因共表达网络
library(WGCNA)
#1.判断数据质量
gsg <- goodSamplesGenes(data, verbose = 3)
gsg$allOK #这个数据结果是TRUE，如果质量不好，可以去百度下怎么处理；
#2.按文章中的设置从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
这里可以发现软阈值结果与文章中一致，选择20；
#4.构建共表达网络
#设置一下power和maxBlockSize两个数值，power就是上面得到的最先达到0.8的数值20；
net = blockwiseModules(dat, power = 20, maxBlockSize = 5439,
TOMType ='unsigned', minModuleSize = 100, #minModuleSize这个参数需要自己设置模块最小基因数
reassignThreshold = 0, mergeCutHeight = 0.25,
numericLabels = TRUE, pamRespectsDendro = FALSE,
saveTOMs=TRUE, saveTOMFileBase = "drought",
verbose = 3)
table(net$colors)
8.png
这里有10个模块，文章中是11个，这可能是由于minModuleSize设置的数值不一样，文章中好像没有标明这个参数的数值；
#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
三 模块与临床特征的相关性分析
这一步是WGCNA最重要的结果；
#1.导入临床信息（需要整理好的GSE140797_clinical.csv数据可以留言给我）
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.模块与性状相关性
nGenes = ncol(dat)
nSamples = nrow(dat)
# 获取eigengenes，用颜色标签计算ME值
MEs0 = moduleEigengenes(dat, moduleColors)$eigengenes
MEs = orderMEs(MEs0) ##不同颜色的模块的ME值矩 (样本vs模块)
#计算每个模块和每个性状之间的相关性
moduleTraitCor = cor(MEs, data_clin , use = "p")
moduleTraitPvalue = corPvalueStudent(moduleTraitCor, nSamples)
# 可视化相关性和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
我们的结果是选取最正相关的yellow模块和最负相关的turquoise模块，分别为0.95和-0.84；
#3.对模块里的具体基因分析
#跟tumor最相关的是yellow模块和turquoise模块；
nSamples=nrow(dat)
#计算MM值和GS值
modNames = substring(names(MEs), 3) ##切割，从第三个字符开始保存
## 算出每个模块跟基因的皮尔森相关系数矩阵
geneModuleMembership = as.data.frame(cor(dat, MEs, use = "p"))
##计算MM值对应的P值
MMPvalue = as.data.frame(
corPvalueStudent(as.matrix(geneModuleMembership), nSamples)
)
names(geneModuleMembership) = paste("MM", modNames, sep="")
names(MMPvalue) = paste("p.MM", modNames, sep="")
##计算基因与每个性状的显著性（相关性）及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)
小结：
通过WGCNA分析，我们找到了与肿瘤表型密切相关的yellow模块和turquoise模块，并得到了这两个模块里的基因yellow_hub_gene.csv和turquoise_hub_gene.csv表格。
往期文章复现：
