Step 1: Collect and store pathways from KEGG API

First we will use the get_kegg function to access the KEGG REST API and download the data sets required to perform our downstream analysis. This function takes two arguments. The first, ‘species’ is required. Currently, diffEnrich supports three species, and the argument is a character string using the KEGG code https://www.pnas.org/doi/10.1073/pnas.0806162105: Homo sapiens (human), use ‘hsa’; Mus musculus (mouse), use ‘mmu’; and Rattus norvegicus (rat), use ‘rno’. The second, ‘path’ is also passed as a character string, and is the path of the directory in which the user would like to write the data sets downloaded from the KEGG REST API. If the user does not provide a path, the data sets will be automatically written to the current working directory using the here::here() functionality. These data sets will be tab delimited files with a name describing the data, and for reproducibility, the date they were generated and the version of KEGG when the API was accessed. In addition to these flat files, get_kegg will also create a named list in R with the three relevant KEGG data sets. The names of this list will describe the data set. For a detailed description of list elements use ?get_kegg.

## run get_kegg() using rat
kegg_rno <- get_kegg('rno')
#> 3 data sets will be written as tab delimited text files
#> File location: /private/var/folders/ms/g80_v_b157b6tryj42yxgfjm0000gq/T/Rtmpzqq8Ou/Rbuildf55625004d35/diffEnrich
#> Kegg Release: Release_102.0+_06-27_Jun_22

Here are examples of the output files:

ncbi_to_kegg2019-10-17Release_92.0+_10-17_Oct_19.txt
kegg_to_pathway2019-10-17Release_92.0+_10-17_Oct_19.txt
pathway_to_species2019-10-17Release_92.0+_10-17_Oct_19.txt

Note: Because it is assumed that the user might want to use the data sets generated by get_kegg, it is careful not to overwrite data sets with exact names. get_kegg checks the path provided for data sets generated ‘same-day/same-version’, and if it finds even one of the three, it will not re-write any of the data sets. It will still however, let the user know it is not writing out new data sets and still generate the named list object. Users can generate ‘same-day/same-version’ data sets in different directories if they so choose.

## run get_kegg() using rat
kegg_rno <- get_kegg('rno')
#> These files already exist in your working directory. New files will not be generated.
#> Kegg Release: Release_102.0+_06-27_Jun_22

## run get_kegg() using rat
kegg_rno <- get_kegg(read = TRUE, path = "/path/to/files", date = "2019-10-17", release = "92")

Here is an example of the output:

Reading in the following files:
ncbi_to_kegg2019-10-17Release_92.0+_10-17_Oct_19.txt
kegg_to_pathway2019-10-17Release_92.0+_10-17_Oct_19.txt
pathway_to_species2019-10-17Release_92.0+_10-17_Oct_19.txt
File location: ~/Desktop

Step 2: Perform enrichment analysis of individual gene sets

In this step we will use the pathEnrich function to identify KEGG pathways that are enriched (i.e. over-represented) based on a gene list of interest. User gene lists must be character vectors and be formatted as ENTREZ gene IDs. The clusterProfiler package offers a nice function (bitr) that maps gene symbols and Ensembl IDs to ENTREZ gene IDs, and an example can be seen in their vignette.

## View sample gene lists from package data
head(geneLists$list1)
#> [1] "361692"    "293654"    "293655"    "500974"    "100361529" "171434"
head(geneLists$list2)
#> [1] "315547" "315548" "315549" "315550" "50938"  "58856"

This function may not always use the complete list of genes provided by the user. Specifically, it will only use the genes from the list provided that are also in the most current species list pulled from the KEGG REST API using get_kegg, or from the older KEGG data loaded by the user from a previous get_kegg call. Users can also decide which KEGG pathways should be tested based on how many genes from their gene list are contained in the KEGG pathway. Users can set this parameter by changing the ‘N’ argument. The default is N = 2. The pathEnrich function should be run once for the genes of interest in list 1 and once for the genes of interest in list2. Each pathEnrich call generates a data frame summarizing the results of an enrichment analysis in which a Fisher’s Exact test is used to identify which KEGG pathways are enriched for the user’s list of genes compared to all genes annotated to a KEGG pathway. Users can choose a multiple correction option from those supported by stats::p.adjust. The default is the False Discovery Rate ( Benjamini and Hochberg, 1995), and the default threshold to reach significance is 0.05.

# run pathEnrich using kegg_rno
## List 1
list1_pe <- pathEnrich(gk_obj = kegg, gene_list = geneLists$list1, cutoff = 0.05, N = 2)
## list2
list2_pe <- pathEnrich(gk_obj = kegg, gene_list = geneLists$list2, cutoff = 0.05, N = 2) 

pathEnrich generates a list object that contains a data frame with 9 columns described below as well as other metadata from the function call. Details also provided in pathEnrich documentation. S3 generic functions for print and summary are provided. The print function prints the results table as a tibble, and the summary function returns the number of pathways that reached statistical significance as well as their descriptions, the number of genes used from the KEGG data base, the KEGG species, and the method used for multiple testing correction.

