Contents

1 The geneList dataset

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.

2 Over-representation analysis

Over-representation test1 is a widely used approach to identify biological themes. DOSE implements hypergeometric model to assess whether the number of selected genes associated with disease is larger than expected.

To determine whether any terms annotate a specified list of genes at frequency greater than that would be expected by chance, DOSE calculates a p-value using the hypergeometric distribution:

\(p = 1 - \displaystyle\sum_{i = 0}^{k-1}\frac{{M \choose i}{{N-M} \choose {n-i}}} {{N \choose n}}\)

In this equation, N is the total number of genes in the background distribution, M is the number of genes within that distribution that are annotated (either directly or indirectly) to the node of interest, n is the size of the list of genes of interest and k is the number of genes within that list which are annotated to the node. The background distribution by default is all the genes that have annotation. User can set the background via universe parameter.

P-values were adjusted for multiple comparison, and q-values were also calculated for FDR control.

2.1 enrichDO function

In the following example, we selected fold change above 1 as the differential genes and analyzing their disease association.

library(DOSE)
data(geneList)
gene <- names(geneList)[abs(geneList) > 1.5]
head(gene)
## [1] "4312"  "8318"  "10874" "55143" "55388" "991"
x <- enrichDO(gene          = gene,
              ont           = "DO", 
              pvalueCutoff  = 0.05,
              pAdjustMethod = "BH",
              universe      = names(geneList), 
              minGSSize     = 5,
              maxGSSize     = 500,
              qvalueCutoff  = 0.05,
              readable      = FALSE)
head(x)
##                    ID                    Description GeneRatio  BgRatio
## DOID:170     DOID:170         endocrine gland cancer    48/331 472/6268
## DOID:10283 DOID:10283                prostate cancer    40/331 394/6268
## DOID:3459   DOID:3459               breast carcinoma    37/331 357/6268
## DOID:3856   DOID:3856 male reproductive organ cancer    40/331 404/6268
## DOID:824     DOID:824                  periodontitis    16/331 109/6268
## DOID:3905   DOID:3905                 lung carcinoma    43/331 465/6268
##                  pvalue    p.adjust      qvalue
## DOID:170   5.662129e-06 0.004784499 0.003826407
## DOID:10283 3.859157e-05 0.013921739 0.011133923
## DOID:3459  4.942629e-05 0.013921739 0.011133923
## DOID:3856  6.821467e-05 0.014410349 0.011524689
## DOID:824   1.699304e-04 0.018859464 0.015082872
## DOID:3905  1.749754e-04 0.018859464 0.015082872
##                                                                                                                                                                                                                                                   geneID
## DOID:170   10874/7153/1381/6241/11065/10232/332/6286/2146/10112/891/9232/4171/993/5347/4318/3576/1515/4821/8836/3159/7980/5888/333/898/9768/4288/3551/2152/9590/185/7043/3357/2952/5327/3667/1634/1287/4582/7122/3479/4680/6424/80310/652/8839/9547/1524
## DOID:10283                                          4312/6280/6279/597/3627/332/6286/2146/4321/4521/891/5347/4102/4318/701/3576/79852/10321/6352/4288/3551/2152/247/2952/3487/367/3667/4128/4582/563/3679/4117/7031/3479/6424/10451/80310/652/4036/10551
## DOID:3459                                                          4312/6280/6279/7153/4751/890/4085/332/6286/6790/891/9232/10855/4171/5347/4318/701/2633/3576/9636/898/8792/4288/2952/4982/4128/4582/7031/3479/771/4250/2066/3169/10647/5304/5241/10551
## DOID:3856                                           4312/6280/6279/597/3627/332/6286/2146/4321/4521/891/5347/4102/4318/701/3576/79852/10321/6352/4288/3551/2152/247/2952/3487/367/3667/4128/4582/563/3679/4117/7031/3479/6424/10451/80310/652/4036/10551
## DOID:824                                                                                                                                                                   4312/6279/820/7850/4321/3595/4318/4069/3576/1493/6352/8842/185/2952/5327/4982
## DOID:3905                          4312/6280/2305/9133/6279/7153/6278/6241/55165/11065/8140/10232/332/6286/3002/9212/4521/891/4171/9928/8061/4318/3576/1978/1894/7980/7083/898/6352/8842/4288/2152/2697/2952/3572/4582/7049/563/3479/1846/3117/2532/2922
##            Count
## DOID:170      48
## DOID:10283    40
## DOID:3459     37
## DOID:3856     40
## DOID:824      16
## DOID:3905     43

