一文掌握单基因GSEA富集分析

    https://www.gsea-msigdb.org/gsea/index.jsp

    https://www.gsea-msigdb.org/gsea/downloads.jsp

    https://www.gsea-msigdb.org/gsea/downloads.jsp

    https://www.gsea-msigdb.org/gsea/msigdb/human/collections.jsp'

    ###'@导入所需的R包
    library(org.Hs.eg.db)
    library(stringr)
    library(BiocGenerics)
    library(clusterProfiler)
    library(enrichplot)
    library(future)
    library(future.apply)

    exprSet <- read.csv("Input_ExpData.csv",row.names = 1)
    exprSet[1:10,1:5]

    > exprSet[1:10,1:5]
           sample01   sample02  sample03    sample04    sample05
    VWA1     2.26049224 4.08491679 3.4792744  1.99819783  2.36426617
    ATAD3C   0.17494291 0.23927633 0.1331588  0.26317521  0.11831695
    ATAD3B   0.45914335 0.62300446 0.7261913  0.23297821  0.33270716
    ATAD3A   1.21857160 1.53510784 1.3646565  1.50753047  1.03695410
    TMEM240  0.08058546 0.27554979 0.0000000  0.00000000  0.05109502
    SSU72    6.95506819 7.85755571 8.3072225 10.27444893 11.24621189
    C1orf233 0.08617899 0.07366901 0.1341729  0.02700902  0.00000000
    MIB2     0.96564412 0.91907833 1.0022802  0.56159762  0.63120040
    MMP23B   0.00000000 0.00000000 0.0000000  0.00000000  0.00000000
    CDK11B   0.70893552 0.92686092 1.3526327  1.18934155  1.29561449

    ##'@计算相关性
    ##'@定义函数:获得目标基因与其他基因之间的P值和Cor值
    batch_cor <- function(gene){
      y <- as.numeric(exprSet[gene,])
      rownames <- rownames(exprSet)
      do.call(rbind,future_lapply(rownames, function(x){
    dd  <- cor.test(as.numeric(exprSet[x,]),y,type="spearman")
    data.frame(gene=gene,mRNAs=x,cor=dd$estimate,p.value=dd$p.value )
      }))
    }

    dd <- batch_cor("HES5")

    gene <- dd$mRNAs

    gene = bitr(gene, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db")
    gene_df <- data.frame(cor=dd$cor,SYMBOL = dd$mRNAs)

    gene <- dplyr::distinct(gene,SYMBOL,.keep_all=TRUE)
    geneList <- gene_df$cor
    names(geneList)=gene_df$SYMBOL
    geneList=sort(geneList,decreasing = T)
    dotplot_internal <- function(object, x = "GeneRatio", color = "pvalue",
                             showCategory=10, size=NULL, split = NULL,
                             font.size=12, title = "", orderBy="x", decreasing=TRUE) {
      
      colorBy <- match.arg(color, c("pvalue", "p.adjust", "qvalue"))
      if (x == "geneRatio" || x == "GeneRatio") {
    x <- "GeneRatio"
    if (is.null(size))
      size <- "Count"
      } else if (x == "count" || x == "Count") {
    x <- "Count"
    if (is.null(size))
      size <- "GeneRatio"
      } else if (is(x, "formula")) {
    x <- as.character(x)[2]
    if (is.null(size))
      size <- "Count"
      } else {
    ## message("invalid x, setting to 'GeneRatio' by default")
    ## x <- "GeneRatio"
    ## size <- "Count"
    if (is.null(size))
      size  <- "Count"
      }
      df <- fortify(object, showCategory = showCategory, split=split)
      ## already parsed in fortify
      ## df$GeneRatio <- parse_ratio(df$GeneRatio)
      
      if (orderBy !=  'x' && !orderBy %in% colnames(df)) {
    message('wrong orderBy parameter; set to default `orderBy = "x"`')
    orderBy <- "x"
      }
      
      if (orderBy == "x") {
    df <- dplyr::mutate(df, x = eval(parse(text=x)))
      }
      
      idx <- order(df[[orderBy]], decreasing = decreasing)
      df$Description <- factor(df$Description, levels=rev(unique(df$Description[idx])))
      ggplot(df, aes_string(x=x, y="Description", size=size, color=colorBy)) +
    geom_point() +
    scale_color_continuous(low="red", high="blue", name = color, guide=guide_colorbar(reverse=TRUE)) +
    ## scale_color_gradientn(name = color, colors=sig_palette, guide=guide_colorbar(reverse=TRUE)) +
    ylab(NULL) + ggtitle(title) + DOSE::theme_dose(font.size) + scale_size(range=c(3, 8))
      
    }

    ##'@GSEA GO 富集
    ##'读取go.gmt文件
    GOgmt<-read.gmt("c5.go.v7.2.symbols.gmt")
    ##'@富集分析
    GO <-GSEA(geneList,TERM2GENE = GOgmt)

    dotplot_internal(GO)

    dotplot(GO,color="pvalue")

    dotplot(GO,split=".sign")+facet_grid(~.sign)

    gseaplot2(GO,1:10,color="red",pvalue_table = F)

    KEGGgmt<-read.gmt("c2.cp.kegg.v7.2.symbols.gmt")

    KEGG <-GSEA(geneList,TERM2GENE = KEGGgmt)

    dotplot_internal(KEGG)

    dotplot(KEGG,color="pvalue")

    dotplot(KEGG,split=".sign")+facet_grid(~.sign)

    gseaplot2(KEGG,1:10,color="red",pvalue_table = F)