GSEA analysis requires a ranked gene list, which contains three features:
numeric vector: fold change or other type of numerical variable
named vector: every number has a name, the corresponding gene ID
sorted vector: number should be sorted in decreasing order
If you import your data from a csv file, the file should contains two columns, one for gene ID (no duplicated ID allowed) and another one for fold change. You can prepare your own geneList via the following command:
d =read.csv(your_csv_file)## assume 1st column is ID## 2nd column is FC## feature 1: numeric vectorgeneList = d[,2]## feature 2: named vectornames(geneList) =as.character(d[,1])## feature 3: decreasing ordergeneList =sort(geneList, decreasing =TRUE)
By default, all the visualization methods provided by enrichplot display most significant pathways. If users are interested to show some specific pathways (e.g. excluding some unimportant pathways among the top categories), users can pass a vector of selected pathways to the showCategory parameter in dotplot(), barplot(), treeplot(), cnetplot() and emapplot() etc.
See Figure 18.1 for a side-by-side comparison of default top categories versus user-selected pathways.
library(clusterProfiler)library(enrichplot)data(geneList, package='DOSE')de <-names(geneList)[abs(geneList)>1]x <-enrichKEGG(de)## show top 10 most significant pathways and want to exclude the second one## dotplot(x, showCategory = x$Description[1:10][-2])set.seed(2020-10-27)selected_pathways <-sample(x$Description, 5)selected_pathways
[1] "Protein digestion and absorption"
[2] "Viral protein interaction with cytokine and cytokine receptor"
[3] "Cytokine-cytokine receptor interaction"
[4] "DNA replication"
[5] "Cytoskeleton in muscle cells"
Most of the functions in enrichplot can automatically split long labels across multiple lines. Users can passed a line width to the label_format parameter (default is 30). It also supports user defined function to format label strings.
See Figure 18.2 for examples of wrapping long axis labels using numeric width and a custom labeller.
Figure 18.2: Wrap long axis labels. Passing a numeric value to specify string width (A), a user specifiable labeller function (B).
The label_format option works with barplot(), dotplot(), heatplot(), treeplot and ridgeplot().
18.6 Why many genes are not annotated in KEGG?
Users often find that a significant portion of their gene list is dropped in the KEGG enrichment analysis. For example, a user reported that only 281 genes were retained from a list of over 700 genes.
This is not a software issue but rather a reflection of the database coverage. KEGG pathways mainly focus on metabolic and signaling pathways, and many genes are not annotated in any KEGG pathway.
We can verify the annotation coverage using bitr_kegg():
# gene is a vector of gene IDs (e.g., Entrez IDs)kk <-bitr_kegg(geneID = gene, fromType='ncbi-geneid', toType='Path', organism='hsa')
The warning message will indicate the percentage of genes that failed to map.
To manually verify a specific gene, you can visit the KEGG website using a URL constructed with the organism code and gene ID, for example: http://www.genome.jp/dbget-bin/www_bget?hsa:100506775. If the gene page does not list any pathway, it means the gene has no KEGG pathway annotation.
