? ?

小云之間接觸過WGCNA分析，但是淺嘗輒止，只做到把基因分成模塊，算出這些模塊和基因之間的相關性及P值，就結束了，沒有繼續(xù)挖掘Hub基因，但小云也是個有求知欲的人，于是小云覺得還是要挖一下，于是就整理了下面這些代碼。小云溫馨提示，以下代碼是在掌握了WGCNA分析模塊基因的基礎，即已經(jīng)把基因分成模塊并算出這些模塊和基因之間的相關性及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分析區(qū)分模塊之后形成的變量，里面記錄了模塊的信息。

# 計算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的數(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") #畫散點圖

abline(h=0.2,v=0.8,col="red",lwd=1.5) #添加參考線

? ? 小伙伴們，看懂了沒有，雖然看上去有點復雜，但一步一步的跟著小云走，其實也不難哦，小伙伴們有什么問題歡迎來和小云分享討論喲。