Motivation: at the end of this chapter, the students will be able to address the following points:
Differential expression analysis is univariate - each gene is tested on its own. This probably doesn’t reflect the underlying biology - genes work in conjunction, not in isolation. One wouldn’t expect that the effect of a drug, or a mutation, … would lead to the perturbation of a single gene expression.
The univariate approach expects (or tests for) significant changes in single genes. Moderate effects in many (related) genes cannot, by definition be identified as statistically significant, even it such an effect may be biologically more plausible that one of a few large ones.
The goal of an enrichment analysis is to test for a group of related genes, called gene sets, and test whether the genes within these sets are enriched for differentially expression.
While nothing would stop a user to create their own gene sets, the sets that are generally used stem from community-maintained resources.
The Gene Ontology (Ashburner et al. 2000Ashburner, M, C A Ball, J A Blake, D Botstein, H Butler, J M Cherry, A P Davis, et al. 2000. “Gene Ontology: Tool for the Unification of Biology. The Gene Ontology Consortium.” Nat Genet 25 (1): 25–29. https://doi.org/10.1038/75556.) defines GO terms. These terms are based on a controlled vocabulary (to resuse the same words in a consistent way) and relations11 Example of relations between terms are, for example, is_a (mitosis is_a cell cycle phase) or part_of (mitosis part_of M phase of mitotic cell cycle). that define the directed links between terms. These relations define a hierarchy between GO terms: all term can be represented as a (direct acyclic) graph (DAG). There exist thus very general and very specific terms.
These terms are classed into three categories, called namespaces:
Here’s an example of GO term GO:0007155
, that describes cell adhesion.
A Term from the GO ontology: GO:0007155
Label: cell adhesion
The attachment of a cell, either to another cell or to an underlying
substrate such as the extracellular matrix, via cell adhesion
molecules.
If you look at the GO term GO:0007155
entry on the Gene Ontology
page, you can
find more details and, if you scroll down, example genes that are
annotated with that term.
The whole Gene Ontology is can be accessed in R with the GO.db package.
In the code chunk below, we query the GO.db
package through the
org.Hs.eg.db
interface. This org.db-type of packages for Homo
sapiens enables to perform various queries for human genes, such as
retrieving all gene symbols and ENTREZ identifiers (the columns
below) that are annotated with a GO term (the keys
below) of
interest.
The GO term of interest here is focal adhesion:
A cell-substrate junction that anchors the cell to the extracellular
matrix and that forms a point of termination of actin filaments. In
insects focal adhesion has also been referred to as hemi-adherens
junction (HAJ).
library("GO.db")
library("org.Hs.eg.db")
0005925 <- AnnotationDbi::select(org.Hs.eg.db,
GO_keys = "GO:0005925",
columns = c("ENTREZID", "SYMBOL"),
keytype = "GO") %>%
as_tibble %>%
filter(!duplicated(ENTREZID))
## 'select()' returned 1:many mapping between keys and columns
0005925 GO_
## # A tibble: 424 × 5
## GO EVIDENCE ONTOLOGY ENTREZID SYMBOL
## <chr> <chr> <chr> <chr> <chr>
## 1 GO:0005925 HDA CC 60 ACTB
## 2 GO:0005925 HDA CC 70 ACTC1
## 3 GO:0005925 ISS CC 71 ACTG1
## 4 GO:0005925 HDA CC 81 ACTN4
## 5 GO:0005925 HDA CC 87 ACTN1
## 6 GO:0005925 IMP CC 88 ACTN2
## 7 GO:0005925 IMP CC 89 ACTN3
## 8 GO:0005925 HDA CC 102 ADAM10
## 9 GO:0005925 HDA CC 118 ADD1
## 10 GO:0005925 HDA CC 214 ALCAM
## # ℹ 414 more rows
We have 424 genes matching this GO term. There are thus 424 genes in the GO_0005925 GO set.
► Question
Repeat the code above to extract the genes annotated with the
GO:0005813
term for the centrosome:
A structure comprised of a core structure (in most organisms, a pair
of centrioles) and peripheral material from which a
microtubule-based structure, such as a spindle apparatus, is
organized. Centrosomes occur close to the nucleus during interphase
in many eukaryotic cells, though in animal cells it changes
continually during the cell-division cycle.
► Solution
KEGG pathway is a collection of manually drawn and curated pathway maps representing current knowledge of the molecular interaction, reaction and relation networks.
The figure below shows the pathways for the cell cycle in humans.
The KEGGREST package provides a client interface to the KEGG server.
Alike KEGG patway, Reactome is a free, open-source, curated and peer-reviewed pathway database. The Bioconductor reactome.db package provides access to reactome maps and annotations within R.
MSigDB is a collection of annotated gene sets for use with GSEA software. The MSigDB gene sets are divided into 9 collections:
The msigdbr CRAN package provides the MSigDB gene sets in a standard R data frame with key-value pairs.
To illustrate enrichment analyses, we will use the DESeq2 results
stored in the res_tbl
variable, computed in the previous
chapter.
We will focus on the genes that have an adjusted p-value (those that have been tested) and that have unique ENTREZ gene identifiers.
res_tbl %>%
res_tbl <- filter(!is.na(ENTREZID),
!is.na(padj),
!duplicated(ENTREZID)) %>%
mutate(ENTREZID = as.character(ENTREZID))
res_tbl
## # A tibble: 14,006 × 9
## ENSEMBL gene ENTREZID baseMean log2FoldChange lfcSE stat pvalue
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ENSG00000101856 PGRMC1 10857 1209. -4.67 0.178 -26.2 8.24e-152
## 2 ENSG00000177494 ZBED2 79413 277. 6.25 0.351 17.8 6.25e- 71
## 3 ENSG00000136159 NUDT15 55270 630. 2.24 0.154 14.6 3.32e- 48
## 4 ENSG00000175315 CST6 1474 220. 4.29 0.300 14.3 2.20e- 46
## 5 ENSG00000128245 YWHAH 7533 2692. 2.05 0.144 14.2 5.64e- 46
## 6 ENSG00000113272 THG1L 54974 381. 2.33 0.171 13.6 3.45e- 42
## 7 ENSG00000213853 EMP2 2013 506. 2.53 0.191 13.2 8.45e- 40
## 8 ENSG00000115107 STEAP3 55240 193. 3.44 0.271 12.7 7.18e- 37
## 9 ENSG00000008513 ST3GA… 6482 1419. 3.00 0.240 12.5 6.39e- 36
## 10 ENSG00000105855 ITGB8 3696 6134. 2.18 0.178 12.2 2.26e- 34
## # ℹ 13,996 more rows
## # ℹ 1 more variable: padj <dbl>
An over representation analysis relies on the hypergeometric distribution. The hypergeometric distribution is a discrete probability distribution that describes the probability of \(k\) successes (random draws for which the object drawn has a specified feature) in \(n\) draws, without replacement, from a finite population of size \(N\) that contains exactly \(K\) objects with that feature, wherein each draw is either a success or a failure.
\[ P(X = k) = \frac{\binom{N}{k} \binom{N-K}{n-k}}{\binom{N}{n}} \]
where
The example used to describe the distribution is an urn contain \(N\) marbles of two colours. There are \(K\) green and and \(N-K\) red marbles in the urn. Drawing a green marble is defined as success, and a red marble failure. Using the formula above, we can compute the probability to draw \(k\) green marbles from the urn.
In the frame of an enrichment analysis (Rivals et al. 2007Rivals, I, L Personnaz, L Taing, and M C Potier. 2007. “Enrichment or Depletion of a Go Category Within a Class of Genes: Which Test?” Bioinformatics 23 (4): 401–7. https://doi.org/10.1093/bioinformatics/btl633.), we use the following formulation to calculate a probabiliy that we have more green marbles than we would expect by change, i.e. there to be an over representation of green marbles.
\[ p = 1 - \sum_{i = 0}^{k - 1} \frac{\binom{N}{k} \binom{N-K}{n-k}}{\binom{N}{n}} \]
To perform an over representation analysis, we thus need to define:
And fill out the following table and count the number of DE genes that are in the set of interest, the non-DE that are in the set, and the DE and non-DE genes that are not in the set:
GO | not_GO | |
---|---|---|
DE | n | p |
not_DE | m | q |
## DE and GO
length(intersect(res_tbl$ENTREZID[res_tbl$padj < 0.05],
n <-0005925$ENTREZID))
GO_
## not DE and GO
length(intersect(res_tbl$ENTREZID[res_tbl$padj >= 0.05],
m <-0005925$ENTREZID))
GO_
## DE and not GO
length(setdiff(res_tbl$ENTREZID[res_tbl$padj < 0.05],
p <-0005925$ENTREZID))
GO_
## not DE not not GO
length(setdiff(res_tbl$ENTREZID[res_tbl$padj >= 0.05],
q <-0005925$ENTREZID)) GO_
matrix(c(n, m, p, q), nrow = 2)
cont_mat <-rownames(cont_mat) <- c("DE", "not_DE")
colnames(cont_mat) <- c("GO", "not_GO")
cont_mat
## GO not_GO
## DE 159 4037
## not_DE 210 9600
We can now apply a Fisher’s exact (or hypergeometric test) that will test whether we can identify a statistically enrichment of DE genes in the GO category.
fisher.test(cont_mat, alternative = "greater")
##
## Fisher's Exact Test for Count Data
##
## data: cont_mat
## p-value = 4.253e-08
## alternative hypothesis: true odds ratio is greater than 1
## 95 percent confidence interval:
## 1.501133 Inf
## sample estimates:
## odds ratio
## 1.800321
Note that we could keep the default alternative = "two.sided"
to
test to over or under representation.
► Question
Repeat the ORA analysis above the GO:0005813 term.
This approach is straightfoward and very fast. Its major drawback however is that we need to define a cutoff to differentiate DE from non-DE genes. Setting this threshold might have a effect on the results.
► Question
Try setting different DE genes and check if, in the cases above, this has and effect on the GO terms of interest.
Gene set enrichment analysis refers to a broad family of tests. Here, we will define the principles based on (Subramanian et al. 2005Subramanian, A, P Tamayo, V K Mootha, S Mukherjee, B L Ebert, M A Gillette, A Paulovich, et al. 2005. “Gene Set Enrichment Analysis: A Knowledge-Based Approach for Interpreting Genome-Wide Expression Profiles.” Proc Natl Acad Sci U S A 102 (43): 15545–50. https://doi.org/10.1073/pnas.0506580102.), keeping in mind that the exact implementation will differ in different tools.
The major advantage of GSEA approaches is that they don’t rely on
defining DE genes. The first step is to order the genes of
interest based on the statistics used; here, we will use the
p-values. This is already the case for our res_tbl
table. We also
need to know which genes are in our set of interest.
$inGO <- res_tbl$ENTREZID %in% GO_0005925$ENTREZID
res_tbl::select(res_tbl, ENTREZID, padj, inGO) dplyr
## # A tibble: 14,006 × 3
## ENTREZID padj inGO
## <chr> <dbl> <lgl>
## 1 10857 1.36e-147 FALSE
## 2 79413 5.17e- 67 FALSE
## 3 55270 1.83e- 44 FALSE
## 4 1474 9.11e- 43 FALSE
## 5 7533 1.87e- 42 FALSE
## 6 54974 9.51e- 39 FALSE
## 7 2013 2.00e- 36 FALSE
## 8 55240 1.48e- 33 FALSE
## 9 6482 1.17e- 32 FALSE
## 10 3696 3.75e- 31 TRUE
## # ℹ 13,996 more rows
We are now going to compute a score by traversing to ordered gene list
and count a positive score when we encounter a gene in the gene set,
and a -1 when the gene is not in the gene set. The positive score,
computed by set_ratio()
below, is defined by \(\frac{n_{genes} - n_{genes~in~set}}{n_{genes~in~set}}\) so that the sum of all genes in
the set and those not in the set becomes zero.
function(tbl) {
set_ratio <-nrow(tbl) - sum(tbl$inGO))/sum(tbl$inGO)
( }
The cumulative sum of these scores along the ordered gene list becomes the GSEA path which is plotted below, and the maximun score obtained along this path is the GSEA score.
ifelse(res_tbl$inGO, set_ratio(res_tbl), -1)
gsea_path <- cumsum(gsea_path)
gsea_path <- max(gsea_path)
gsea_score <-plot(gsea_path, type = "l",
xlab = "Ordered genes (padj)",
ylab = "score",
main = paste0("GSEA score: ", round(gsea_score)))
abline(h = gsea_score, col = "steelblue", lty = "dotted")
To be able to compute a statistical significance, we need to compute a null distribution of GSEA scores, i.e. a distribution of scores reflecting the absence of any enrichment. This is done by permuting the samples in the original data, repeating the statistical analysis, reorder the genes according to the new statistics (we used the adjusted p-value above) and compute a new GSEA score.
dds
dds_tmp <-$Condition <- dds_tmp$Condition[sample(ncol(dds_tmp))]
dds_tmpcolData(dds)
## DataFrame with 6 rows and 4 columns
## Cell Type Condition sizeFactor
## <character> <character> <factor> <numeric>
## sample1 Cell1 Epithelial mock 0.726883
## sample2 Cell1 Epithelial mock 1.308664
## sample3 Cell1 Epithelial mock 1.077368
## sample4 Cell1 Epithelial KD 0.966025
## sample5 Cell1 Epithelial KD 1.042111
## sample6 Cell1 Epithelial KD 1.012664
colData(dds_tmp)
## DataFrame with 6 rows and 4 columns
## Cell Type Condition sizeFactor
## <character> <character> <factor> <numeric>
## sample1 Cell1 Epithelial KD 0.726883
## sample2 Cell1 Epithelial mock 1.308664
## sample3 Cell1 Epithelial KD 1.077368
## sample4 Cell1 Epithelial mock 0.966025
## sample5 Cell1 Epithelial mock 1.042111
## sample6 Cell1 Epithelial KD 1.012664
tibble(ENSEMBL = rownames(dds)) %>%
ensembl_to_geneName <- mutate(gene = AnnotationDbi::mapIds(org.Hs.eg.db, ENSEMBL, "SYMBOL",
"ENSEMBL")) %>%
mutate(ENTREZID = AnnotationDbi::mapIds(org.Hs.eg.db, ENSEMBL, "ENTREZID",
"ENSEMBL"))
DESeq(dds_tmp) %>%
res_tmp <- results(name = "Condition_KD_vs_mock") %>%
as_tibble(rownames = "ENSEMBL") %>%
left_join(ensembl_to_geneName) %>%
mutate(inGO = ENTREZID %in% GO_0005925$ENTREZID) %>%
filter(!is.na(ENTREZID),
!is.na(padj),
!duplicated(ENTREZID)) %>%
dplyr::select(ENSEMBL, ENTREZID, padj, inGO) %>%
arrange(padj)
res_tmp %>%
res_tmp <- mutate(score = ifelse(inGO, set_ratio(res_tmp), -1)) %>%
mutate(score = cumsum(score))
res_tmp
## # A tibble: 22,214 × 5
## ENSEMBL ENTREZID padj inGO score
## <chr> <chr> <dbl> <lgl> <dbl>
## 1 ENSG00000278267 102466751 1.00 FALSE -1
## 2 ENSG00000233750 100420257 1.00 FALSE -2
## 3 ENSG00000228463 728481 1.00 FALSE -3
## 4 ENSG00000236679 728517 1.00 FALSE -4
## 5 ENSG00000225972 100887749 1.00 FALSE -5
## 6 ENSG00000225630 100652939 1.00 FALSE -6
## 7 ENSG00000237973 107075141 1.00 FALSE -7
## 8 ENSG00000229344 107075310 1.00 FALSE -8
## 9 ENSG00000248527 106480796 1.00 FALSE -9
## 10 ENSG00000198744 107075270 1.00 FALSE -10
## # ℹ 22,204 more rows
max(res_tmp$score)
## [1] 431.4161
plot(res_tmp$score, type = "l")
This approach is very time consuming, given that the statistical tests for all the genes need to be recomputed at each permutation. In addition, for experiments with limited number of samples, the number of permutations would be limited.
We start by defining a function that will repeat the steps above and return a list containing the GSEA score and path.
##' A function that assigns the permuted Conditions,
##' runs DESeq and a GSEA analysis
##'
##' @param x a DESeq object.
##' @param perm `character()` of length `ncol(x)` indicating the
##' permutation to be tested.
##'
##' @return `numeric()` containing the GSEA path.
function(x, perm) {
gsea_perm <-## Set the Condition based on the permutation
$Condition[perm] <- "mock"
x$Condition[-perm] <- "KD"
x## Run DESeq2, extract and annotate results
suppressMessages(
DESeq(x) %>%
tbl <- results(name = "Condition_KD_vs_mock") %>%
as_tibble(rownames = "ENSEMBL") %>%
left_join(ensembl_to_geneName) %>%
mutate(inGO = ENTREZID %in% GO_0005925$ENTREZID) %>%
filter(!is.na(ENTREZID),
!is.na(padj),
!duplicated(ENTREZID)) %>%
dplyr::select(ENSEMBL, ENTREZID, padj, inGO) %>%
arrange(padj)
)## Compute GSEA results
tbl %>%
tbl <- mutate(score = if_else(inGO, set_ratio(tbl), -1)) %>%
mutate(score = cumsum(score))
return(tbl$score)
}
We will execute this function on all unique permutations, corresponding on the 10 first columns (including the actual design, in the first column). Note that columns 11 to 20 are simply the opposite of the first 10.
combn(6, 3)) (perms <-
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
## [1,] 1 1 1 1 1 1 1 1 1 1 2 2 2 2
## [2,] 2 2 2 2 3 3 3 4 4 5 3 3 3 4
## [3,] 3 4 5 6 4 5 6 5 6 6 4 5 6 5
## [,15] [,16] [,17] [,18] [,19] [,20]
## [1,] 2 2 3 3 3 4
## [2,] 4 5 4 4 5 5
## [3,] 6 6 5 6 6 6
We not apply our function on every permutation.
apply(perms[, 1:10], 2, function(p) gsea_perm(dds, p)) paths <-
And calculate the scores
sapply(paths, max) scores <-
If we compare the actual score to the permutation scores, we see that it is the largest one.
plot(density(scores))
rug(scores)
abline(v = scores[1])
An empirical p-value can be computer by dividing the number of null scores that are greater than the real score divided by the number of permutations (the number of null scores). In our case, given that no null scores are greater, the nominal p-value would be 0.
sum(scores > scores[1])/length(scores)
## [1] 0
Below, we illustrate all GSEA path and confirm that the actual score is indeed the largest one.
plot(paths[[1]], type = "l",
ylim = c(min(unlist(paths)), max(scores)),
col = "steelblue")
grid()
for (i in 2:length(paths))
lines(paths[[i]], type = "l", col = "#00000060")
## Warning in plot.xy(xy.coords(x, y), type = type, ...): semi-transparency is not
## supported on this device: reported only once per page
► Question
Given that we have 10 permutations, how many null scores greater that the actual GSEA score do we need for our results to become non significant (assuming we set an alpha of 0.05)?
► Solution
There exist various approaches to the GSEA analysis, that apply different permutation approaches to increate the number of possible permutations and reduce the running time.
► Question
Repeat the GSEA analysis above the GO:0005813 term.
While there is not need to define a list of differentially expressed genes for the GSEA analysis, the genes need to be ordered, which can be done in different ways and lead to different results.
clusterProfiler
package
I practice, the analyses presented above are executed using any of the very many packages that are available. Here, we will demonstrate clusterProfiler12 https://yulab-smu.top/clusterProfiler-book/, that itself relies on other packages to perform the GSEA-related computations.
library("clusterProfiler")
These dedicated packages will perform either of the two tests (or variations thereof) to all sets, adjust the computed p-values for multiplicity, and provide additional details and visualisations for further exploration of the results.
res_tbl$ENTREZID[res_tbl$padj < 0.05]
de_genes <-
enrichGO(gene = de_genes,
go_ora <-universe = res_tbl$ENTREZID,
OrgDb = org.Hs.eg.db,
ont = "CC",
readable = TRUE) %>%
as_tibble
go_ora
## # A tibble: 21 × 9
## ID Description GeneRatio BgRatio pvalue p.adjust qvalue geneID Count
## <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <chr> <int>
## 1 GO:0022… cytosolic … 70/3862 113/12… 4.79e-11 3.23e-8 3.04e-8 RPL27… 70
## 2 GO:0022… cytosolic … 40/3862 56/120… 1.52e- 9 5.12e-7 4.83e-7 RPL27… 40
## 3 GO:0044… ribosomal … 92/3862 185/12… 3.30e- 7 7.40e-5 6.96e-5 DAP3/… 92
## 4 GO:0005… ribosome 105/3862 226/12… 3.20e- 6 5.39e-4 5.08e-4 DAP3/… 105
## 5 GO:0005… focal adhe… 158/3862 368/12… 4.93e- 6 6.63e-4 6.24e-4 ITGB8… 158
## 6 GO:0030… cell-subst… 159/3862 374/12… 9.01e- 6 1.01e-3 9.52e-4 ITGB8… 159
## 7 GO:0015… large ribo… 57/3862 114/12… 4.48e- 5 4.31e-3 4.06e-3 MRPL1… 57
## 8 GO:0098… neuron to … 103/3862 235/12… 8.04e- 5 6.76e-3 6.37e-3 YWHAH… 103
## 9 GO:0009… cell surfa… 176/3862 435/12… 9.04e- 5 6.76e-3 6.37e-3 EMP2/… 176
## 10 GO:0099… postsynapt… 95/3862 215/12… 1.04e- 4 6.86e-3 6.46e-3 NETO2… 95
## # ℹ 11 more rows
The GeneRatio
and BgRatio
correspond respectively to
the ration between the number of genes of interest in the gene set and the number of (unique) genes-of-interest: \(\frac{|genes\_of\_interest~in~gene\_set|}{|genes\_of\_interest|}\)
the ratio between the number of genes of interest in the gene set and size of the universe, i.e. the number of (unique) genes across all gene sets: \(\frac{|genes~in~gene\_set|}{|universe|}\)
Below, we have 12 genes of interest, "a"
, to "l"
:
## [1] "a" "b" "c" "d" "e" "f" "g" "h" "i" "j" "k" "l"
And two gene sets, "X"
and "Y"
:
## geneset entrez_gene
## 1 X a
## 2 X b
## 3 X c
## 4 X d
## 5 X e
## 6 X f
## 7 X g
## 8 X h
## 9 X i
## 10 X j
## geneset entrez_gene
## 1 Y b
## 2 Y c
## 3 Y d
## 4 Y e
## 5 Y f
## 6 Y g
## 7 Y h
## 8 Y i
## 9 Y j
## 10 Y k
## 11 Y l
## 12 Y m
## 13 Y n
## 14 Y o
## 15 Y p
## 16 Y q
## 17 Y r
## 18 Y s
## 19 Y t
## 20 Y u
## 21 Y v
## 22 Y w
## 23 Y x
## 24 Y y
## 25 Y z
We get:
## GeneRatio BgRatio geneID Count
## X 10/12 10/26 a/b/c/d/e/f/g/h/i/j 10
## Y 11/12 25/26 b/c/d/e/f/g/h/i/j/k/l 11
► Question
Repeat the ORA analysis above using all GO namespaces. See ?enrichGO
for details.
The functions that perform GSEA in clusterProfiler
require the genes
to be ordered in decreasing order. This indicates that the p-values
can’t be used without transformations. One could also use the log2
fold-change or the (absolute value) test statistics.
abs(res_tbl$stat)
ordered_genes <-names(ordered_genes) <- res_tbl$ENTREZID
sort(ordered_genes, decreasing = TRUE)
ordered_genes <-
gseGO(gene = ordered_genes,
go_gsea <-OrgDb = org.Hs.eg.db,
scoreType = "pos") %>%
as_tibble
go_gsea
## # A tibble: 603 × 11
## ID Description setSize enrichmentScore NES pvalue p.adjust qvalue
## <chr> <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 GO:00308… epithelial… 441 0.479 1.40 1.11e-10 5.85e-7 4.46e-7
## 2 GO:00015… blood vess… 493 0.459 1.34 4.34e- 9 1.15e-5 8.73e-6
## 3 GO:00485… blood vess… 434 0.460 1.34 1.99e- 8 3.50e-5 2.66e-5
## 4 GO:00343… cell junct… 499 0.448 1.31 5.93e- 8 7.82e-5 5.96e-5
## 5 GO:00015… angiogenes… 376 0.463 1.35 1.00e- 7 1.06e-4 8.06e-5
## 6 GO:00435… regulation… 394 0.457 1.33 1.21e- 7 1.07e-4 8.13e-5
## 7 GO:00713… cellular r… 472 0.447 1.31 1.53e- 7 1.15e-4 8.79e-5
## 8 GO:00708… response t… 493 0.445 1.30 1.99e- 7 1.31e-4 1.00e-4
## 9 GO:00430… positive r… 383 0.457 1.33 3.33e- 7 1.80e-4 1.37e-4
## 10 GO:00455… negative r… 426 0.450 1.31 3.41e- 7 1.80e-4 1.37e-4
## # ℹ 593 more rows
## # ℹ 3 more variables: rank <int>, leading_edge <chr>, core_enrichment <chr>
► Question
Repeat the GSEA analysis above using all GO namespaces. See
?enrichGO
for details.
► Question
Compare and discuss the use of these different ways to order the genes for the GSEA analysis.
► Question
Compare the GO and GSEA results. Which GO terms are shared or unique to each approach?
► Solution
The enrichKEGG()
and gseKEGG()
functions can be used to perform
the ORA and GSEA analyses againt the KEGG sets.
► Question
Perform ORA and GSEA analysis using the KEGG set.
We will use the hallmark gene set as an example here, but any set described above can be utilised.
library("msigdbr")
msigdbr(species = "Homo sapiens", category = "H") %>%
msig_h <- dplyr::select(gs_name, entrez_gene) %>%
dplyr::rename(ont = gs_name, gene = entrez_gene)
msig_h
## # A tibble: 8,209 × 2
## ont gene
## <chr> <int>
## 1 HALLMARK_ADIPOGENESIS 19
## 2 HALLMARK_ADIPOGENESIS 11194
## 3 HALLMARK_ADIPOGENESIS 10449
## 4 HALLMARK_ADIPOGENESIS 33
## 5 HALLMARK_ADIPOGENESIS 34
## 6 HALLMARK_ADIPOGENESIS 35
## 7 HALLMARK_ADIPOGENESIS 47
## 8 HALLMARK_ADIPOGENESIS 50
## 9 HALLMARK_ADIPOGENESIS 51
## 10 HALLMARK_ADIPOGENESIS 112
## # ℹ 8,199 more rows
The msig_h
table can now be used with the enricher()
and GSEA()
functions from the clusterProfiler
package to perform ORA and GSEA
analyses on the MSigDB hallmark gene set.
enricher(gene = de_genes,
msig_ora <-universe = res_tbl$ENTREZID,
TERM2GENE = msig_h) %>%
as_tibble
msig_ora
## # A tibble: 3 × 9
## ID Description GeneRatio BgRatio pvalue p.adjust qvalue geneID Count
## <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl> <chr> <int>
## 1 HALLMARK_… HALLMARK_H… 82/1212 164/33… 1.46e-4 0.00732 0.00647 10370… 82
## 2 HALLMARK_… HALLMARK_M… 34/1212 58/3347 3.79e-4 0.00948 0.00838 1844/… 34
## 3 HALLMARK_… HALLMARK_A… 53/1212 106/33… 2.17e-3 0.0362 0.0321 3566/… 53
GSEA(ordered_genes, TERM2GENE = msig_h, scoreType = "pos") %>%
msig_gsea <- as_tibble
msig_gsea
## # A tibble: 19 × 11
## ID Description setSize enrichmentScore NES pvalue p.adjust qvalue
## <chr> <chr> <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 HALLMARK_… HALLMARK_E… 128 0.554 1.55 7.29e-8 3.64e-6 1.92e-6
## 2 HALLMARK_… HALLMARK_H… 164 0.509 1.44 2.63e-6 6.58e-5 3.46e-5
## 3 HALLMARK_… HALLMARK_M… 58 0.605 1.62 9.10e-6 1.52e-4 7.98e-5
## 4 HALLMARK_… HALLMARK_P… 178 0.489 1.39 1.58e-5 1.97e-4 1.04e-4
## 5 HALLMARK_… HALLMARK_G… 175 0.477 1.35 6.15e-5 6.15e-4 3.24e-4
## 6 HALLMARK_… HALLMARK_T… 160 0.482 1.37 8.18e-5 6.82e-4 3.59e-4
## 7 HALLMARK_… HALLMARK_U… 122 0.492 1.37 3.90e-4 2.78e-3 1.46e-3
## 8 HALLMARK_… HALLMARK_A… 106 0.491 1.36 6.09e-4 3.05e-3 1.60e-3
## 9 HALLMARK_… HALLMARK_K… 124 0.486 1.36 6.09e-4 3.05e-3 1.60e-3
## 10 HALLMARK_… HALLMARK_X… 134 0.478 1.34 5.91e-4 3.05e-3 1.60e-3
## 11 HALLMARK_… HALLMARK_E… 152 0.470 1.33 7.46e-4 3.39e-3 1.79e-3
## 12 HALLMARK_… HALLMARK_A… 123 0.471 1.32 1.91e-3 7.94e-3 4.18e-3
## 13 HALLMARK_… HALLMARK_T… 50 0.531 1.41 3.83e-3 1.37e-2 7.20e-3
## 14 HALLMARK_… HALLMARK_I… 112 0.472 1.31 3.61e-3 1.37e-2 7.20e-3
## 15 HALLMARK_… HALLMARK_M… 116 0.464 1.29 4.50e-3 1.50e-2 7.89e-3
## 16 HALLMARK_… HALLMARK_I… 147 0.446 1.26 6.95e-3 2.17e-2 1.14e-2
## 17 HALLMARK_… HALLMARK_E… 159 0.436 1.24 1.26e-2 3.71e-2 1.95e-2
## 18 HALLMARK_… HALLMARK_I… 53 0.498 1.33 1.50e-2 3.99e-2 2.10e-2
## 19 HALLMARK_… HALLMARK_C… 68 0.479 1.30 1.52e-2 3.99e-2 2.10e-2
## # ℹ 3 more variables: rank <int>, leading_edge <chr>, core_enrichment <chr>
The clusterProfiler
documentation provides a chapter on the
visualization of functional enrichment
results.
Another useful visualisation, that links the enrichment results back to the whole set of results is to highlight the genes in a particular set of interest on the volcano plot.
res_tbl$ENTREZID %in% GO_0005925$ENTREZID
sel <-plot(res_tbl$log2FoldChange, -log10(res_tbl$padj))
points(res_tbl$log2FoldChange[sel],
-log10(res_tbl$padj)[sel],
col = "red")
grid()
► Question
Extract the Entrez ids from the top hit in the msig_gsea
results
above (available in core_enrichment
) and visualise the results on a
volcano plot.
► Solution
Enrichment analyses are powerful techniques, as they allow to integrate thousands of univariate results into biologically driven sets.
Their interpretation isn’t however always straightforward when different methods (ORA and GSEA) and implementations (different GSEA permutation strategies) or user-set parameters (DE cutoffs in ORA, ordering criteria in GSEA) provide different results.
Page built: 2024-10-04 using R version 4.4.1 (2024-06-14)