#### 使用clusterProfiler进行通路富集分析

#' @param gene 字符向量，感兴趣的基因集合, 如"FUT7"         "AIRE"         "IFI35"

#' @param geneType 字符型向量，指明输入的geneList中name属性类型，即基因id的类型

#' @param db.type 字符型向量，指明GSEA针对的数据库

#' @param saveDir gsea输入目录

#' @param pAdjustMethod, pvalueCutoff p值矫正方法和阈值

#' @param MsigDB.category 指定特定的MsigDB中的数据库进行GSEA分析

#' MsigDB.pathway, msigdbr(species = "Homo sapiens", category = i)， 指定i为H, C1, C2, C3, C4, C5, C6, C7, C8

# 通过gs_cat gs_subcat，来选择库

#' @param GO.ont 指定GSEA分析GO数据库时使用的条目类型，支持"BP", "MF", and "CC"，以及"ALL"

#' 功能富集分析的原理是基于hypergeometric test

#' clusterProfiler主要基于EntrezID做分析

#' 分析过程中的universe，即背景基因的数目，是跟你选择的db中的subcategory有关系的，而不是固定不变的

cluterProfiler.enricher <- function(gene, geneType = c("ENTREZID", "SYMBOL", "ENSEMBL"),

                                    db.type = c("MsigDB", "GO", "KEGG"),

                                    saveDir, title,

                                    pAdjustMethod = "BH",

                                    qvalueCutoff = 0.2,

                                    pvalueCutoff = 0.05,

                                    MsigDB.category = list(H = c("All"), C2 = c("BIOCARTA", "KEGG", "REACTOME"), C5 = c("BP")),

                                    GO.ont = "BP", simplify = T,

                                    minGSSize = 10, maxGSSize = 500,

                                    showCategory = 20

                                    ){

  require(enrichplot)

  require(clusterProfiler)

  require(tidyverse)

  geneType <- match.arg(geneType)

  db.type <- match.arg(db.type)  

  # 如果基因名字类型不是SYMBOL，则需要进行ID转换

  if(geneType == "ENTREZID"){

    cat("不需要转换\n")

    geneList <- gene

  }

  if(geneType %in% c("SYMBOL", "ENSEMBL")){

    genes <- IDConvert(genes = gene, fromType = geneType, toType = "ENTREZID")

    cat("转换前基因数目：", length(gene), "\n")

    index <- match(gene, genes[,1])

    geneList <- as.character(genes[na.omit(index),2])

    cat("转换后基因数目：", length(geneList), "\n")    

  }

  ## 第一种富集策略：MsigDB基因集合

  if(db.type == "MsigDB"){

      library(msigdbr)

      # 通路提取

      MsigDB.pathway <- tibble(gs_name = NULL, entrez_gene = NULL)

      # 将MsigDB所有基因作为背景基因

      #MsigDB.gene <- unique(msigdbr(species = "Homo sapiens") %>% dplyr::select(gs_name,entrez_gene))

      if(class(MsigDB.category) == "character"){

        MsigDB.pathway <- msigdbr(species = "Homo sapiens") %>% dplyr::select(gs_name, entrez_gene)

      }else{

        for(i in names(MsigDB.category)){

          if(MsigDB.category[[i]][1] == "All"){

            MsigDB.pathway <- rbind(MsigDB.pathway, msigdbr(species = "Homo sapiens", category = i) %>% dplyr::select(gs_name, entrez_gene))

          }else{

            for(j in MsigDB.category[[i]]){

              MsigDB.pathway <- rbind(MsigDB.pathway, msigdbr(species = "Homo sapiens", category = i, subcategory = j) %>% dplyr::select(gs_name, entrez_gene))

            }

          }

        }

      }

      em.res <- enricher(geneList, TERM2GENE = MsigDB.pathway, qvalueCutoff = qvalueCutoff, pvalueCutoff = pvalueCutoff, pAdjustMethod = pAdjustMethod, minGSSize = minGSSize, maxGSSize = maxGSSize)

  }

  if(db.type == "GO"){

    em.res <- enrichGO(gene = geneList,

                    OrgDb = org.Hs.eg.db,

                    ont = GO.ont,minGSSize = minGSSize, maxGSSize = maxGSSize,

                    pAdjustMethod = pAdjustMethod,

                    pvalueCutoff  = pvalueCutoff,

                    qvalueCutoff  = qvalueCutoff)

    #基于simplify method, we can remove redundant terms

    if(simplify){

      em.res <- simplify(em.res, cutoff = 0.7, by="p.adjust", select_fun = min)

    }

  }

  if(db.type == "KEGG"){

    em.res <- enrichKEGG(gene = geneList, organism = "hsa", keyType = "kegg",

                       pAdjustMethod = pAdjustMethod, minGSSize = minGSSize, maxGSSize = maxGSSize,

                       pvalueCutoff  = pvalueCutoff,

                       qvalueCutoff  = qvalueCutoff)

  }

  #转换为geneSymbol

  em.res.geneSymbol <- setReadable(em.res, org.Hs.eg.db, keyType = "ENTREZID")

  # if(nrow(em.res@result) > 0){

  #   write.csv(em.res.geneSymbol, file = file.path(saveDir, paste0(title, "_", "enricherResult.csv")))

  #   tryCatch({

  #     p1 <- barplot(em.res.geneSymbol, showCategory = showCategory, main = title) + theme(axis.text.x = element_text(size=4), axis.text.y = element_text(size=4))

  #     print(p1)

  #     p2 <- dotplot(em.res.geneSymbol, showCategory = showCategory, title = title) + theme(axis.text.x = element_text(size=4), axis.text.y = element_text(size=4))

  #     print(p2)

  #     p3 <- heatplot(em.res.geneSymbol, showCategory = showCategory) + theme(axis.text.x = element_text(size=4), axis.text.y = element_text(size=4))

  #     print(p3)

  #   }, warning = function(w){

  #     # 这里是出现warning状态时，应该怎么做，可以用print打印出来，可以执行其它命令

  #     cat("Warning in plot!\n")

  #   }, error = function(e){

  #     # 这里时出现Error状态时，应该怎么做，可以用print打印出来，也可以执行其它命令

  #     cat("Error in plot!\n")

  #   },finally = {

  #     # 这是运行正常时，应该怎么做，可以用print打印出来，也可以执行其它命令

  #     cat("Plot complete!\n")

  #   })

  # }

  return(em.res = list(em.res.entrezid = em.res, em.res.genesymbol = em.res.geneSymbol))

}

