david_enrichment_util.R 2.02 KB
Newer Older
1 2 3 4 5 6 7 8

########################################################################################################################
### enrichment analysis


## http://www.bioconductor.org/packages/release/bioc/vignettes/RDAVIDWebService/inst/doc/RDavidWS-vignette.pdf
## e.g. getClusterReport --> plot2D

9
## complete list can be fetched using  getAllAnnotationCategoryNames(david)
10 11 12 13 14 15 16 17 18 19 20 21
DEF_DAVID_ONTOLOGIES=ontologies=c("GOTERM_CC_FAT", "GOTERM_MF_FAT", "GOTERM_BP_FAT", "PANTHER_PATHWAY", "PANTHER_FAMILY", "PANTHER_PATHWAY", "KEGG_PATHWAY", "REACTOME_PATHWAY")

davidAnnotationChart <- function( someGenes, ontologies=DEF_DAVID_ONTOLOGIES ){

    require_auto(RDAVIDWebService) ## just works if installed on non-network-drive (e.g. /tmp/)

    ## expexted to have a column with gene_id
    #    echo("processing list with", length(someGenes), "genes")
    #    someGenes <- degs$ensembl_gene_id[1:100]


    if(length(someGenes)>1500){
22 23
        # stop("it's up to you to reduce your input to 1500 genes")
        warning("trimming david query list to 1500 entries")
24 25 26
        someGenes <- sample(someGenes) %>% head(1500)
    }

27
    david<-DAVIDWebService$new(email="bioinformatics@mpi-cbg.de")
28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55

    #    ## list all ontologies
    #    getAllAnnotationCategoryNames(david)


    #    getTimeOut(david)
    setTimeOut(david, 80000) ## http://www.bioconductor.org/packages/release/bioc/vignettes/RDAVIDWebService/inst/doc/RDavidWS-vignette.pdf

    result<-addList(david, someGenes, idType="ENSEMBL_GENE_ID", listName=paste0("list_", sample(10000)[1]), listType="Gene")

    david %>% setAnnotationCategories(ontologies)

    annoChart <-getFunctionalAnnotationChart(david)

    #    clusterReport <-getClusterReport(david)

    unloadNamespace('RDAVIDWebService')

    ## remove gene colum
    #    browser()
    annoChart <- as.data.frame(unclass(annoChart))

    # http://stackoverflow.com/questions/25271856/cannot-coerce-class-typeof-is-double-to-a-data-frame
    #    if(nrow(annoChart) >0) annoChart <-  annoChart %>%  dplyr::select(select=-Genes)

    return(annoChart)
}