R/pathway_analysis2.R

In LittleHeronCodes/Lazy2: Lazy Package

#' Gene list pathway analysis by hypergeometric test
#'
#' Run pathway analysis (hypergeometric test) for gene list. Wrapper function for hypergeoTestForGeneset.
#' @param glist List of query genes. (DEG list)
#' @param refgmt Reference gene set list. (Pathway set in list form)
#' @param tglist List of background gene space, same length as glist.
#' @param minGeneSet Minimum number of genes in gene set. Default 10.
#' @param ncore Number of cores to use.
#' @param ef.psc Pseudocount option for calculating enrichment factors. Default 0.
#' @return List of enricher Object: up, down
#' @export

Gen_enrichment <- function(glist, refgmt, tglist, minGeneSet = 10, ncore = 1, ef.psc = 0) {
    # .Deprecated("glist_enrichment")
    bgspace <- unique(unlist(refgmt))
    glist <- lapply(glist, function(gg) intersect(gg, bgspace))
    tglist <- lapply(tglist, function(gg) intersect(gg, bgspace))

    if (ncore <= 1) {
        enrobj <- lapply(names(glist), function(aid) {
			hypergeoTestForGeneset(glist[[aid]], refgmt, tglist[[aid]], minGeneSet, ef.psc = ef.psc)
			})
    }
    if (ncore > 1) {
        enrobj <- mclapply(names(glist), function(aid) {
			hypergeoTestForGeneset(glist[[aid]], refgmt, tglist[[aid]], minGeneSet, ef.psc = ef.psc)
		}, mc.cores = ncore)
    }
    names(enrobj) <- names(glist)
    return(enrobj)
}


#' Enrichment object list to matrix
#'
#' Turn enrichment object list (generated by ClusterProfiler enricher or Gen_enrichment) into matrix of needed value for drawing heatmap.
#' @param enrobj Enrichment result list.
#' @param val.col Column to output.
#' @param log Should the output be log10 transformed? (used if p value is input)
#' @return List of enricher Object: up, down
#' @importFrom dplyr rename
#' @importFrom reshape2 acast
#' @export


enrobj2Matrix <- function(enrobj, val.col = "pvalue", log = TRUE) {
    LS <- lapply(names(enrobj), function(set) {
        dff <- data.frame(enrobj[[set]])
        if ("Description" %in% names(dff)) dff <- dff %>% dplyr::rename("termID" = "ID", "ID" = "Description")
        dff$set <- set
        dff <- dff[order(dff$set), ]
        return(dff)
    })
    hmplot <- do.call(rbind, LS)

    # detect log values
    is_log <- FALSE
    if (any(quantile(hmplot[, val.col], na.rm = TRUE) > 1)) is_log <- TRUE
    if (log & is_log) cat("val.col seems to already be in log values.\n")
    log <- FALSE
    if (log & !is_log) {
        hmplot$logV <- -log10(hmplot[, val.col])
        plotMat <- reshape2::acast(hmplot, ID ~ set, value.var = "logV", fill = NA)
    }

    if (!log) plotMat <- reshape2::acast(hmplot, ID ~ set, value.var = val.col, fill = NA)
    plotMat <- plotMat[order(apply(plotMat, 1, sum, na.rm = TRUE), decreasing = TRUE), ]
    return(plotMat)
}