wrapText <- function(x, len) {
    sapply(x, function(y) paste(strwrap(y, len), collapse = "\n"), USE.NAMES = FALSE)
}
Great.enrich <- function(gr, bg, title, species = "hg38", fdr.threshold = 0.05, Type = "GO Biological Process", n_terms = 20){
    require(rGREAT)
    require(tidyverse)
    require(cowplot)
    theme_set(theme_cowplot())
    set.seed(101)
    job <- submitGreatJob(gr = gr, bg = bg, species = species)
    tb <- getEnrichmentTables(job)
    region.graph <- plotRegionGeneAssociationGraphs(job)

    # "GO Molecular Function" "GO Biological Process" "GO Cellular Component"
    sig.tb <- lapply(names(tb), function(x){
        enr <- tb[[x]]
        enr$Term <- rep(x, nrow(enr))
        sig <- which(enr$Hyper_Adjp_BH < fdr.threshold)
        if(length(sig>0)){
            return(enr[sig,])
        }else{
            return(data.frame())
        }
    })
    rGREAT.res <- Reduce(rbind, sig.tb)
    if(nrow(rGREAT.res) > 0){
        if(!is.null(Type)){
            rGREAT.res <- rGREAT.res[which(rGREAT.res$Term == Type),]
        }

        if(nrow(rGREAT.res) > 0){
            if(nrow(rGREAT.res) >= n_terms){
            rGREAT.res <- rGREAT.res[1:n_terms,]
            }
            rGREAT.res$wrap <- wrapText(rGREAT.res$name, 45)
            p <- ggplot(rGREAT.res, aes(y = log(Hyper_Fold_Enrichment), x = reorder(wrap, log(Hyper_Fold_Enrichment)))) +
                geom_bar(stat='identity', aes(color = name, fill = name)) +
                xlab('') + theme(legend.position="none") + ggtitle(title) + coord_flip()
            print(p)
        }
    }
    return(rGREAT.res)
}