用clusterProfiler做GSEA（一）

生信编程日常

发布于 2020-04-01 15:45:10

2.6K0

发布于 2020-04-01 15:45:10

文章被收录于专栏：生物信息学、python、R、linux生物信息学、python、R、linux

image.png

GSEA的介绍：https://www.omicsclass.com/article/230 GSEA有相应的软件，其实clusterProfiler除了做go term 富集，也可以做GSEA。首先介绍GSEA需要的文件： 1.GSEA输入的geneList要求是数值型向量，可以是fold change，或者logFC，数值型向量的名字是基因ID，数字从高到低排序,如：

image.png

2.其次是需要富集的go term 及基因，形式为两列一列pathway，一列gene,如：

image.png

如果是需要Molecular Signature Database (MSigDB) 的数据的话，可以安装msigdf 包。

devtools::install_github("ToledoEM/msigdf")
library(msigdf)
library(dplyr)
c2 <- msigdf.human %>% 
    filter(category_code == "c2") %>% select(geneset, symbol) %>% as.data.frame

如果是需要自己从go term里选择合适的数据，如大鼠的wnt 信号通路，可以

source('GetGoTerm.R')#clusterProfiler里有这个GetGoTerm.R
GO_DATA_NRVC <- get_GO_data("org.Rn.eg.db", "ALL", "SYMBOL") 　
Wnt_NRVCGO<-names(GO_DATA$PATHID2NAME[grep("Wnt", GO_DATA$PATHID2NAME)])
write.csv(GO_DATA$PATHID2NAME[Wnt_NRVCGO],"NRVC_Wnt_NRVC_related.csv")
Wnt_NRVCgo<-unlist(GO_DATA$PATHID2EXTID["GO:0016055"])
length(Wnt_NRVCgo)
#需要基因id转换的话，不转化id这一行可以忽略
#Wnt_NRVC<-bitr(Wnt_NRVC, fromType="SYMBOL", toType=c("ENTREZID"), OrgDb="org.Hs.eg.db"); head(Wnt_NRVC)
Wnt_NRVCgo<-cbind(rep("Wnt signaling pathway",368),as.data.frame(Wnt_NRVCgo))
colnames(Wnt_NRVCgo)<-c("ont","gene")
head(Wnt_NRVCgo)

image.png

自己制作的go term就做好了

然后进行GSEA富集：

library(dplyr)
#NRVC_NM_tTag_count为edgeR输出的data.frame,其他来源均可，只要包含Fold change和基因名即可
geneList_NRVC <-dplyr::select(NRVC_tTag_count, Row.names, logFC) 
colnames(geneList_NRVC)[1]<-c("SYMBOL")
geneList_NRVC.sort <- arrange(geneList_NRVC, desc(logFC)); head(geneList_NRVC.sort)
#按FC降序
geneList_NRVC<-geneList_NRVC.sort$logFC
names(geneList_NRVC)<-geneList_NRVC.sort$SYMBOL

#GSEA富集Wnt信号通路
gsea_Wnt_NRVC <- GSEA(geneList_NRVC, TERM2GENE = Wnt_NRVCgo, verbose=FALSE, pvalueCutoff = 0.05); head(gsea_Wnt_NRVC)
library(DOSE)
DOSE::gseaplot(gsea_Wnt_NRVC, 1)

image.png

#GSEA富集人的c2通路
c2 <- msigdf.human %>% 
    filter(category_code == "c2") %>% dplyr::select(geneset, symbol) %>% as.data.frame
head(c2)
colnames(c2)<-c("ont","gene")
head(geneList_human)
gsea_c2_human <- GSEA(geneList_human, TERM2GENE = c2, verbose=FALSE, pvalueCutoff = 1)
library(DOSE)
DOSE::gseaplot(gsea_c2_human, 1)

image.png

image.png

附GetGoTerm.R的代码