The enrichDO function requires an entrezgene ID vector as input, mostly is the differential gene list of gene expression profile studies. If user needs to convert other gene ID type to entrezgene ID, we recommend using bitr function provided by clusterProfiler.

The ont parameter can be “DO” or “DOLite”, DOLite2 was constructed to aggregate the redundant DO terms. The DOLite data is not updated, we recommend user use ont="DO". pvalueCutoff setting the cutoff value of p value and p value adjust; pAdjustMethod setting the p value correction methods, include the Bonferroni correction (“bonferroni”), Holm (“holm”), Hochberg (“hochberg”), Hommel (“hommel”), Benjamini & Hochberg (“BH”) and Benjamini & Yekutieli (“BY”) while qvalueCutoff is used to control q-values.

The universe setting the background gene universe for testing. If user do not explicitly setting this parameter, enrichDO will set the universe to all human genes that have DO annotation.

The minGSSize (and maxGSSize) indicates that only those DO terms that have more than minGSSize (and less than maxGSSize) genes annotated will be tested.

The readable is a logical parameter, indicates whether the entrezgene IDs will mapping to gene symbols or not.

We also implement setReadable function that helps the user to convert entrezgene IDs to gene symbols.

x <- setReadable(x, 'org.Hs.eg.db')
head(x)
##                    ID                    Description GeneRatio  BgRatio
## DOID:170     DOID:170         endocrine gland cancer    48/331 472/6268
## DOID:10283 DOID:10283                prostate cancer    40/331 394/6268
## DOID:3459   DOID:3459               breast carcinoma    37/331 357/6268
## DOID:3856   DOID:3856 male reproductive organ cancer    40/331 404/6268
## DOID:824     DOID:824                  periodontitis    16/331 109/6268
## DOID:3905   DOID:3905                 lung carcinoma    43/331 465/6268
##                  pvalue    p.adjust      qvalue
## DOID:170   5.662129e-06 0.004784499 0.003826407
## DOID:10283 3.859157e-05 0.013921739 0.011133923
## DOID:3459  4.942629e-05 0.013921739 0.011133923
## DOID:3856  6.821467e-05 0.014410349 0.011524689
## DOID:824   1.699304e-04 0.018859464 0.015082872
## DOID:3905  1.749754e-04 0.018859464 0.015082872
##                                                                                                                                                                                                                                                                                             geneID
## DOID:170   NMU/TOP2A/CRABP1/RRM2/UBE2C/MSLN/BIRC5/S100P/EZH2/KIF20A/CCNB1/PTTG1/MCM2/CDC25A/PLK1/MMP9/CXCL8/CTSV/NKX2-2/GGH/HMGA1/TFPI2/RAD51/APLP1/CCNE1/KIAA0101/MKI67/IKBKB/F3/AKAP12/AGTR1/TGFB3/HTR2B/GSTT1/PLAT/IRS1/DCN/COL4A5/MUC1/CLDN5/IGF1/CEACAM6/SFRP4/PDGFD/BMP4/WISP2/CXCL14/CX3CR1
## DOID:10283                                                      MMP1/S100A9/S100A8/BCL2A1/CXCL10/BIRC5/S100P/EZH2/MMP12/NUDT1/CCNB1/PLK1/MAGEA3/MMP9/BUB1B/CXCL8/EPHX3/CRISP3/CCL5/MKI67/IKBKB/F3/ALOX15B/GSTT1/IGFBP4/AR/IRS1/MAOA/MUC1/AZGP1/ITGA7/MAK/TFF1/IGF1/SFRP4/VAV3/PDGFD/BMP4/LRP2/AGR2
## DOID:3459                                                              MMP1/S100A9/S100A8/TOP2A/NEK2/CCNA2/MAD2L1/BIRC5/S100P/AURKA/CCNB1/PTTG1/HPSE/MCM2/PLK1/MMP9/BUB1B/GBP1/CXCL8/ISG15/CCNE1/TNFRSF11A/MKI67/GSTT1/TNFRSF11B/MAOA/MUC1/TFF1/IGF1/CA12/SCGB2A2/ERBB4/FOXA1/SCGB1D2/PIP/PGR/AGR2
## DOID:3856                                                       MMP1/S100A9/S100A8/BCL2A1/CXCL10/BIRC5/S100P/EZH2/MMP12/NUDT1/CCNB1/PLK1/MAGEA3/MMP9/BUB1B/CXCL8/EPHX3/CRISP3/CCL5/MKI67/IKBKB/F3/ALOX15B/GSTT1/IGFBP4/AR/IRS1/MAOA/MUC1/AZGP1/ITGA7/MAK/TFF1/IGF1/SFRP4/VAV3/PDGFD/BMP4/LRP2/AGR2
## DOID:824                                                                                                                                                                                           MMP1/S100A8/CAMP/IL1R2/MMP12/IL12RB2/MMP9/LYZ/CXCL8/CTLA4/CCL5/PROM1/AGTR1/GSTT1/PLAT/TNFRSF11B
## DOID:3905                               MMP1/S100A9/FOXM1/CCNB2/S100A8/TOP2A/S100A7/RRM2/CEP55/UBE2C/SLC7A5/MSLN/BIRC5/S100P/GZMB/AURKB/NUDT1/CCNB1/MCM2/KIF14/FOSL1/MMP9/CXCL8/EIF4EBP1/ECT2/TFPI2/TK1/CCNE1/CCL5/PROM1/MKI67/F3/GJA1/GSTT1/IL6ST/MUC1/TGFBR3/AZGP1/IGF1/DUSP4/HLA-DQA1/ACKR1/GRP
##            Count
## DOID:170      48
## DOID:10283    40
## DOID:3459     37
## DOID:3856     40
## DOID:824      16
## DOID:3905     43