summary(list1_pe)
#> 219 KEGG pathways were tested. 
#>  Only KEGG pathways that contained at least 2 genes from gene_list were tested. 
#>  KEGG pathway species: Rattus norvegicus (rat)
#>  8856 genes from gene_list were in the KEGG data pull. 
#>  p-value adjustment method: BH
#>  36 pathways reached statistical significance after multiple testing correction at a cutoff of 0.05. 
#>  
#> Significant pathways: 
#>  Tight junction
#> Yersinia infection
#> Colorectal cancer
#> Choline metabolism in cancer
#> Endometrial cancer
#> Endocytosis
#> Neurotrophin signaling pathway
#> Thermogenesis
#> Oocyte meiosis
#> VEGF signaling pathway
#> Thyroid hormone signaling pathway
#> Hippo signaling pathway
#> T cell receptor signaling pathway
#> Apoptosis
#> Hepatocellular carcinoma
#> MAPK signaling pathway
#> Focal adhesion
#> Salmonella infection
#> Non-alcoholic fatty liver disease (NAFLD)
#> ErbB signaling pathway
#> Sphingolipid signaling pathway
#> Pancreatic cancer
#> Progesterone-mediated oocyte maturation
#> Alzheimer disease
#> Endocrine resistance
#> Adrenergic signaling in cardiomyocytes
#> IL-17 signaling pathway
#> Chronic myeloid leukemia
#> Dopaminergic synapse
#> Prostate cancer
#> EGFR tyrosine kinase inhibitor resistance
#> Hepatitis C
#> Ras signaling pathway
#> Acute myeloid leukemia
#> Insulin signaling pathway
#> Fc epsilon RI signaling pathway

Step 3: Identify differentially enriched KEGG pathways

The diffEnrich function will merge two results from the pathEnrich calls generated above. Specifically, the data frame ‘list1_pe’ and the data frame ‘list2_pe’ will be merged by the following columns: “KEGG_PATHWAY_ID”, “KEGG_PATHWAY_description”, “KEGG_PATHWAY_cnt”, “KEGG_DATABASE_cnt”. This merged data set will then be used to perform differential enrichment using the same method and p-value calculation as described above. Users do have the option of choosing a method for multiple testing adjustment. Users can choose from those supported by stats::p.adjust. The default is the False Discovery Rate ( Benjamini and Hochberg, 1995). KEGG pathways that do not contain any genes from either gene list (i.e., list1_pe\(enrich_table\)KEGG_PATHWAY_in_list for ‘rno04530’ = 0 AND list2_pe\(enrich_table\)KEGG_PATHWAY_in_list for ‘rno04530’ = 0) will be removed as these cannot be tested. If this is the case a warning will be printed that tells the user how many pathways were removed. This can be avoided by setting the ‘N’ parameter to a value > 0 in the pathEnrich calls.

## Perform differential enrichment
diff_enrich <- diffEnrich(list1_pe = list1_pe, list2_pe = list2_pe, method = 'none', cutoff = 0.05)

The result of the diffEnrich call is a list object that contains a data frame with the estimated odds ratio generated by the Fisher’s Exact test and the associated p-value as well as other metadata from the function call. S3 generic functions for print and summary are provided. The print function prints the results table as a tibble, and the summary function returns the number of pathways that reached statistical significance as well as their descriptions, the number of genes used from the KEGG database, the KEGG species, the number of pathways that were shared (i.e. had at least 1 gene from each gene list present in the pathway based on what the user chose for N in pathEnrich) between the gene lists and the method used for multiple testing correction.

summary(diff_enrich)
#> 219 KEGG pathways were shared between gene lists and were tested. 
#>  KEGG pathway species: Rattus norvegicus (rat)
#>  8856 genes from gene_list were in the KEGG data pull. 
#>  p-value adjustment method: none
#>  19 pathways reached statistical significance after multiple testing correction at a cutoff of 0.05. 
#>  
#> Significant pathways: 
#>  Tight junction
#> Yersinia infection
#> Colorectal cancer
#> Endometrial cancer
#> T cell receptor signaling pathway
#> IL-17 signaling pathway
#> Salmonella infection
#> Choline metabolism in cancer
#> Endocytosis
#> VEGF signaling pathway
#> Oocyte meiosis
#> Thermogenesis
#> Hippo signaling pathway
#> Neurotrophin signaling pathway
#> Apoptosis
#> Hepatocellular carcinoma
#> Fc epsilon RI signaling pathway
#> Thyroid hormone signaling pathway
#> Alzheimer disease

Step 4: Plot fold enrichment

plotFoldEnrichment generates a grouped bar plot using ggplot2 and the ggnewscale package. There are 3 arguments: 1) de_res is the dataframe generated from the diffEnrich function, 2) pval is the threshold for the adjusted p-value associated with differential enrichment that will filter which KEGG pathways to plot, and 3) after filtering based on pval N tells the function how many pathways to plot. It is important to make a note that the significance of the fold change is associated with the number of genes in the gene list. Notice that in this example the pathways in gene list 2 have smaller fold changes (shorter bars) than those in list 1, but that many of them are more significant (darker blue). This is because there are more genes in gene list 2 compared to gene list 1.

## Plot fold enrichment
plotFoldEnrichment(de_res = diff_enrich, pval = 0.05, N = 5)

Figure 2. Example of a differential enrichment graphic. KEGG pathways are plotted on the y-axis and fold enrichment is plotted on the x-axis. Each KEGG pathway has a bar depicting its fold enrichment in list 1 (red) and its fold enrichment in list 2 (blue). The transparency of the bars correspond to the unadjusted p-value for the pathway’s enrichment in the given list. The p-value presented as text to the right of each pair of bars is the adjusted p-value (user defined: default is FDR) associated with the differential enrichment of the pathway between the two lists, and the pathways are ordered from top to bottom by this p-value (i.e. smallest p-value on top of plot, and largest p-value on bottom of plot). The dotted line represents a fold enrichment of 1. Finally, the number of genes used for analysis from each gene list (recall that this number may not be the same as the number of genes in the user’s original list) are reported below their respective p-values in the legend.