爾雲間 一個專門做科研的團隊

原創 小果 生信果

小果之前接觸過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) #添加參考線

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

推薦閱讀

• GEO2R分析R代碼學習之差異分析結果可視化
• 使用R語言完成序列比對及進化樹美化
• 你不知道的PCA及在R中的實現
• 柱狀圖-腫瘤某一指標的比較和GSVA結果展示
• JASPAR——可預測轉錄因子DNA結合蛋白結合識別位點的資料庫

關注小果，小果將會持續為你帶來更多生信乾貨哦。

生信果生信入門、R語言、生信圖解讀與繪製、軟體操作、代碼復現、生信硬核知識技能、伺服器等原創內容