2.2 enrichNCG function

Network of Cancer Gene (NCG)3 is a manually curated repository of cancer genes. NCG release 5.0 (Aug. 2015) collects 1,571 cancer genes from 175 published studies. DOSE supports analyzing gene list and determine whether they are enriched in genes known to be mutated in a given cancer type.

gene2 <- names(geneList)[abs(geneList) < 3]
ncg <- enrichNCG(gene2)
head(ncg)
##                                        ID          Description GeneRatio
## soft_tissue_sarcomas soft_tissue_sarcomas soft_tissue_sarcomas   28/1172
## bladder                           bladder              bladder   61/1172
## glioma                             glioma               glioma   68/1172
##                      BgRatio       pvalue    p.adjust      qvalue
## soft_tissue_sarcomas 28/1571 0.0002517511 0.008056037 0.006360029
## bladder              67/1571 0.0005108168 0.008173069 0.006452423
## glioma               76/1571 0.0008511747 0.009079196 0.007167787
##                                                                                                                                                                                                                                                                                                                                                                  geneID
## soft_tissue_sarcomas                                                                                                                                                                                                       1029/999/6850/4914/4342/2185/55294/2041/4851/23512/2044/4058/5290/8726/4486/5297/5728/3815/2324/7403/5925/4763/1499/7157/5159/2045/3667/2066
## bladder                                            9700/2175/9603/1029/8997/688/1026/896/677/6256/55294/8085/4851/3265/1999/3845/8243/10605/8295/4854/5290/2033/4780/23224/23217/2064/23385/55252/10735/4853/387/288/30849/9794/7403/287/463/472/4297/2065/2262/8289/9611/5925/2068/4763/7157/2186/1387/3910/2261/7248/23037/23345/7832/79633/10628/22906/388/4036/3169
## glioma               4603/4609/1029/3418/8877/1019/7027/4613/1030/1956/1106/2264/3417/6597/4914/55359/896/894/2321/3954/5335/5781/8439/673/9444/4851/8087/2050/8493/3845/3482/667/56999/5290/2033/4233/577/5894/5156/80036/9407/3020/1021/5598/5728/8621/1828/63035/23592/8880/2260/54880/4916/2263/1639/90/546/8289/4763/7157/23152/5295/4602/595/2261/6938/4915/26137
##                      Count
## soft_tissue_sarcomas    28
## bladder                 61
## glioma                  68

2.3 enrichDGN and enrichDGNv functions

DisGeNET4 is an integrative and comprehensive resources of gene-disease associations from several public data sources and the literature. It contains gene-disease associations and snp-gene-disease associations.

The enrichment analysis of disease-gene associations is supported by the enrichDGN function and analysis of snp-gene-disease associations is supported by the enrichDGNv function.

