data(geneList, package="DOSE")
head(geneList) 4312 8318 10874 55143 55388 991
4.572613 4.514594 4.418218 4.144075 3.876258 3.677857 11 Universal enrichment analysis
The clusterProfiler package (Yu et al. 2012) supports both hypergeometric test and gene set enrichment analyses of many ontologies/pathways, but it’s still not enough as users may want to analyze their data with unsupported organisms, slim versions of GO, novel functional annotations (e.g. GO via BlastGO or KEGG via KAAS), unsupported ontologies/pathways, or customized annotations.
The clusterProfiler package provides enricher() function for hypergeometric test and GSEA() function for gene set enrichment analysis that are designed to accept user-defined annotations. They accept two additional parameters TERM2GENE and TERM2NAME. As indicated in the parameter names, TERM2GENE is a data.frame with the first column of term ID and the second column of corresponding mapped genes, and TERM2NAME is a data.frame with the first column of term ID and the second column of corresponding term names. TERM2NAME is optional.
11.1 Input data
For over representation analysis, all we need is a gene vector, that is a vector of gene IDs. These gene IDs can be obtained by differential expression analysis (e.g. with the DESeq2 package).
For gene set enrichment analysis, we need a ranked list of genes. DOSE provides an example dataset geneList which was derived from R package breastCancerMAINZ that contained 200 samples, including 29 samples in grade I, 136 samples in grade II and 35 samples in grade III. We computed the ratios of geometric means of grade III samples versus geometric means of grade I samples. Logarithm of these ratios (base 2) were stored in geneList dataset. If you want to prepare your own geneList, please refer to the FAQ.
We can load the sample data into R via:
Suppose we define fold change greater than 2 as DEGs:
gene <- names(geneList)[abs(geneList) > 2]
head(gene)[1] "4312" "8318" "10874" "55143" "55388" "991" 11.2 Cell Marker
# not executed: it downloads an ~8 MB file, and the URL changes between
# CellMarker releases (see the note below).
## The original download endpoint no longer resolves:
## http://bio-bigdata.hrbmu.edu.cn/CellMarker/download/Human_cell_markers.txt
## The file is now served from the CellMarker download page
## http://117.50.127.228/CellMarker/CellMarker_download.html
## This direct link worked when this chapter was last checked:
url <- "http://www.bio-bigdata.center/CellMarker_download_files/file/Cell_marker_Human.xlsx"
f <- tempfile(fileext = ".xlsx")
download.file(url, f, mode = "wb")
cell_marker_data <- readxl::read_excel(f, 1)
## Column names differ between CellMarker releases (`Cell name` in older
## exports, `cell_name` in the current one), so check before selecting:
## names(cell_marker_data)
## instead of `cell_name`, users can use other features (e.g. `cancer_type`)
cells <- cell_marker_data %>%
dplyr::select(cell_name, GeneID) %>%
dplyr::mutate(GeneID = strsplit(GeneID, ", ")) %>%
tidyr::unnest(cols = c(GeneID))The CellMarker download URL is not stable. The host has moved more than once (bio-bigdata.hrbmu.edu.cn now refuses connections; the mirror that used to be cited here, yikedaxue.slwshop.cn, no longer resolves), and the column headers have been renamed between releases — older exports used Cell name and cellType, the current one uses cell_name and cell_type. If the code above fails, download the file from the CellMarker download page by hand and run names(cell_marker_data) to see which names your copy uses, then adjust the dplyr::select() call accordingly.
11.2.1 Cell Marker over-representation analysis
x <- enricher(gene, TERM2GENE = cells)
head(x)11.2.2 Cell Marker gene set enrichment analysis
y <- GSEA(geneList, TERM2GENE = cells)
head(y)11.3 MSigDb analysis
Molecular Signatures Database is a collection of annotated gene sets. It contains 8 major collections:
- H: hallmark gene sets
- C1: positional gene sets
- C2: curated gene sets
- C3: motif gene sets
- C4: computational gene sets
- C5: GO gene sets
- C6: oncogenic signatures
- C7: immunologic signatures
Users can download GMT files from Broad Institute and use the read.gmt() function to parse the file to be used in enricher() and GSEA().
There is an R package, msigdbr, that already packed the MSigDB gene sets in tidy data format that can be used directly with clusterProfiler (Yu et al. 2012).
It supports several species:
## Must not be cached: attaching the package is a side effect, and the chunks
## below call msigdbr() unqualified (they are cache: false themselves). With a
## warm cache this chunk would be skipped, the package would never be attached,
## and those chunks would fail with "could not find function msigdbr".
library(msigdbr)
msigdbr_species()# A tibble: 20 × 2
species_name species_common_name
<chr> <chr>
1 Anolis carolinensis Carolina anole, green anole
2 Bos taurus bovine, cattle, cow, dairy cow, domestic cat…
3 Caenorhabditis elegans <NA>
4 Canis lupus familiaris dog, dogs
5 Danio rerio leopard danio, zebra danio, zebra fish, zebr…
6 Drosophila melanogaster fruit fly
7 Equus caballus domestic horse, equine, horse
8 Felis catus cat, cats, domestic cat
9 Gallus gallus bantam, chicken, chickens, Gallus domesticus
10 Homo sapiens human
11 Macaca mulatta rhesus macaque, rhesus macaques, Rhesus monk…
12 Monodelphis domestica gray short-tailed opossum
13 Mus musculus house mouse, mouse
14 Ornithorhynchus anatinus duck-billed platypus, duckbill platypus, pla…
15 Pan troglodytes chimpanzee
16 Rattus norvegicus brown rat, Norway rat, rat, rats
17 Saccharomyces cerevisiae baker's yeast, brewer's yeast, S. cerevisiae
18 Schizosaccharomyces pombe 972h- <NA>
19 Sus scrofa pig, pigs, swine, wild boar
20 Xenopus tropicalis tropical clawed frog, western clawed frog We can retrieve human gene sets with msigdbr(). Since the complete human MSigDB table is large, it is usually better to request only the collection you need in analysis workflows.
msigdbr() downloads the gene set archive from a remote host, and that download can fail transiently, so these calls are wrapped in retry() (defined in _common.R); feel free to drop the wrapper in your own analysis.
retry(msigdbr(species = "Homo sapiens", collection = "C6")) |>
head(2) |>
as.data.frame() gene_symbol ncbi_gene ensembl_gene db_gene_symbol db_ncbi_gene
1 ACKR3 57007 ENSG00000144476 ACKR3 57007
2 ADGRL1 22859 ENSG00000072071 ADGRL1 22859
db_ensembl_gene source_gene gs_id gs_name gs_collection gs_subcollection
1 ENSG00000144476 CXCR7 M2666 AKT_UP.V1_DN C6
2 ENSG00000072071 LPHN1 M2666 AKT_UP.V1_DN C6
gs_collection_name
1 Oncogenic Signature
2 Oncogenic Signature
gs_description
1 Genes down-regulated in mouse prostate by transgenic expression of human AKT1 gene [Gene ID=207] vs controls.
2 Genes down-regulated in mouse prostate by transgenic expression of human AKT1 gene [Gene ID=207] vs controls.
gs_source_species gs_pmid gs_geoid
1 MM 15156201 GSE1413
2 MM 15156201 GSE1413
gs_exact_source gs_url
1 AKT_PLACEBO vs WT_PLACEBO; bottom 150 genes (diff. of means)
2 AKT_PLACEBO vs WT_PLACEBO; bottom 150 genes (diff. of means)
db_version db_target_species
1 2026.1.Hs HS
2 2026.1.Hs HSOr specific collection. Here we use C6, oncogenic gene sets as an example:
m_t2g <- retry(msigdbr(species = "Homo sapiens", collection = "C6")) |>
dplyr::select(gs_name, ncbi_gene)
head(m_t2g)# A tibble: 6 × 2
gs_name ncbi_gene
<chr> <chr>
1 AKT_UP.V1_DN 57007
2 AKT_UP.V1_DN 22859
3 AKT_UP.V1_DN 137872
4 AKT_UP.V1_DN 249
5 AKT_UP.V1_DN 271
6 AKT_UP.V1_DN 51129 11.3.1 MSigDb over-representation analysis
em <- enricher(gene, TERM2GENE=m_t2g)
head(em) ID Description GeneRatio
RPS14_DN.V1_DN RPS14_DN.V1_DN RPS14_DN.V1_DN 22/183
GCNP_SHH_UP_LATE.V1_UP GCNP_SHH_UP_LATE.V1_UP GCNP_SHH_UP_LATE.V1_UP 16/183
PRC2_EZH2_UP.V1_DN PRC2_EZH2_UP.V1_DN PRC2_EZH2_UP.V1_DN 15/183
VEGF_A_UP.V1_DN VEGF_A_UP.V1_DN VEGF_A_UP.V1_DN 15/183
RB_P107_DN.V1_UP RB_P107_DN.V1_UP RB_P107_DN.V1_UP 10/183
E2F1_UP.V1_UP E2F1_UP.V1_UP E2F1_UP.V1_UP 12/183
BgRatio RichFactor FoldEnrichment oddsRatio zScore
RPS14_DN.V1_DN 186/10927 0.11827957 7.062518 8.815331 10.883293
GCNP_SHH_UP_LATE.V1_UP 181/10927 0.08839779 5.278266 6.142769 7.574549
PRC2_EZH2_UP.V1_DN 192/10927 0.07812500 4.664874 5.330408 6.686237
VEGF_A_UP.V1_DN 193/10927 0.07772021 4.640703 5.299960 6.659725
RB_P107_DN.V1_UP 130/10927 0.07692308 4.593106 5.117534 5.378527
E2F1_UP.V1_UP 188/10927 0.06382979 3.811301 4.213716 5.074316
pvalue p.adjust qvalue
RPS14_DN.V1_DN 4.615656e-13 7.800459e-11 7.800459e-11
GCNP_SHH_UP_LATE.V1_UP 5.729300e-08 4.841259e-06 4.841259e-06
PRC2_EZH2_UP.V1_DN 7.531662e-07 3.400322e-05 3.400322e-05
VEGF_A_UP.V1_DN 8.048099e-07 3.400322e-05 3.400322e-05
RB_P107_DN.V1_UP 6.370479e-05 2.005313e-03 2.005313e-03
E2F1_UP.V1_UP 7.445084e-05 2.005313e-03 2.005313e-03
geneID
RPS14_DN.V1_DN 9787/5105/54821/79901/55388/79733/10112/7021/51659/4174/1062/2532/10874/9493/9055/9133/83461/23397/9319/55872/991/4605
GCNP_SHH_UP_LATE.V1_UP 9787/55388/79733/6241/3169/1111/4174/51203/3833/6790/1580/9212/7153/9055/983/9319
PRC2_EZH2_UP.V1_DN 55355/9787/55388/3887/11182/6362/10460/2146/24137/9212/81620/23397/8318/7272/4605
VEGF_A_UP.V1_DN 9787/891/3832/6241/10234/10403/1062/2167/9493/9133/8318/7272/23362/332/4085
RB_P107_DN.V1_UP 79733/6241/4174/24137/9055/23397/8208/8318/1307/4085
E2F1_UP.V1_UP 9787/55388/79733/1111/3833/2146/9212/7153/9055/23397/2842/8208
Count
RPS14_DN.V1_DN 22
GCNP_SHH_UP_LATE.V1_UP 16
PRC2_EZH2_UP.V1_DN 15
VEGF_A_UP.V1_DN 15
RB_P107_DN.V1_UP 10
E2F1_UP.V1_UP 1211.3.2 MSigDb gene set enrichment analysis
In over-representation analysis, we use oncogenic gene sets (i.e. C6) to test whether the DE genes are involved in the process that leads to cancer. In this example, we will use the C3 category to test whether genes are up/down-regulated by sharing specific motifs using the GSEA approach.
C3_t2g <- retry(msigdbr(species = "Homo sapiens", collection = "C3")) |>
dplyr::select(gs_name, ncbi_gene)
head(C3_t2g)# A tibble: 6 × 2
gs_name ncbi_gene
<chr> <chr>
1 AAACCAC_MIR140 10257
2 AAACCAC_MIR140 23172
3 AAACCAC_MIR140 81
4 AAACCAC_MIR140 90
5 AAACCAC_MIR140 8754
6 AAACCAC_MIR140 11096 em2 <- GSEA(geneList, TERM2GENE = C3_t2g)
head(em2) ID Description setSize
HSD17B8_TARGET_GENES HSD17B8_TARGET_GENES HSD17B8_TARGET_GENES 400
SGCGSSAAA_E2F1DP2_01 SGCGSSAAA_E2F1DP2_01 SGCGSSAAA_E2F1DP2_01 129
E2F1_Q3 E2F1_Q3 E2F1_Q3 187
E2F1_Q6 E2F1_Q6 E2F1_Q6 185
E2F_Q4 E2F_Q4 E2F_Q4 193
E2F4DP1_01 E2F4DP1_01 E2F4DP1_01 189
enrichmentScore NES pvalue p.adjust
HSD17B8_TARGET_GENES 0.6375083 3.065600 1.000000e-10 1.650500e-07
SGCGSSAAA_E2F1DP2_01 0.5408184 2.294678 2.775259e-10 2.290283e-07
E2F1_Q3 0.4983913 2.224661 1.000000e-10 1.650500e-07
E2F1_Q6 0.4906658 2.174662 2.015871e-10 2.218130e-07
E2F_Q4 0.4685394 2.094634 2.997241e-09 1.413413e-06
E2F4DP1_01 0.4646480 2.074758 1.698289e-09 1.121211e-06
qvalue rank leading_edge
HSD17B8_TARGET_GENES 1.228259e-07 1042 tags=36%, list=8%, signal=34%
SGCGSSAAA_E2F1DP2_01 1.704368e-07 1797 tags=38%, list=14%, signal=33%
E2F1_Q3 1.228259e-07 2045 tags=36%, list=16%, signal=30%
E2F1_Q6 1.650674e-07 1819 tags=34%, list=15%, signal=30%
E2F_Q4 1.051825e-06 1663 tags=29%, list=13%, signal=26%
E2F4DP1_01 8.343753e-07 1797 tags=33%, list=14%, signal=29%
core_enrichment
HSD17B8_TARGET_GENES 55143/55388/991/2305/9493/1062/4605/9833/9133/10403/7153/23397/79733/259266/6241/55165/9787/11065/55355/9582/220134/55872/83461/10460/4751/10635/55839/890/9415/983/54821/4085/9837/5080/332/7272/64151/2842/9212/51659/9319/9055/3833/146909/10112/51514/4174/9232/1033/9928/3161/11004/993/4603/79801/990/5347/55215/701/55723/51512/55635/9156/11130/10024/57405/10615/1894/2491/8438/9700/5888/7083/898/3149/56992/9768/4288/10733/1163/4175/5307/2237/29899/55010/59336/56942/5984/29028/3838/1058/54478/3015/699/6491/81611/1063/27346/64785/9401/26271/51026/641/1869/10535/1029/28998/1763/8970/54892/55159/864/93594/1789/4176/7111/3148/116832/9585/9735/55732/81624/23310/1871/1031/79915/995/10051/30012/1104/80178/78995/7468/284403/5230/8549/7884/5558/4172/5424/79866/3182/56902/83990
SGCGSSAAA_E2F1DP2_01 79733/6241/4171/993/990/3159/9768/4175/4173/29028/3015/4609/85236/5111/64785/26271/51053/1869/5427/4176/5902/4436/1871/79915/9088/6632/284403/7290/4172/5424/124222/1633/7374/2189/4678/2289/50/63967/56886/23234/6839/1786/204/51087/23636/1174/83463/10492/79173
E2F1_Q3 79733/6241/983/5080/2146/8715/4171/993/990/56992/9768/4998/10733/7037/4175/54962/54954/29028/3015/2643/85236/5111/64785/26271/51053/1869/3925/114/5427/4176/7112/4436/1871/79915/9088/6632/284403/7290/2597/7398/4172/5424/124222/1633/10849/9824/4678/6646/29902/63967/56886/6059/23234/1786/204/51087/23636/2048/83463/10393/79173/688/3007/79677/9462/56946/9221
E2F1_Q6 79733/6241/983/5080/81620/2146/4171/993/990/3159/9768/4998/4175/4173/29028/3015/4609/85236/5111/64785/26271/51053/1869/3925/5427/4176/7112/5902/4436/1871/79915/9088/6632/284403/7290/2597/4172/5424/124222/1633/10849/7374/9824/2189/4678/2289/6646/29902/50/63967/56886/23234/6839/1786/204/51087/23636/11335/1174/83463/10492/79173/688
E2F_Q4 8318/79733/6241/983/5080/81620/2146/4171/993/990/5888/9768/4175/4173/29028/3015/85236/5111/64785/26271/51053/1869/3925/5427/23649/4176/7112/5902/4436/1871/79915/9088/6632/284403/7290/2597/4172/5424/124222/1633/10849/7473/7374/9824/7353/4678/6646/50/63967/56886/23234/1786/204/6427/51087/23636
E2F4DP1_01 79733/6241/983/5080/2146/4171/993/990/3159/9768/4175/4173/29028/3015/4609/85236/5111/64785/26271/51053/1869/3925/5427/57103/4176/7112/5902/4436/1871/79915/9088/64211/6632/284403/7290/2597/4172/5424/124222/1633/10849/7374/9824/2189/4678/2289/6646/50/63967/5457/56886/23234/6839/1786/204/51087/23636/11335/1174/83463/10492/79173
log2err
HSD17B8_TARGET_GENES NaN
SGCGSSAAA_E2F1DP2_01 0.8140358
E2F1_Q3 NaN
E2F1_Q6 0.8266573
E2F_Q4 0.7749390
E2F4DP1_01 0.788186811.4 Proteomics data and UniProt identifiers
Proteomics experiments normally report UniProt accessions rather than Entrez gene IDs, and the organism is frequently not one of the model organisms that have an OrgDb package. This is the situation where the usual enrichGO()-style entry points cannot be used directly, so it is worth spelling out the two ways out.
11.4.1 Model organisms: convert the identifiers
If an OrgDb exists for your organism and carries the UNIPROT keytype, the simplest solution is to translate the accessions and then use the standard pipeline:
library(org.Hs.eg.db)
up <- c("P04637", "P00533", "Q9NQ11", "P38398", "O14746")
bitr(up, fromType = "UNIPROT", toType = "ENTREZID", OrgDb = org.Hs.eg.db) UNIPROT ENTREZID
1 P04637 7157
2 P00533 1956
3 Q9NQ11 23400
4 P38398 672
5 O14746 7015Use keytypes(org.Hs.eg.db) to check what is available for your organism — not every OrgDb provides UNIPROT.
11.4.2 Any organism: build TERM2GENE from the EBI GOA files
For a non-model organism there is no OrgDb to convert with, but you do not need one. The Gene Ontology Annotation (GOA) group at EMBL-EBI publishes GO annotation for complete proteomes in the standard GAF format, keyed by UniProt accession:
- Browse the available proteomes: https://www.ebi.ac.uk/GOA/proteomes
- Download directory: http://ftp.ebi.ac.uk/pub/databases/GO/goa/proteomes/
A GAF file is plain tab-delimited text, so it can be read with read.delim(). In the GAF 2.x layout the columns we need are: 2 = UniProt accession, 4 = qualifier, 5 = GO ID, and 9 = ontology aspect (P/C/F).
## example: Solanum lycopersicum (tomato), proteome UP000004994
url <- "http://ftp.ebi.ac.uk/pub/databases/GO/goa/proteomes/284179.S_lycopersicum_genome_cv_Heinz_1706.goa"
goa_file <- "tomato.goa"
download.file(url, goa_file)library(GO.db)
library(AnnotationDbi)
gaf <- read.delim(goa_file, header = FALSE, comment.char = "!", sep = "\t",
quote = "", fill = TRUE, stringsAsFactors = FALSE)
colnames(gaf) <- c("DB", "UniProt", "Symbol", "Qualifier", "GO", "Reference",
"Evidence", "With", "Aspect", "Name", "Synonym", "Type",
"Taxon", "Date", "AssignedBy", "Extension", "GeneProductForm")
## drop negative annotations, which are written as "NOT|..."
gaf <- gaf[!grepl("^NOT", gaf$Qualifier), ]
## TERM2GENE: (GO id, UniProt accession)
term2gene <- unique(gaf[, c("GO", "UniProt")])
## TERM2NAME: GO term names come from GO.db
goterm <- AnnotationDbi::select(GO.db, keys = unique(term2gene$GO),
columns = c("TERM", "ONTOLOGY"), keytype = "GOID")
term2name <- goterm[, c("GOID", "TERM")]These two data frames plug straight into enricher() and GSEA(), exactly like the MSigDb example above — the only difference is that the gene identifiers are UniProt accessions:
res <- enricher(gene = my_proteins,
universe = all_proteins,
TERM2GENE = term2gene,
TERM2NAME = term2name)Two details are easy to get wrong:
universeshould be the proteins detected in your experiment, not the entire proteome. The GOA file annotates the full proteome, but the relevant background for a mass-spectrometry experiment is the set of proteins you could have quantified. Passing the whole proteome instead inflates the significance of everything.- Aspect filtering. If you only want one ontology, filter
term2geneby theONTOLOGYcolumn ofgotermbefore callingenricher():
bp_ids <- goterm$GOID[goterm$ONTOLOGY == "BP"]
term2gene_bp <- term2gene[term2gene$GO %in% bp_ids, ]If you need curated annotations only, filter on the evidence code as well — the GOA files include electronically inferred (IEA) annotations, and gaf[gaf$Evidence != "IEA", ] keeps just the manually reviewed ones.
11.4.3 KEGG analysis with UniProt identifiers
KEGG does accept UniProt accessions directly, so no conversion is needed for organisms that KEGG covers:
kk <- retry(enrichKEGG(gene = my_proteins, organism = "sly", keyType = "uniprot"))Use bitr_kegg() to check the coverage first — for most proteomes a substantial fraction of accessions has no KEGG entry, and those proteins are silently dropped:
uni2kegg <- bitr_kegg(unique(gaf$UniProt), fromType = "uniprot",
toType = "kegg", organism = "sly")References
Yu, Guangchuang, Le-Gen Wang, Yanyan Han, and Qing-Yu He. 2012. “clusterProfiler: An r Package for Comparing Biological Themes Among Gene Clusters.” OMICS: A Journal of Integrative Biology 16 (5): 284–87. https://doi.org/10.1089/omi.2011.0118.