TCGA

再做GO+KEGG+GSEA

2021-01-27  本文已影响0人  晓颖_9b6f

1、再做GO

代码如下

DEG1 <- DEG_symbolid2
gene_up <- DEG1[DEG1$g == "up", "ENTREZID"] 
gene_down <- DEG1[DEG1$g == "down","ENTREZID" ]
gene_diff <- c(gene_up,gene_down)
gene_all <- as.character(DEG1[,"ENTREZID"])
geneList <- DEG1$log2FoldChange
names(geneList)=DEG1$ENTREZID
geneList=sort(geneList,decreasing = T)

if(F){
  go_enrich_results <- lapply( g_list , function(gene) {
    lapply( c('BP','MF','CC') , function(ont) {
      cat(paste('Now process ',ont ))
      ego <- enrichGO(gene          = gene,
                      universe      = gene_all,
                      OrgDb         = org.Hs.eg.db,
                      ont           = ont ,
                      pAdjustMethod = "BH",
                      pvalueCutoff  = 0.99,
                      qvalueCutoff  = 0.99,
                      readable      = TRUE)
      
      print( head(ego) )
      return(ego)
    })
  })
  save(go_enrich_results,file = 'go_enrich_results.Rdata')  
}


load(file = 'go_enrich_results.Rdata')
n1= c('gene_up','gene_down','gene_diff')
n2= c('BP','MF','CC') 
for (i in 1:3){
  for (j in 1:3){
    fn=paste0('dotplot_',n1[i],'_',n2[j],'.png')
    cat(paste0(fn,'\n'))
    png(fn,res=150,width = 1080)
    print( dotplot(go_enrich_results[[i]][[j]] ))
    dev.off()
  }
}
文章里面是做了生物过程和细胞 dotplot_gene_down_BP.png

的两种分析。


dotplot_gene_down_CC.png

2、再做KEGG分析

代码如下:

kk1_gse <- gseKEGG(geneList     = geneList,
                   organism     = 'hsa',
                   nPerm        = 1000,
                   minGSSize    = 120,
                   pvalueCutoff = 0.9,
                   verbose      = FALSE)
browseKEGG(kk1_gse, "hsa04662")     #B cell signaling
browseKEGG(kk1_gse, "hsa04062")    #chemokine signal transduction
browseKEGG(kk1_gse, "hsa04666")    #FcγR-mediated phagocytosis

input:但是没有红蓝色的标点。


image.png

所以我直接倒出文件,在网站上面做了。
https://www.genome.jp/kegg/tool/map_pathway2.html

library(openxlsx)
write.xlsx(DEGs,file = "DEGS")

input:


image.png image.png image.png

3、GSEA分析

代码:

kk_gse <- gseKEGG(geneList     = geneList,
                  organism     = 'hsa',
                  nPerm        = 1000,
                  minGSSize    = 120,
                  pvalueCutoff = 0.9,
                  verbose      = FALSE)
head(kk_gse)[,1:6]
gseaplot(kk_gse, geneSetID = rownames(kk_gse[1,]))

input:


image.png

想要做KEGG之后的热图分析,但是找富集到通路的gene-list,不知道怎么找。

上一篇下一篇

猜你喜欢

热点阅读