dgn <- enrichDGN(gene)
head(dgn)
##                          ID                      Description GeneRatio
## umls:C1134719 umls:C1134719 Invasive Ductal Breast Carcinoma    28/476
## umls:C0032460 umls:C0032460        Polycystic Ovary Syndrome    38/476
## umls:C0206698 umls:C0206698               Cholangiocarcinoma    36/476
## umls:C0007138 umls:C0007138     Carcinoma, Transitional Cell    35/476
## umls:C0031099 umls:C0031099                    Periodontitis    28/476
## umls:C0005695 umls:C0005695                 Bladder Neoplasm    36/476
##                 BgRatio       pvalue     p.adjust       qvalue
## umls:C1134719 231/17381 4.312190e-11 1.225524e-07 9.164539e-08
## umls:C0032460 434/17381 2.819624e-10 3.521620e-07 2.633487e-07
## umls:C0206698 399/17381 3.717403e-10 3.521620e-07 2.633487e-07
## umls:C0007138 389/17381 7.093837e-10 5.040171e-07 3.769068e-07
## umls:C0031099 270/17381 1.634417e-09 9.290027e-07 6.947133e-07
## umls:C0005695 442/17381 5.871618e-09 2.781190e-06 2.079789e-06
##                                                                                                                                                                                                         geneID
## umls:C1134719                                                 9133/7153/6241/55165/11065/51203/22974/4751/5080/332/2568/3902/6790/891/24137/9232/10855/79801/4318/55635/5888/1493/9768/3070/4288/367/4582/5241
## umls:C0032460 4312/6280/6279/7153/259266/6241/55165/55872/4085/6286/7272/366/891/4171/7941/1164/3161/4603/990/29127/4318/53335/3294/3070/2952/5327/367/3667/4582/563/27324/3479/114899/9370/2167/652/5346/5241
## umls:C0206698             4312/2305/55872/4751/8140/10635/10232/5918/332/6286/2146/4521/891/10855/2921/7941/1164/4318/3576/1978/79852/8842/4485/214/65982/6863/1036/6935/4128/3572/4582/7031/7166/4680/80310/9
## umls:C0007138                       4312/991/6280/6241/55165/10460/6373/8140/890/10232/4085/332/6286/2146/4171/1033/6364/5347/4318/3576/8836/9700/898/4288/2952/367/8382/2947/3479/9338/23158/2167/2066/2625/9
## umls:C0031099                                                       4312/6279/3669/820/7850/332/4321/6364/3595/4318/3576/3898/8792/1493/4485/10472/185/6863/2205/2952/5327/4982/23261/2200/3572/2006/1308/2625
## umls:C0005695                   4312/10874/6280/3868/6279/597/7153/6241/9582/10460/4085/5080/332/2146/6790/10855/4171/5347/4318/3576/8836/9636/9700/898/4288/214/2952/367/2947/4582/3479/6424/9338/2066/1580/9
##               Count
## umls:C1134719    28
## umls:C0032460    38
## umls:C0206698    36
## umls:C0007138    35
## umls:C0031099    28
## umls:C0005695    36
snp <- c("rs1401296", "rs9315050", "rs5498", "rs1524668", "rs147377392",
         "rs841", "rs909253", "rs7193343", "rs3918232", "rs3760396",  
         "rs2231137", "rs10947803", "rs17222919", "rs386602276", "rs11053646", 
         "rs1805192", "rs139564723", "rs2230806", "rs20417", "rs966221")
dgnv <- enrichDGNv(snp)
head(dgnv)
##                          ID                       Description GeneRatio
## umls:C3272363 umls:C3272363 Ischemic Cerebrovascular Accident     20/20
## umls:C0948008 umls:C0948008                   Ischemic stroke     20/20
## umls:C0038454 umls:C0038454          Cerebrovascular accident      7/20
## umls:C0027051 umls:C0027051             Myocardial Infarction      6/20
## umls:C0010054 umls:C0010054         Coronary Arteriosclerosis      6/20
## umls:C0010068 umls:C0010068            Coronary heart disease      6/20
##                 BgRatio       pvalue     p.adjust       qvalue
## umls:C3272363 141/46589 1.014503e-51 1.379725e-49 1.922217e-50
## umls:C0948008 148/46589 2.867870e-51 1.950151e-49 2.716929e-50
## umls:C0038454 243/46589 7.045680e-12 3.194042e-10 4.449903e-11
## umls:C0027051 163/46589 6.222154e-11 1.889883e-09 2.632964e-10
## umls:C0010054 166/46589 6.948100e-11 1.889883e-09 2.632964e-10
## umls:C0010068 314/46589 3.198889e-09 7.250815e-08 1.010175e-08
##                                                                                                                                                                                                              geneID
## umls:C3272363 rs1401296/rs9315050/rs5498/rs1524668/rs147377392/rs841/rs909253/rs7193343/rs3918232/rs3760396/rs2231137/rs10947803/rs17222919/rs386602276/rs11053646/rs1805192/rs139564723/rs2230806/rs20417/rs966221
## umls:C0948008 rs1401296/rs9315050/rs5498/rs1524668/rs147377392/rs841/rs909253/rs7193343/rs3918232/rs3760396/rs2231137/rs10947803/rs17222919/rs386602276/rs11053646/rs1805192/rs139564723/rs2230806/rs20417/rs966221
## umls:C0038454                                                                                                                              rs1524668/rs147377392/rs2231137/rs10947803/rs386602276/rs2230806/rs20417
## umls:C0027051                                                                                                                                              rs5498/rs147377392/rs909253/rs11053646/rs1805192/rs20417
## umls:C0010054                                                                                                                                             rs5498/rs147377392/rs11053646/rs1805192/rs2230806/rs20417
## umls:C0010068                                                                                                                                             rs5498/rs147377392/rs11053646/rs1805192/rs2230806/rs20417
##               Count
## umls:C3272363    20
## umls:C0948008    20
## umls:C0038454     7
## umls:C0027051     6
## umls:C0010054     6
## umls:C0010068     6

