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
#...3 enriched terms found
'data.frame': 3 obs. of 12 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 : num 0.00517 0.0163 0.04528
$ geneID : chr "332/6279/7153/4312/6280/6241/6790/983/2146" "6373/3627/4283/3002" "3627/4312/4283/4321"
$ Count : int 9 4 4
#...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 : num 0.00517 0.0163 0.04528
$ geneID : chr "332/6279/7153/4312/6280/6241/6790/983/2146" "6373/3627/4283/3002" "3627/4312/4283/4321"
$ 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 0.00516878 9
DOID:2799 21.665782 8.924992 3.045959e-05 0.01629588 0.01629588 4
DOID:14004 15.224604 7.333375 1.269586e-04 0.04528190 0.04528190 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 : num 0.00517 0.0163 0.04528
$ geneID : chr "332/6279/7153/4312/6280/6241/6790/983/2146" "6373/3627/4283/3002" "3627/4312/4283/4321"
$ 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')16.7 Alternative filtering approaches
While this chapter focuses on using dplyr verbs for manipulating enrichment results, users can also employ base R operations for filtering. For detailed information on using [, [[, and $ operators with the asis parameter to preserve enrichment result objects, please refer to the Data frame interface section in the Useful utilities chapter.
The dplyr approach provides a more readable and chainable syntax for complex filtering operations, while base R operators offer familiarity for users accustomed to traditional data frame manipulation.