39  Frequently asked questions

39.1 How to prepare your own geneList

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 vector
geneList = d[,2]

## feature 2: named vector
names(geneList) = as.character(d[,1])

## feature 3: decreasing order
geneList = sort(geneList, decreasing = TRUE)

The ranking must also be finite. Check for NA, NaN, or Inf before calling GSEA; enrichit stops rather than assigning a meaningful rank to a non-finite statistic:

stopifnot(is.numeric(geneList), !is.null(names(geneList)))
stopifnot(all(is.finite(geneList)))

39.1.1 Inspecting pathways with NA multilevel p-values

Some multilevel GSEA pathways can have NA p-values when the relevant null side is too unbalanced for the estimate. gsea_gson() excludes these rows by default. Set reportNA = TRUE when the goal is to inspect the complete tested set rather than only pathways with reportable p-values; retain the NA values as an estimation status, not as evidence of non-significance. Using pvalueCutoff = 1 is a convenient way to inspect all retained rows.

all_tested <- gsea_gson(
  geneList = geneList,
  gson = gson,
  reportNA = TRUE,
  pvalueCutoff = 1
)
sum(is.na(all_tested@result$pvalue))

39.2 No gene can be mapped

39.2.1 First, check that your identifiers are still identifiers

Before blaming the annotation database, look at the IDs themselves. A very common cause of No gene can be mapped is that the ID column passed through a spreadsheet on its way into R. Excel quietly turns long numbers into scientific notation, and when the column is too narrow or the cell is formatted as a number it writes out the displayed text rather than the value. A gene list then arrives looking like this:

25081
84604
1.01E+08     <- should be an Entrez ID
1E+08        <- should be an Entrez ID
24614

None of these can be mapped, and the failure is easy to misread as “clusterProfiler does not support my organism” — particularly when pasting the same list into DAVID seems to work, because DAVID re-parses the text differently.

Note that 1E+08 is not merely mis-formatted, it is lossy: the original identifier is gone and cannot be recovered. The fix is therefore not to repair the strings in R, but to re-export the IDs from their source.

A quick check in R — any ID that is not a plain integer is suspect:

ids <- c("25081", "84604", "1.01E+08", "1E+08", "24614")
ids[!grepl("^[0-9]+$", ids)]
[1] "1.01E+08" "1E+08"

How to avoid it. Export the ID column as text: in Excel, format the column as “Text” before pasting the IDs in, or import via “Data → From Text/CSV” and set the column type to text. Then read it into R as character, e.g. readr::read_tsv(file, col_types = readr::cols(ENTREZID = "c")). Never round-trip identifiers through a spreadsheet’s default number formatting.

39.2.2 Other causes

If the IDs are clean integers and the problem persists, the usual suspects are a mismatched keyType, an organism code that does not correspond to your IDs, or an outdated annotation source — see the two links above.

39.2.3 Do not confuse a failed download with a failed mapping

No gene can be mapped means the mapping step found nothing, which invites you to question your identifiers. But a download that never delivered usable annotation used to produce the same symptom, because an empty gene-set table simply maps nothing. Since clusterProfiler 4.21.2 the two are reported separately, so read the message before editing your gene list:

  • unreachable network — Failed to reach the KEGG REST API at '<url>'. The network is unreachable (connection failed, timed out or DNS lookup failed). ... Your IDs are fine; the machine could not reach KEGG (proxy, firewall, or a timeout on a large response).
  • HTTP error status — ... the request returned HTTP status <code>. The KEGG REST API may have changed or is temporarily unavailable. Either KEGG retired that endpoint (this is how a KEGG API change shows up) or the service is down; updating the package or using a local database is the fix.
  • empty response — ... the response was empty. The requested data may not exist in KEGG ...
  • empty annotation — No KEGG pathway annotation could be assembled for species '<x>': the downloaded gene/pathway table shares no ID with the pathway-name table. The download succeeded but the two halves did not line up, which again usually means the KEGG response format changed.

Only if none of those appear should you suspect the identifiers themselves.

For an unstable or restricted connection, build the annotation once and reuse it locally instead of downloading on every run:

library(clusterProfiler)
createKEGGdb::create_kegg_db("hsa") # builds a local KEGG.db package
x <- enrichKEGG(gene, organism = "hsa", use_internal_data = TRUE)

39.3 Showing specific pathways

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 39.1 for a side-by-side comparison of default top categories versus user-selected pathways.

enrichKEGG() and enrichPathway() in this chapter reach the remote KEGG and Reactome services, which can fail transiently, so the calls are wrapped in the retry() helper defined in _common.R; readers who prefer a single attempt can simply drop the wrapper in their own analysis.

