小果手把手教你用GS-MM散點(diǎn)圖挖掘WGCNA的hub基因
爾云間? 一個(gè)專門做科研的團(tuán)隊(duì)
??
小果之前接觸過WGCNA分析，但是淺嘗輒止，只做到把基因分成模塊，算出這些模塊和基因之間的相關(guān)性及P值，就結(jié)束了，沒有繼續(xù)挖掘Hub基因，但小果也是個(gè)有求知欲的人，于是小果覺得還是要挖一下，于是就整理了下面這些代碼。
溫馨提示：以下代碼是在掌握了WGCNA分析模塊基因的基礎(chǔ)，即已經(jīng)把基因分成模塊并算出這些模塊和基因之間的相關(guān)性及P值的前提下學(xué)習(xí)的。
代碼如下：
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]) #倒置表達(dá)矩陣，行為樣本，列為基因
datExpr=read.table(“性狀t”,header=T,row.names=1,sep=”\t”) #讀取性狀文件
nGenes = ncol(datExpr)
nSamples = nrow(datExpr)
# 指定datTrait中感興趣的一個(gè)性狀，這里選擇SubB
SubB = as.data.frame(datTraits$SubB)
names(SubB) = "SubB"
#? 各基因模塊的名字（顏色）
modNames = substring(names(MEs), 3) #這里的MEs是在WGCNA分析區(qū)分模塊之后形成的變量，里面記錄了模塊的信息。
# 計(jì)算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="")
# 計(jì)算性狀和基因表達(dá)量之間的相關(guān)性（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的數(shù)據(jù)，可以保留一下
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") #畫散點(diǎn)圖
abline(h=0.2,v=0.8,col="red",lwd=1.5) #添加參考線
小伙伴們，看懂了沒有，雖然看上去有點(diǎn)復(fù)雜，但一步一步的跟著小果走，其實(shí)也不難哦，小伙伴們有什么問題歡迎來和小果分享討論喲。
推薦閱讀
JASPAR——可預(yù)測(cè)轉(zhuǎn)錄因子DNA結(jié)合蛋白結(jié)合識(shí)別位點(diǎn)的數(shù)據(jù)庫(kù)
關(guān)注小果，小果將會(huì)持續(xù)為你帶來更多生信干貨哦。
生信果
生信入門、R語(yǔ)言、生信圖解讀與繪制、軟件操作、代碼復(fù)現(xiàn)、生信硬核知識(shí)技能、服務(wù)器等
原創(chuàng)內(nèi)容
