library(DOSE)
data(geneList)
de = names(geneList)[1:100]
x = enrichDO(de)16 dplyr verbs for manipulating enrichment result
16.1 filter
filter(x, p.adjust < .05, qvalue < 0.2)#
# over-representation test
#
#...@organism Homo sapiens
#...@ontology HDO
#...@keytype ENTREZID
#...@gene chr [1:100] "4312" "8318" "10874" "55143" "55388" "991" "6280" "2305" ...
#...pvalues adjusted by 'BH' with cutoff < 0.05
#...0 enriched terms found
#...Citation
Guangchuang Yu, Li-Gen Wang, Guang-Rong Yan, Qing-Yu He. DOSE: an R/Bioconductor package for Disease Ontology Semantic and Enrichment analysis. Bioinformatics. 2015, 31(4):608-609 16.2 arrange
mutate(x, geneRatio = parse_ratio(GeneRatio)) %>%
arrange(desc(geneRatio))#
# over-representation test
#
#...@organism Homo sapiens
#...@ontology HDO
#...@keytype ENTREZID
#...@gene chr [1:100] "4312" "8318" "10874" "55143" "55388" "991" "6280" "2305" ...
#...pvalues adjusted by 'BH' with cutoff < 0.05
#...3 enriched terms found
'data.frame': 3 obs. of 13 variables:
$ ID : chr "DOID:11054" "DOID:2799" "DOID:14004"
$ Description : chr "urinary bladder cancer" "bronchiolitis obliterans" "thoracic aortic aneurysm"
$ GeneRatio : chr "9/58" "4/58" "4/58"
$ BgRatio : chr "183/8168" "26/8168" "37/8168"
$ RichFactor : num 0.0492 0.1538 0.1081
$ FoldEnrichment: num 6.93 21.67 15.22
$ zScore : num 6.86 8.92 7.33
$ pvalue : num 4.83e-06 3.05e-05 1.27e-04
$ p.adjust : num 0.00517 0.0163 0.04528
$ qvalue : logi NA NA NA
$ geneID : chr "6790/6279/2146/4312/983/6280/332/6241/7153" "3002/6373/3627/4283" "4321/4312/3627/4283"
$ Count : int 9 4 4
$ geneRatio : num 0.155 0.069 0.069
#...Citation
Guangchuang Yu, Li-Gen Wang, Guang-Rong Yan, Qing-Yu He. DOSE: an R/Bioconductor package for Disease Ontology Semantic and Enrichment analysis. Bioinformatics. 2015, 31(4):608-609 16.3 select
select(x, -geneID) %>% head ID Description GeneRatio BgRatio RichFactor
DOID:11054 DOID:11054 urinary bladder cancer 9/58 183/8168 0.04918033
DOID:2799 DOID:2799 bronchiolitis obliterans 4/58 26/8168 0.15384615
DOID:14004 DOID:14004 thoracic aortic aneurysm 4/58 37/8168 0.10810811
FoldEnrichment zScore pvalue p.adjust qvalue Count
DOID:11054 6.925947 6.856156 4.830635e-06 0.00516878 NA 9
DOID:2799 21.665782 8.924992 3.045959e-05 0.01629588 NA 4
DOID:14004 15.224604 7.333375 1.269586e-04 0.04528190 NA 416.4 mutate
# k/M
y <- mutate(x, richFactor = Count / as.numeric(sub("/\\d+", "", BgRatio)))
y#
# over-representation test
#
#...@organism Homo sapiens
#...@ontology HDO
#...@keytype ENTREZID
#...@gene chr [1:100] "4312" "8318" "10874" "55143" "55388" "991" "6280" "2305" ...
#...pvalues adjusted by 'BH' with cutoff < 0.05
#...3 enriched terms found
'data.frame': 3 obs. of 13 variables:
$ ID : chr "DOID:11054" "DOID:2799" "DOID:14004"
$ Description : chr "urinary bladder cancer" "bronchiolitis obliterans" "thoracic aortic aneurysm"
$ GeneRatio : chr "9/58" "4/58" "4/58"
$ BgRatio : chr "183/8168" "26/8168" "37/8168"
$ RichFactor : num 0.0492 0.1538 0.1081
$ FoldEnrichment: num 6.93 21.67 15.22
$ zScore : num 6.86 8.92 7.33
$ pvalue : num 4.83e-06 3.05e-05 1.27e-04
$ p.adjust : num 0.00517 0.0163 0.04528
$ qvalue : logi NA NA NA
$ geneID : chr "6790/6279/2146/4312/983/6280/332/6241/7153" "3002/6373/3627/4283" "4321/4312/3627/4283"
$ Count : int 9 4 4
$ richFactor : num 0.0492 0.1538 0.1081
#...Citation
Guangchuang Yu, Li-Gen Wang, Guang-Rong Yan, Qing-Yu He. DOSE: an R/Bioconductor package for Disease Ontology Semantic and Enrichment analysis. Bioinformatics. 2015, 31(4):608-609 library(ggplot2)
library(forcats)
library(enrichplot)
ggplot(y, showCategory = 20,
aes(richFactor, fct_reorder(Description, richFactor))) +
geom_segment(aes(xend=0, yend = Description)) +
geom_point(aes(color=p.adjust, size = Count)) +
scale_color_viridis_c(guide=guide_colorbar(reverse=TRUE)) +
scale_size_continuous(range=c(2, 10)) +
theme_minimal() +
xlab("rich factor") +
ylab(NULL) +
ggtitle("Enriched Disease Ontology")
A very similar concept is Fold Enrichment, which is defined as the ratio of two proportions, (k/n) / (M/N). Using mutate to add the fold enrichment variable is also easy:
mutate(x, FoldEnrichment = parse_ratio(GeneRatio) / parse_ratio(BgRatio))Here, the calculation of rich factor and fold enrichment is only for demonstration purposes. The enrichplot package provides the dotplot function that can directly visualize these two values without adding them to the enrichment result.
16.5 slice
We can use slice to choose rows by their ordinal position in the enrichment result. Grouped result use the ordinal position with the group.
In the following example, a GSEA result of Reactome pathway was sorted by the absolute values of NES and the result was grouped by the sign of NES. We then extracted first 5 rows of each groups. The result was displayed in Figure 16.2.
library(ReactomePA)
x <- gsePathway(geneList)
y <- arrange(x, abs(NES)) %>%
group_by(sign(NES)) %>%
slice(1:5)
library(forcats)
library(ggplot2)
library(enrichplot)
ggplot(y, aes(NES, fct_reorder(Description, NES), fill=qvalue), showCategory=10) +
geom_col(orientation='y') +
scale_fill_continuous(low='red', high='blue', guide=guide_colorbar(reverse=TRUE)) +
theme_minimal() + ylab(NULL)
16.6 summarise
library(ggplot2)
pbar <- function(x) {
pi=seq(0, 1, length.out=11)
mutate(x, pp = cut(p.adjust, pi)) |>
group_by(pp) |>
summarise(cnt = n()) |>
ggplot(aes(pp, cnt)) + geom_col() +
theme_minimal() +
xlab("p value intervals") +
ylab("Frequency") +
ggtitle("p value distribution")
}
x <- enrichDO(de, pvalueCutoff=1, qvalueCutoff=1)
set.seed(2020-09-10)
random_genes <- sample(names(geneList), 100)
y <- enrichDO(random_genes, pvalueCutoff=1, qvalueCutoff=1)
p1 <- pbar(x)
p2 <- pbar(y)
aplot::plot_list(p1, p2, ncol=1, tag_levels = 'A')