library(DOSE)
library(GOSemSim)
library(clusterProfiler)
library(org.Hs.eg.db)
library(org.Mm.eg.db)
library(org.Rn.eg.db)
library(dplyr)
library(GO.db)
#
get_GO_data <- function(OrgDb, ont, keytype) {
  GO_Env <- get_GO_Env()
  use_cached <- FALSE
  
  if (exists("organism", envir=GO_Env, inherits=FALSE) &&
      exists("keytype", envir=GO_Env, inherits=FALSE)) {
    
    org <- get("organism", envir=GO_Env)
    kt <- get("keytype", envir=GO_Env)
    
    if (org == DOSE:::get_organism(OrgDb) &&
        keytype == kt &&
        exists("goAnno", envir=GO_Env, inherits=FALSE)) {
      ## https://github.com/GuangchuangYu/clusterProfiler/issues/182
      ## && exists("GO2TERM", envir=GO_Env, inherits=FALSE)){
      
      use_cached <- TRUE
    }
  }
  
  if (use_cached) {
    goAnno <- get("goAnno", envir=GO_Env)
  } else {
    OrgDb <- GOSemSim:::load_OrgDb(OrgDb)
    kt <- keytypes(OrgDb)
    if (! keytype %in% kt) {
      stop("keytype is not supported...")
    }
    
    kk <- keys(OrgDb, keytype=keytype)
    goAnno <- suppressMessages(
      AnnotationDbi::select(OrgDb, keys=kk, keytype=keytype,
             columns=c("GOALL", "ONTOLOGYALL")))
    
    goAnno <- unique(goAnno[!is.na(goAnno$GOALL), ])
    
    assign("goAnno", goAnno, envir=GO_Env)
    assign("keytype", keytype, envir=GO_Env)
    assign("organism", DOSE:::get_organism(OrgDb), envir=GO_Env)
  }
  
  if (ont == "ALL") {
    GO2GENE <- unique(goAnno[, c(2,1)])
  } else {
    GO2GENE <- unique(goAnno[goAnno$ONTOLOGYALL == ont, c(2,1)])
  }
  
  GO_DATA <- DOSE:::build_Anno(GO2GENE, get_GO2TERM_table())
  
  goOnt.df <- goAnno[, c("GOALL", "ONTOLOGYALL")] %>% unique
  goOnt <- goOnt.df[,2]
  names(goOnt) <- goOnt.df[,1]
  assign("GO2ONT", goOnt, envir=GO_DATA)
  return(GO_DATA)
}

get_GO_Env <- function () {
  if (!exists(".GO_clusterProfiler_Env", envir = .GlobalEnv)) {
    pos <- 1
    envir <- as.environment(pos)
    assign(".GO_clusterProfiler_Env", new.env(), envir=envir)
  }
  get(".GO_clusterProfiler_Env", envir = .GlobalEnv)
}

get_GO2TERM_table <- function() {
  GOTERM.df <- get_GOTERM()
  GOTERM.df[, c("go_id", "Term")] %>% unique
}

get_GOTERM <- function() {
  pos <- 1
  envir <- as.environment(pos)
  if (!exists(".GOTERM_Env", envir=envir)) {
    assign(".GOTERM_Env", new.env(), envir)
  }
  GOTERM_Env <- get(".GOTERM_Env", envir = envir)
  if (exists("GOTERM.df", envir = GOTERM_Env)) {
    GOTERM.df <- get("GOTERM.df", envir=GOTERM_Env)
  } else {
    GOTERM.df <- toTable(GOTERM)
    assign("GOTERM.df", GOTERM.df, envir = GOTERM_Env)
  }
  return(GOTERM.df)
}

参考：https://guangchuangyu.github.io/cn/2018/11/msigdf_clusterprofiler/

欢迎关注～

本文参与腾讯云自媒体分享计划，分享自作者个人站点/博客。

如有侵权请联系 cloudcommunity@tencent.com 删除