8 ORA analysis
As the paper “Urgent need for consistent standards in functional enrichment analysis” mentioned, popular functional enrichment tools can be classified into two main categories: (1) over-representation analysis (ORA) and (2) functional class scoring (FCS).
Kindly Note: No matter which method you choose, the first thing is always the preparation of gene sets.
8.1 ORA method introduction
According to paper “Over-representation of correlation analysis (ORCA): a method for identifying associations between variable sets”:
Over-representation analysis (ORA) is a simple method for objectively deciding whether a set of variables of known or suspected biological relevance, such as a gene set or pathway, is more prevalent in a set of variables of interest than we expect by chance.
ORA usually has these procedures:
- Identify features, such as mRNA, proteins or microarray probes, considered “differentially expressed” in two or more conditions
- Determine the number of differentially expressed features in each pathway
- For each pathway, using a hypergeometric distribution to calculate a probability value (P-value) of obtaining the number of differentially expressed features against the background list of all features
8.2 Basic usage
The simplest arguments are:
-
id: gene id of Entrez (numeric id is also supported), Ensembl, Symbol (alias is also accepted) or Uniprot -
geneset: gene set is a two-column data frame with term id and gene id. (It’s recommended to use geneset) -
p_cutoff: numeric of cutoff for both pvalue and adjusted pvalue, default is 0.05. -
q_cutoff: numeric of cutoff for qvalue, default is 0.15.
8.2.1 1st step: prepare input IDs
To show the various input ID types, we choose Entrez ID, Symbol from human and Uniprot from rat protein example data.
data(geneList, package = "genekitr")
entrez_id <- names(geneList)[abs(geneList) > 2]
head(entrez_id, 5)## [1] "948" "1638" "158471" "10610" "6947"
symbol_id <- c("TP53","BRCA1","TET2") # all official symbols
mix_symbol_id <- c("BCC7","BRCA1","TET2") # mixture of official symbol and alias
unip_id <- rat_prodata$Acc[1:100]
head(unip_id, 5)## [1] "P27364" "P17475" "P14046" "P00507" "P15650"8.2.2 2nd step: prepare gene set
For more details, please refer to chapter 7
Here we prepare GO MF gene set:
8.2.3 3rd step: ORA analysis
For non-symbol genes (e.g. Entrez IDs):
## Hs_MF_ID Description
## GO:0000217 GO:0000217 DNA secondary structure binding
## GO:0005041 GO:0005041 low-density lipoprotein particle receptor activity
## GO:0030228 GO:0030228 lipoprotein particle receptor activity
## GeneRatio BgRatio RichFactor FoldEnrich oddsRatio zScore
## GO:0000217 0.032 37/18845 0.1081081 16.29838 18.71976 7.611383
## GO:0005041 0.024 15/18845 0.2000000 30.15200 38.33607 9.229485
## GO:0030228 0.024 18/18845 0.1666667 25.12667 30.66393 8.368192
## pvalue p.adjust qvalue geneID
## GO:0000217 0.0001028428 0.02244063 0.001208787 641/4208/63967/7516
## GO:0005041 0.0001222922 0.02244063 0.001208787 53353/948/3949
## GO:0030228 0.0002161518 0.02644257 0.001424355 53353/948/3949
## geneID_symbol Count
## GO:0000217 BLM/MEF2C/CLSPN/XRCC2 4
## GO:0005041 LRP1B/CD36/LDLR 3
## GO:0030228 LRP1B/CD36/LDLR 3For symbol genes without alias:
## Hs_MF_ID Description
## GO:0002039 GO:0002039 p53 binding
## GO:0140097 GO:0140097 catalytic activity, acting on DNA
## GO:0031625 GO:0031625 ubiquitin protein ligase binding
## GO:0044389 GO:0044389 ubiquitin-like protein ligase binding
## GO:0001094 GO:0001094 TFIID-class transcription factor complex binding
## GO:0140666 GO:0140666 annealing activity
## GeneRatio BgRatio RichFactor FoldEnrich oddsRatio zScore
## GO:0002039 0.6666667 71/18845 0.028169014 176.94836 544.1449 18.742161
## GO:0140097 0.6666667 270/18845 0.007407407 46.53086 138.6119 9.508391
## GO:0031625 0.6666667 305/18845 0.006557377 41.19126 122.3696 8.929154
## GO:0044389 0.6666667 325/18845 0.006153846 38.65641 114.6687 8.640598
## GO:0001094 0.3333333 10/18845 0.100000000 628.16667 1046.2778 25.031325
## GO:0140666 0.3333333 11/18845 0.090909091 571.06061 941.6000 23.863260
## pvalue p.adjust qvalue geneID Count
## GO:0002039 4.188388e-05 0.002387381 0.0003967946 BRCA1/TP53 2
## GO:0140097 6.077577e-04 0.012531517 0.0020828007 TP53/TET2 2
## GO:0031625 7.748981e-04 0.012531517 0.0020828007 BRCA1/TP53 2
## GO:0044389 8.794047e-04 0.012531517 0.0020828007 BRCA1/TP53 2
## GO:0001094 1.591174e-03 0.015546393 0.0025838880 TP53 1
## GO:0140666 1.750198e-03 0.015546393 0.0025838880 TP53 1For symbol genes with alias:
## Hs_MF_ID Description
## GO:0002039 GO:0002039 p53 binding
## GO:0140097 GO:0140097 catalytic activity, acting on DNA
## GO:0031625 GO:0031625 ubiquitin protein ligase binding
## GO:0044389 GO:0044389 ubiquitin-like protein ligase binding
## GO:0001094 GO:0001094 TFIID-class transcription factor complex binding
## GO:0140666 GO:0140666 annealing activity
## GeneRatio BgRatio RichFactor FoldEnrich oddsRatio zScore
## GO:0002039 0.6666667 71/18845 0.028169014 176.94836 544.1449 18.742161
## GO:0140097 0.6666667 270/18845 0.007407407 46.53086 138.6119 9.508391
## GO:0031625 0.6666667 305/18845 0.006557377 41.19126 122.3696 8.929154
## GO:0044389 0.6666667 325/18845 0.006153846 38.65641 114.6687 8.640598
## GO:0001094 0.3333333 10/18845 0.100000000 628.16667 1046.2778 25.031325
## GO:0140666 0.3333333 11/18845 0.090909091 571.06061 941.6000 23.863260
## pvalue p.adjust qvalue geneID geneID_symbol Count
## GO:0002039 4.188388e-05 0.002387381 0.0003967946 BRCA1/BCC7 BRCA1/TP53 2
## GO:0140097 6.077577e-04 0.012531517 0.0020828007 BCC7/TET2 TP53/TET2 2
## GO:0031625 7.748981e-04 0.012531517 0.0020828007 BRCA1/BCC7 BRCA1/TP53 2
## GO:0044389 8.794047e-04 0.012531517 0.0020828007 BRCA1/BCC7 BRCA1/TP53 2
## GO:0001094 1.591174e-03 0.015546393 0.0025838880 BCC7 TP53 1
## GO:0140666 1.750198e-03 0.015546393 0.0025838880 BCC7 TP53 1For Uniprot proteins:
## Rn_MF_ID
## GO:0004497 GO:0004497
## GO:0020037 GO:0020037
## GO:0046906 GO:0046906
## GO:0016712 GO:0016712
## GO:0016705 GO:0016705
## GO:0008395 GO:0008395
## Description
## GO:0004497 monooxygenase activity
## GO:0020037 heme binding
## GO:0046906 tetrapyrrole binding
## GO:0016712 oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen, reduced flavin or flavoprotein as one donor, and incorporation of one atom of oxygen
## GO:0016705 oxidoreductase activity, acting on paired donors, with incorporation or reduction of molecular oxygen
## GO:0008395 steroid hydroxylase activity
## GeneRatio BgRatio RichFactor FoldEnrich oddsRatio zScore
## GO:0004497 0.1978022 169/20570 0.10650888 24.07569 33.19459 20.07914
## GO:0020037 0.1868132 183/20570 0.09289617 20.99862 28.11145 18.11431
## GO:0046906 0.1868132 192/20570 0.08854167 20.01431 26.65390 17.64507
## GO:0016712 0.1318681 75/20570 0.16000000 36.16703 49.22483 20.33833
## GO:0016705 0.1868132 257/20570 0.06614786 14.95232 19.37292 15.00372
## GO:0008395 0.0989011 37/20570 0.24324324 54.98366 80.16507 21.90836
## pvalue p.adjust qvalue
## GO:0004497 3.251908e-20 9.983359e-18 6.195741e-18
## GO:0020037 4.175578e-18 6.409512e-16 3.977787e-16
## GO:0046906 9.502882e-18 9.724616e-16 6.035164e-16
## GO:0016712 5.516482e-16 4.233900e-14 2.627588e-14
## GO:0016705 1.305185e-15 8.013833e-14 4.973440e-14
## GO:0008395 4.860771e-14 2.487094e-12 1.543508e-12
## geneID
## GO:0004497 P10633/P13107/P17178/P04800/P05179/P20816/P08683/P24483/P00176/P10634/P11715/P12939/P15149/P19225/P11510/P20814/P23457/P18125
## GO:0020037 P10633/P04762/P13107/P17178/P05179/P20816/P08683/P01946/P00176/P10634/P11715/P12939/P15149/P19225/P11510/P20814/P18125
## GO:0046906 P10633/P04762/P13107/P17178/P05179/P20816/P08683/P01946/P00176/P10634/P11715/P12939/P15149/P19225/P11510/P20814/P18125
## GO:0016712 P10633/P13107/P17178/P05179/P20816/P08683/P00176/P11715/P15149/P11510/P20814/P18125
## GO:0016705 P10633/P13107/P17178/P05179/P20816/P08683/P24483/P00176/P10634/P11715/P12939/P15149/P19225/P11510/P20814/P23457/P18125
## GO:0008395 P13107/P17178/P04800/P08683/P24483/P00176/P11715/P15149/P18125
## geneID_symbol
## GO:0004497 Cyp2d1/Cyp2b3/Cyp27a1/Cyp3a1/Cyp2c7/Cyp4a2/Cyp2c11/Fdx1/Cyp2b1/Cyp2d26/Cyp17a1/Cyp2d10/Cyp2a2/Cyp2c70/Cyp2c12/Cyp2c13/Akr1c9/Cyp7a1
## GO:0020037 Cyp2d1/Cat/Cyp2b3/Cyp27a1/Cyp2c7/Cyp4a2/Cyp2c11/Hba1/Cyp2b1/Cyp2d26/Cyp17a1/Cyp2d10/Cyp2a2/Cyp2c70/Cyp2c12/Cyp2c13/Cyp7a1
## GO:0046906 Cyp2d1/Cat/Cyp2b3/Cyp27a1/Cyp2c7/Cyp4a2/Cyp2c11/Hba1/Cyp2b1/Cyp2d26/Cyp17a1/Cyp2d10/Cyp2a2/Cyp2c70/Cyp2c12/Cyp2c13/Cyp7a1
## GO:0016712 Cyp2d1/Cyp2b3/Cyp27a1/Cyp2c7/Cyp4a2/Cyp2c11/Cyp2b1/Cyp17a1/Cyp2a2/Cyp2c12/Cyp2c13/Cyp7a1
## GO:0016705 Cyp2d1/Cyp2b3/Cyp27a1/Cyp2c7/Cyp4a2/Cyp2c11/Fdx1/Cyp2b1/Cyp2d26/Cyp17a1/Cyp2d10/Cyp2a2/Cyp2c70/Cyp2c12/Cyp2c13/Akr1c9/Cyp7a1
## GO:0008395 Cyp2b3/Cyp27a1/Cyp3a1/Cyp2c11/Fdx1/Cyp2b1/Cyp17a1/Cyp2a2/Cyp7a1
## Count
## GO:0004497 18
## GO:0020037 17
## GO:0046906 17
## GO:0016712 12
## GO:0016705 17
## GO:0008395 98.2.4 Some features
We can draw some conclusions about the following features based on the examples above:
- No need to specify input ID type and organism. genekitr will automatically detect the organism name and convert IDs to the suitable type.
- Returns a data.frame object which is easy to filter or export
- Non-symbol IDs will be automatically converted to the latest gene symbol for easier interpretation
- If input symbols are mixed with aliases, a new column “geneID_symbol” will be added; if all symbols are official, only the “geneID” column will be returned
-
genORA()includes five gene-level measurements: “Count”, “BgRatio”, “GeneRatio”, “FoldEnrich” and “RichFactor” (the “oddsRatio” and “zScore” columns from clusterProfiler >= 4.14 are kept as well).-
BgRatio: Number of all genes in a specific term / Number of universal genes (if not specified, this is all genes in the gene set) -
GeneRatio: Number of genes enriched in a specific term / Number of input genes -
FoldEnrich: GeneRatio / BgRatio -
RichFactor: Number of genes enriched in a specific term / Number of all genes in that specific term
-
8.3 Advanced usage
8.3.1 Additional arguments
-
padj_method: choose from “holm”, “hochberg”, “hommel”, “bonferroni”, “BH” (default), “BY”, “fdr”, “none” -
min_gset_size: Minimal size of each gene set for analysis. Default is 10. Pathways with fewer than ten genes will be omitted (e.g. GO:0062196 is a child term with 8 genes at the bottom of the directed acyclic graph). -
max_gset_size: Maximal size. Default is 500. Pathways with more than 500 genes will be omitted (e.g. GO:0007049 has 1796 genes as parent nodes at the top of the directed acyclic graph). -
universe: If users have their own background genes, set them asuniverse. If missing, all genes in the geneset will be used as the background universe.
8.3.2 Group enrichment
If there is a multi-group comparison (e.g. multiple experimental conditions or multiple time points), users can specify gene groups using the group_list argument:
For example, we have 100 genes:
- 50/100 genes are up-regulated and 50/100 are down-regulated
- 40/100 genes are detected at time point 1 and 60/100 are detected at time point 2
id_100 <- c(head(names(geneList), 50), tail(names(geneList), 50))
two_groups <- list(
exp_group = c(rep("up", 50), rep("down", 50)),
time_group = c(rep("time1", 40), rep("time2", 60))
)
gora <- genORA(id_100, geneset = hg_gs, group_list = two_groups)
head(gora)## Hs_MF_ID Cluster Description GeneRatio
## 1 GO:0140612 down.time2 DNA damage sensor activity 0.06122449
## 2 GO:0140097 down.time2 catalytic activity, acting on DNA 0.12244898
## 3 GO:0000217 down.time2 DNA secondary structure binding 0.06122449
## 4 GO:0003697 down.time2 single-stranded DNA binding 0.08163265
## 5 GO:0032050 down.time2 clathrin heavy chain binding 0.04081633
## 6 GO:0008094 down.time2 ATP-dependent activity, acting on DNA 0.08163265
## BgRatio RichFactor FoldEnrich oddsRatio zScore pvalue p.adjust
## 1 25/18845 0.12000000 46.151020 55.654150 11.533992 3.649624e-05 0.006755819
## 2 270/18845 0.02222222 8.546485 9.794926 6.376966 6.824059e-05 0.006755819
## 3 37/18845 0.08108108 31.183122 35.988491 9.383084 1.206178e-04 0.007960774
## 4 123/18845 0.03252033 12.507052 13.951074 6.537209 2.917713e-04 0.012596885
## 5 11/18845 0.18181818 69.925788 88.827423 11.675048 3.588643e-04 0.012596885
## 6 132/18845 0.03030303 11.654298 12.963889 6.271797 3.817238e-04 0.012596885
## qvalue geneID geneID_symbol
## 1 0.004417681 7516/5888/29128 XRCC2/RAD51/UHRF1
## 2 0.004417681 9156/5427/7516/5888/641/54821 EXO1/POLE2/XRCC2/RAD51/BLM/ERCC6L
## 3 0.005205610 7516/641/63967 XRCC2/BLM/CLSPN
## 4 0.008237198 5888/641/8318/55388 RAD51/BLM/CDC45/MCM10
## 5 0.008237198 57214/3949 CEMIP/LDLR
## 6 0.008237198 7516/5888/641/54821 XRCC2/RAD51/BLM/ERCC6L
## Count
## 1 3
## 2 6
## 3 3
## 4 4
## 5 2
## 6 4
table(gora$Cluster)##
## down.time2 up.time1 up.time2
## 15 3 158.3.3 Add background genes
With the help of transId(), background genes can be of a different type from the input IDs.
For example, the input genes are symbols while the background genes are Entrez.
## [1] "948" "1638" "158471" "10610" "6947" "100133941"8.3.4 [GO specific] Simplify GO ORA result
To improve GO term simplifying performance, genekitr extracts species-specific GO term information from Bioconductor organism-level annotation packages for fifteen species, including human, mouse, rat, fly, thale cress, yeast, zebrafish, worm, bovine, pig, chicken, mosquito, dog, frog, and chimpanzee.
Then genekitr obtains the ancestors of GO terms and their relations from GO.db
To calculate semantic similarity for GO BP, CC and MF, genekitr utilizes five algorithms (‘Resnik’, ‘Lin’, ‘Jiang’, ‘Rel’ and ‘Wang’) from GOSemSim
## [1] 43 13## [1] 21 13Let’s look at which terms are considered redundant:
## [1] "GO:0031625" "GO:0140666" "GO:0140852" "GO:0140036" "GO:0001091"
## [6] "GO:0001223"
Figure 8.1: Term GO:0044389
The first record GO:0044389 has the child term GO:0031625, so we only keep the child term.
8.4 Comparison with other tools
8.4.1 genekitr::genORA vs clusterProfiler::enrichGO/enrichKEGG/enrichDO...
The function genORA does not require you to load the orgdb in advance and does not require specifying the input ID type.
Next, let’s take GO analysis as an example:
The input IDs are a mixture of official symbols and aliases. All five genes are also very popular in biomedical research.
check_gene <- c("BCC7", "PDL1", "PD1", "TET2", "BRCA1")
# Genekitr
start = Sys.time()
gs <- geneset::getGO(org = "human",ont = "bp")
genekitr_res <- genekitr::genORA(check_gene, gs, p_cutoff = 0.01, q_cutoff = 0.01)
end = Sys.time()
(end-start)## Time difference of 33.77259 secs
# clusterProfiler
start = Sys.time()
require(org.Hs.eg.db)
clustp_res <- clusterProfiler::enrichGO(check_gene,
OrgDb = "org.Hs.eg.db",
keyType = "SYMBOL",
pvalueCutoff = 0.01,
qvalueCutoff = 0.01,
ont = "BP")
end = Sys.time()
(end-start)## Time difference of 1.073757 minsThen check the result:
class(genekitr_res)## [1] "data.frame"
class(clustp_res)## [1] "enrichResult"
## attr(,"package")
## [1] "enrichit"
isS4(clustp_res)## [1] TRUE
# we need to convert the S4 object to a data frame to view it
clustp_res2 <- as.data.frame(clustp_res)
dim(genekitr_res)## [1] 69 14
dim(clustp_res2)## [1] 0 13
head(genekitr_res,5)## Hs_BP_ID
## GO:0002835 GO:0002835
## GO:0002838 GO:0002838
## GO:0070234 GO:0070234
## GO:0070230 GO:0070230
## GO:0002840 GO:0002840
## Description
## GO:0002835 negative regulation of response to tumor cell
## GO:0002838 negative regulation of immune response to tumor cell
## GO:0070234 positive regulation of T cell apoptotic process
## GO:0070230 positive regulation of lymphocyte apoptotic process
## GO:0002840 regulation of T cell mediated immune response to tumor cell
## GeneRatio BgRatio RichFactor FoldEnrich oddsRatio zScore
## GO:0002835 0.4 11/19073 0.1818182 693.5636 1411.7778 37.20509
## GO:0002838 0.4 11/19073 0.1818182 693.5636 1411.7778 37.20509
## GO:0070234 0.4 11/19073 0.1818182 693.5636 1411.7778 37.20509
## GO:0070230 0.4 13/19073 0.1538462 586.8615 1154.9697 34.21651
## GO:0002840 0.4 20/19073 0.1000000 381.4600 705.5556 27.56594
## pvalue p.adjust qvalue geneID geneID_symbol Count
## GO:0002835 3.021116e-06 0.0004249704 0.0001038840 PD1/PDL1 PDCD1/CD274 2
## GO:0002838 3.021116e-06 0.0004249704 0.0001038840 PD1/PDL1 PDCD1/CD274 2
## GO:0070234 3.021116e-06 0.0004249704 0.0001038840 BCC7/PDL1 TP53/CD274 2
## GO:0070230 4.283594e-06 0.0004519191 0.0001104716 BCC7/PDL1 TP53/CD274 2
## GO:0002840 1.042673e-05 0.0007986270 0.0001952244 PD1/PDL1 PDCD1/CD274 2
head(clustp_res,5)## [1] ID Description GeneRatio BgRatio RichFactor
## [6] FoldEnrichment oddsRatio zScore pvalue p.adjust
## [11] qvalue geneID Count
## <0 rows> (or 0-length row.names)-
clustp_resonly contains two genes:TET2andBRCA1 -
genekitr_resincludes all five genes and provides mapping for aliases
Let’s also check whether some terms are found only in clustp_res (this depends on the GO data version; with older releases there were a few):
clustp_res$ID[!clustp_res$ID%in%genekitr_res$Hs_BP_ID ]## character(0)Question 1: Why could some terms be missing from the genORA result?
To get the adjusted p-value, the critical parameter is the number of statistical tests performed.
For a simple explanation, adj.p = f * p (f>1) and f is related to the testing method and the number of tests. genORA recognizes five genes, so we need to perform more tests (for example, five tests), while enrichGO only includes two genes, so the number of tests will be lower. This directly reduces the adj.p in the enrichGO result, allowing it to remain under the threshold of 0.01.
So, if such terms exist, a slightly higher cutoff (e.g. 0.05) lets genekitr retrieve them as well:
genekitr_res2 <- genekitr::genORA(check_gene, gs, p_cutoff = 0.05, q_cutoff = 0.05)
table(clustp_res$ID%in%genekitr_res2$Hs_BP_ID)## < table of extent 0 >Question 2: Why does enrichGO only recognize two genes?
It’s hard to imagine that such essential genes (“BCC7”, “PD1” or “PDL1”) would have no enriched terms at all.
The core reason is that these three genes are actually aliases, while enrichGO sets the argument keyType = "SYMBOL", which causes enrichGO to only accept official symbols. This function is based on an organism annotation package (e.g. org.Hs.eg.db) and it will accept different types of data from org.Hs.eg.db:
## [1] "ACCNUM" "ALIAS" "ENSEMBL" "ENSEMBLPROT" "ENSEMBLTRANS"
## [6] "ENTREZID" "ENZYME" "EVIDENCE" "EVIDENCEALL" "GENENAME"
## [11] "GENETYPE" "GO" "GOALL" "IPI" "MAP"
## [16] "OMIM" "ONTOLOGY" "ONTOLOGYALL" "PATH" "PFAM"
## [21] "PMID" "PROSITE" "REFSEQ" "SYMBOL" "UCSCKG"
## [26] "UNIPROT"Only by specifying keyType = 'ALIAS' can the program recognize gene aliases; however, it still cannot distinguish a mixture of symbols and aliases.
