Skip to content

YuLab-SMU/Background-bias-of-functional-enrichment-analysis

Folders and files

NameName
Last commit message
Last commit date

Latest commit

 

History

5 Commits
 
 
 
 
 
 
 
 
 
 

Repository files navigation

Load packages and example data

library(ggplot2)
library(clusterProfiler)
library(org.Hs.eg.db)

allgenes <- keys(org.Hs.eg.db)
data(geneList, package="DOSE")
de <- names(geneList)[1:200]

Functional enrichment analysis

options(enrichment_force_universe=TRUE)
bp_allgene <- enrichGO(de, OrgDb = 'org.Hs.eg.db', ont="BP", universe=allgenes)
bp_bg1 <- enrichGO(de, OrgDb = 'org.Hs.eg.db', ont="BP", universe = names(geneList))

bp_allgene_kegg <- enrichKEGG(de, universe=allgenes)
bp_bg1_kegg <- enrichKEGG(de, universe = names(geneList))

options(enrichment_force_universe=FALSE)
bp_allgene2 <- enrichGO(de, OrgDb = 'org.Hs.eg.db', ont="BP", universe=allgenes)
bp_bg2 <- enrichGO(de, OrgDb = 'org.Hs.eg.db', ont="BP", universe = names(geneList))
bp_allgene2_kegg <- enrichKEGG(de, universe = allgenes)
bp_bg2_kegg <- enrichKEGG(de, universe = names(geneList))

Visualization

myplot <- function(res1, res2) {
    x <- merge(res1@result, res2@result, by="ID")
    ggplot(x, aes(pvalue.x, pvalue.y)) + 
        geom_point(color='steelblue') + 
        geom_smooth(method='lm', color='firebrick') +
        geom_abline(slope=1, intercept=0, linetype='dashed') +
        xlim(0, 1) + ylim(0, 1) + 
        DOSE::theme_dose()
}

p1 <- myplot(bp_allgene, bp_bg1) + 
    xlab('All genes in the genome as universe') +
    ylab("All expressed genes as universe")

p2 <- myplot(bp_bg2, bp_bg1) +
    xlab("All expressed genes as universe\n(intersected with gene sets)") +
    ylab("All expressed genes as universe")

p3 <- myplot(bp_allgene2, bp_bg2) + 
    xlab('All genes in the genome as universe\n(intersected with gene sets)') +
    ylab("All expressed genes as universe\n(intersected with gene sets)")

p4 <- myplot(bp_allgene_kegg, bp_bg1_kegg) + 
    xlab('All genes in the genome as universe') +
    ylab("All expressed genes as universe")

p5 <- myplot(bp_bg2_kegg, bp_bg1_kegg) +
    xlab("All expressed genes as universe\n(intersected with gene sets)") +
    ylab("All expressed genes as universe")

p6 <- myplot(bp_allgene2_kegg, bp_bg2_kegg) + 
    xlab('All genes in the genome as universe\n(intersected with gene sets)') +
    ylab("All expressed genes as universe\n(intersected with gene sets)")

aplot::plot_list(p1, p2, p3, p4, p5, p6, ncol=3, tag_levels='A')

About

No description, website, or topics provided.

Resources

Stars

Watchers

Forks

Releases

No releases published

Packages

No packages published