小果手把手教你用GS-MM散点图挖掘WGCNA的hub基因
尔云间 一个专门做科研的团队
小果之前接触过WGCNA分析，但是浅尝辄止，只做到把基因分成模块，算出这些模块和基因之间的相关性及P值，就结束了，没有继续挖掘Hub基因，但小果也是个有求知欲的人，于是小果觉得还是要挖一下，于是就整理了下面这些代码。
温馨提示：以下代码是在掌握了WGCNA分析模块基因的基础，即已经把基因分成模块并算出这些模块和基因之间的相关性及P值的前提下学习的。
代码如下：
library('WGCNA')
fpkm=read.table(“after_meger_mdd.txt”,header=T,row.names=1,sep=”\t”, comment.char="",check.names=F)
datExpr=as.data.frame(t(fpkm[,-1]))
names(datExpr)=fpkm$Tag #第一行第一列为Tag
rownames(datExpr)=names(fpkm[,-1]) #倒置表达矩阵，行为样本，列为基因
datExpr=read.table(“性状t”,header=T,row.names=1,sep=”\t”) #读取性状文件
nGenes = ncol(datExpr)
nSamples = nrow(datExpr)
# 指定datTrait中感兴趣的一个性状，这里选择SubB
SubB = as.data.frame(datTraits$SubB)
names(SubB) = "SubB"
# 各基因模块的名字（颜色）
modNames = substring(names(MEs), 3) #这里的MEs是在WGCNA分析区分模块之后形成的变量，里面记录了模块的信息。
# 计算MM的P值
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="")
# 计算性状和基因表达量之间的相关性（GS）
geneTraitSignificance = as.data.frame(cor(datExpr, SubB, use = "p"))
GSPvalue = as.data.frame(corPvalueStudent(as.matrix(geneTraitSignificance),
nSamples))
names(geneTraitSignificance) = paste("GS.", names(SubB), sep="")
names(GSPvalue) = paste("p.GS.", names(SubB), sep="")
module = "green" #选择模块
column = match(module, modNames)
moduleGenes = moduleColors==module
green_module<-as.data.frame(dimnames(data.frame(datExpr))[[2]][moduleGenes])
names(green_module)="genename"
MM<-abs(geneModuleMembership[moduleGenes,column])
GS<-abs(geneTraitSignificance[moduleGenes, 1])
c<-as.data.frame(cbind(MM,GS)) #包含了MM和GS的数据，可以保留一下
rownames(c)=green_module$genename
green_hub <-abs(c$MM)>0.8&abs(c$GS)>0.2 #筛选hub基因
write.csv(green_hub, "hubgene_MMGS_green.csv")
verboseScatterplot(abs(geneModuleMembership[moduleGenes,
column]),abs(geneTraitSignificance[moduleGenes, 1]), xlab =
paste("Module Membership in", module, "module"), ylab = "Gene
significance for SubB", main = paste("Module membership
vs. gene significance"), pch = 20,col="grey") #画散点图
abline(h=0.2,v=0.8,col="red",lwd=1.5) #添加参考线
小伙伴们，看懂了没有，虽然看上去有点复杂，但一步一步的跟着小果走，其实也不难哦，小伙伴们有什么问题欢迎来和小果分享讨论哟。
推荐阅读
- 关注小果，小果将会持续为你带来更多生信干货哦。
生信果
生信入门、R语言、生信图解读与绘制、软件操作、代码复现、生信硬核知识技能、服务器等
原创内容