library(clusterProfiler)
library(enrichplot)
data(geneList, package='DOSE')
de <- names(geneList)[abs(geneList)>1]

x <- retry(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"                                 
p1 <- dotplot(x, showCategory = 10, font.size=14)
p2 <- dotplot(x, showCategory = selected_pathways, font.size=14)


aplot::plot_list(p1, p2, tag_levels = "A")
Figure 39.1: Showing specific pathways. Top ten most significant pathways (A), selected ten pathways (B).

Note: Another solution is using the filter verb to extract a subset of the result as described in Manipulating enrichment result objects.

39.4 How to extract genes of a specific term/pathway

id <- x$ID[1:3]
id
[1] "hsa04110" "hsa04061" "hsa04974"
x[[id[1]]]
 [1] "8318"  "991"   "9133"  "10403" "890"   "983"   "4085"  "81620" "7272" 
[10] "9212"  "1111"  "9319"  "891"   "4174"  "9232"  "4171"  "993"   "990"  
[19] "5347"  "701"   "9700"  "898"   "23594" "4998"  "9134"  "4175"  "4173" 
[28] "10926" "6502"  "994"   "699"   "4609"  "5111"  "26271" "1869"  "1029" 
[37] "7043" 
geneInCategory(x)[id]
$hsa04110
 [1] "8318"  "991"   "9133"  "10403" "890"   "983"   "4085"  "81620" "7272" 
[10] "9212"  "1111"  "9319"  "891"   "4174"  "9232"  "4171"  "993"   "990"  
[19] "5347"  "701"   "9700"  "898"   "23594" "4998"  "9134"  "4175"  "4173" 
[28] "10926" "6502"  "994"   "699"   "4609"  "5111"  "26271" "1869"  "1029" 
[37] "7043" 

$hsa04061
 [1] "3627"  "10563" "6373"  "4283"  "6362"  "6355"  "2921"  "6364"  "3576" 
[10] "6352"  "3559"  "1230"  "6347"  "3561"  "6351"  "1237"  "1236"  "11009"
[19] "57007" "6387"  "53832" "3572"  "9547"  "1524" 

$hsa04974
 [1] "59272" "5645"  "1299"  "1281"  "50509" "1290"  "477"   "1294"  "1360" 
[10] "1289"  "1292"  "1296"  "23428" "1359"  "1300"  "1287"  "6505"  "2006" 
[19] "7373"  "1307"  "1308"

39.5 Wrap long axis labels

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 39.2 for examples of wrapping long axis labels using numeric width and a custom labeller.

library(ReactomePA)
y <- retry(enrichPathway(de))

p1 <- dotplot(y, label_format = 20) 
p2 <- dotplot(y, label_format = function(x) stringr::str_wrap(x, width=20))
cowplot::plot_grid(p1, p2, ncol=2, labels=c("A", "B"))
Figure 39.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().

39.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: https://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.

39.7 External services, downloads, and version boundaries

Network-backed examples are optional evidence or knowledge acquisition steps; they are not substitutes for checking the local input and statistical design. A failure to reach a service should not be reported as a failed enrichment test.

Service or resource Typical use Record Local or offline fallback
KEGG REST / enrichKEGG() Pathway and module mappings, ID conversion Organism code, endpoint, access date, response status, package version, and mapping coverage Build and reuse a local KEGG database with createKEGGdb, or use a saved GSON resource where appropriate.
Pathview / KEGG diagrams Render a pathway diagram with gene-level values KEGG release or endpoint, Pathview version, license constraints, and input IDs Keep the enrichment table and exported input values; render again when the service is available.
Reactome, WikiPathways, g:Profiler, enrichR, WebGestalt External knowledge or result import Source URL, query parameters, database release, access date, and returned schema Reuse a saved result table or convert a complete local table with the import constructors.
biomaRt, AnnotationHub, and GOA downloads Identifier or annotation acquisition Host/release, query, access date, species, and package version Cache the downloaded mapping or freeze it as TERM2GENE/TERM2NAME or GSON.

For any remote call, distinguish these cases:

  1. Connection failure or timeout: the machine did not obtain the resource; inspect proxy, firewall, DNS, or service availability.
  2. HTTP or API schema failure: the endpoint or response format may have changed; retain the status and response, then check the package release or use a local resource.
  3. Empty or low-coverage result: the request succeeded but the organism, identifier namespace, release, or biological scope may not match.

The reproducible record should include the exact downloaded or serialized resource when redistribution permits it. If the resource cannot be redistributed, retain its accession, release, access date, query, checksum where possible, and the mapping counts used in the analysis. The reproducible report provides the delivery checklist; the KEGG chapter and external-result import chapter show the corresponding APIs.

39.8 Next steps