WGCNA的使用( 五 ) _基因

<- chooseTopHubInEachModule(datExpr,#WGCNA分析输入的表达矩阵moduleColors)#模块颜色信息# 保存hub基因结果write.table (HubGenes,file = "data/Step04-HubGenes_of_each_module.xls",quote=F,sep='\t',col.names = F)#### 与某种特征相关的hub基因NS = networkScreening(datTraits$M_stage,#M分期MEs,#datExpr)#WGCNA分析输入的表达矩阵# 将结果写入到文件write.table(NS,file="data/Step04-Genes_for_M_stage.xls",quote=F,sep='\t')## 模块GO/KEGG分析# 加载R包library(anRichment)library(clusterProfiler)##### GO分析# 构建GO背景基因集GOcollection = buildGOcollection(organism = "human")geneNames = colnames(datExpr)# 将基因SYMBOL转换为ENTREZID基因名geneID = bitr(geneNames,fromType = "SYMBOL", toType = "ENTREZID", OrgDb = "org.Hs.eg.db", drop = FALSE)# 将基因名对应结果写入文件中write.table(geneID, file = "data/Step04-geneID_map.xls", sep = "\t", quote = TRUE, row.names = FALSE)# 进行GO富集分析GOenr = enrichmentAnalysis(classLabels = moduleColors,#基因所在的模块信息identifiers = geneID$ENTREZID,refCollection = GOcollection,useBackground = "given",threshold = 1e-4,thresholdType = "Bonferroni",getOverlapEntrez = TRUE,getOverlapSymbols = TRUE,ignoreLabels = "grey");# 提取结果，并写入结果到文件tab = GOenr$enrichmentTablenames(tab)write.table(tab, file = "data/Step04-GOEnrichmentTable.xls", sep = "\t", quote = TRUE, row.names = FALSE)# 提取主要结果，并写入文件keepCols = c(1, 3, 4, 6, 7, 8, 13)screenTab = tab[, keepCols]# 小数位为2位numCols = c(4, 5, 6)screenTab[, numCols] = signif(apply(screenTab[, numCols], 2, as.numeric), 2)# 给结果命名colnames(screenTab) = c("module", "GOID", "term name", "p-val", "Bonf", "FDR", "size")rownames(screenTab) = NULL# 查看结果head(screenTab)# 写入文件中write.table(screenTab, file = "data/Step04-GOEnrichmentTableScreen.xls", sep = "\t", quote = TRUE, row.names = FALSE)##### KEGG富集分析# AnRichment没有直接提供KEGG数据的背景集# 这里使用MSigDBCollection构建C2通路数据集KEGGCollection = MSigDBCollection("data/msigdb_v7.1.xml", MSDBVersion = "7.1",organism = "human",excludeCategories = c("h","C1","C3","C4","C5","C6","C7")) # KEGG分析KEGGenr = enrichmentAnalysis(classLabels = moduleColors,identifiers = geneID$ENTREZID,refCollection = KEGGCollection,useBackground = "given",threshold = 1e-4,thresholdType = "Bonferroni",getOverlapEntrez = TRUE,getOverlapSymbols = TRUE,ignoreLabels = "grey")# 提取KEGG结果，并写入文件tab = KEGGenr$enrichmentTablenames(tab)write.table(tab, file = "data/Step04-KEGGEnrichmentTable.xls", sep = "\t", quote = TRUE, row.names = FALSE)# 提取主要结果并写入文件keepCols = c(1, 3, 4, 6, 7, 8, 13)screenTab = tab[, keepCols]# 取两位有效数字numCols = c(4, 5, 6)screenTab[, numCols] = signif(apply(screenTab[, numCols], 2, as.numeric), 2)# 对结果表格进行重命名colnames(screenTab) = c("module", "ID", "term name", "p-val", "Bonf", "FDR", "size")rownames(screenTab) = NULL# 查看结果head(screenTab)# 写入文件中write.table(screenTab, file = "data/Step04-KEGGEnrichmentTableScreen.xls", sep = "\t", quote = TRUE, row.names = FALSE)### 输出cytoscape可视化# 重新计算TOM ， power值设置为前面选择好的TOM = TOMsimilarityFromExpr(datExpr, power = 7)# 输出全部网络模块cyt = exportNetworkToCytoscape(TOM,edgeFile = "data/Step04-CytoscapeInput-edges-all.txt",#基因间的共表达关系nodeFile = "data/Step04-CytoscapeInput-nodes-all.txt",#weighted = TRUE,threshold = 0.1,nodeNames = geneID$SYMBOL,altNodeNames = geneID$ENTREZID,nodeAttr = moduleColors)# 输出感兴趣网络模块modules = c("brown", "red")# 选择上面模块中包含的基因inModule = is.finite(match(moduleColors, modules))modGenes = geneID[inModule,]# 选择指定模块的TOM矩阵modTOM = TOM[inModule, inModule]# 输出为Cytoscape软件可识别格式cyt = exportNetworkToCytoscape(modTOM,edgeFile = paste("data/Step04-CytoscapeInput-edges-", paste(modules, collapse="-"), ".txt", sep=""),nodeFile = paste("data/Step04-CytoscapeInput-nodes-", paste(modules, collapse="-"), ".txt", sep=""),weighted = TRUE,threshold = 0.2,nodeNames = modGenes$SYMBOL,altNodeNames = modGenes$ENTREZID,nodeAttr = moduleColors[inModule])