##### ID转换

#author: Zhilin long

#time: 2021.1.19

#' @author Zhilin Long

#' @param gene 字符型向量，需要转换的ID向量

#' @param fromType 字符向量，初始的ID类型

#' @param toType 字符向量，需要转换的ID类型

#' 

#' @return data.frame 已经去除了重复和NA，可直接使用

IDConvert <- function(genes, 

                      method = c("clusterProfiler", "org.Hs.eg.db"),

                      fromType = c("ENTREZID", "SYMBOL", "ENSEMBL"), 

                      toType = c("ENTREZID", "SYMBOL", "ENSEMBL")

                      ){

    require(org.Hs.eg.db)

    require(clusterProfiler)

    #规范参数输入

    method <- match.arg(method)

    fromType <- match.arg(fromType)

    toType <- match.arg(toType)

    if(method == "clusterProfiler"){

        gene.df <- bitr(genes, fromType = fromType, toType = toType, OrgDb = "org.Hs.eg.db")

        #默认会过滤没有map到的id，即条目会减少

    }else{

        gene.df <- mapIds(org.Hs.eg.db, 

                          keys = genes, 

                          column = toType, 

                          keytype = fromType,

                          multiVals="first")

        gene.df <- data.frame(ENSEMBL = names(gene.df), SYMBOL = gene.df)

        colnames(gene.df) <- c(fromType, toType)

        #没有map，则为NA，即条目不变

        gene.df <- na.omit(gene.df) #去除NA的行

    }

    #去除重复, 

    gene.df <- gene.df[!duplicated(gene.df[,2]),]

    gene.df <- gene.df[!duplicated(gene.df[,1]),]

    return(gene.df)

}