3 Visualze enrichment result

To help interpreting enrichment result, we implemented barplot, dotplot, cnetplot (category-gene-network) upsetplot and enrichMap for visualization.

3.1 barplot

barplot(x, showCategory=10)

3.2 dotplot

dotplot is a good alternative to barplot.

dotplot(x)

3.3 cnetplot

In order to consider the potentially biological complexities in which a gene may belong to multiple annotation categories, we developed cnetplot function to extract the complex association between genes and diseases.

cnetplot(x, categorySize="pvalue", foldChange=geneList)

3.4 upsetplot

upsetplot is an alternative to cnetplot for visualizing the complex association between genes and diseases.

upsetplot(x)

3.5 enrichMap

Enrichment Map can be visualized by enrichMap function. It’s designed to summarize enriched result.

enrichMap(x)

4 Disease analysis of NGS data

Disease analysis using NGS data (eg, RNA-Seq and ChIP-Seq) can be performed by linking coding and non-coding regions to coding genes via ChIPseeker package, which can annotates genomic regions to their nearest genes, host genes, and flanking genes respectivly. In addtion, it provides a function, seq2gene, that simultaneously considering host genes, promoter region and flanking gene from intergenic region that may under control via cis-regulation. This function maps genomic regions to genes in a many-to-many manner and facilitate functional analysis. For more details, please refer to ChIPseeker5.

5 Disease association comparison

We have developed an R package clusterProfiler6 for comparing biological themes among gene clusters. DOSE works fine with clusterProfiler and can compare biological themes at disease perspective.

library(clusterProfiler)
data(gcSample)
cdo <- compareCluster(gcSample, fun="enrichDO")
plot(cdo)

6 Other enrichment analysis tools

We provide enrichment analysis in clusterProfiler6 for GO, KEGG, MeSH, DAVID, Molecular Signatures Database and others (user’s annotation), and Reactome pathway enrichment analysis in ReactomePA7 package. Both hypergeometric test and GSEA are supported.

References

1. Boyle, E. I. et al. GO::TermFinder–open source software for accessing gene ontology information and finding significantly enriched gene ontology terms associated with a list of genes. Bioinformatics (Oxford, England) 20, 3710–3715 (2004).

2. Du, P. et al. From disease ontology to disease-ontology lite: Statistical methods to adapt a general-purpose ontology for the test of gene-ontology associations. Bioinformatics 25, i63–i68 (2009).

3. A., O., D., G. M., M., T. P. & C., F. D. NCG 5.0: Updates of a manually curated repository of cancer genes and associated properties from cancer mutational screenings. Nucleic Acids Research 44, D992–D999 (2016).

4. Janet, P. et al. DisGeNET: A discovery platform for the dynamical exploration of human diseases and their genes. Database 2015, bav028 (2015).

5. Yu, G., Wang, L.-G. & He, Q.-Y. ChIPseeker: An r/bioconductor package for chip peak annotation, comparison and visualization. Bioinformatics 31, 2382–2383 (2015).

6. Yu, G., Wang, L.-G., Han, Y. & He, Q.-Y. clusterProfiler: an r package for comparing biological themes among gene clusters. OMICS: A Journal of Integrative Biology 16, 284–287 (2012).

7. Yu, G. & He, Q.-Y. ReactomePA: An r/bioconductor package for reactome pathway analysis and visualization. Mol. BioSyst. 12, 477–479 (2016).