# How to Do Enrichment Analysis for Non-Model Organisms in R By the end of this guide you will be able to take a list of differentially expressed genes (or proteins, or metabolite-linked genes) from any species that lacks an official annotation package, build your own functional term-to-gene map, and run over-representation analysis (ORA) and gene set enrichment analysis (GSEA) in R using clusterProfiler. You will also be able to sanity-check the results, visualize them as plots and term-similarity networks, and export them to Cytoscape. R for bioinformatics work at this level requires three things installed before you start: 1. A recent release of R (4.3 or later works well with current Bioconductor). 2. The annotation outputs for your species. This guide assumes you have (or can produce) a BLAST or orthology result file and, ideally, an eggNOG-mapper or InterProScan output. 3. Your differential expression table with gene identifiers, log2 fold changes, and adjusted p values. Everything below runs on a normal laptop. The workflow is deliberately general: the same code works for a bacterium, a mollusk, or an insect as long as you supply the right term-to-gene table. ## Why Non-Model Organisms Break Standard Enrichment Workflows Tools like clusterProfiler, g:Profiler, DAVID, and GOrilla ship curated annotations for well-studied organisms. For everything else, curated functional annotations are scarce, and researchers commonly fall back on homology-based transfer, using BLAST to borrow annotations from a related species. That strategy works often enough to be standard practice, but it introduces annotation errors and overlooks functions that are unique to your species [1]. The scale of the gap varies enormously by taxon. In an early transcriptome project for the Manila clam, Ruditapes philippinarum, only 9,747 of 32,606 unique transcripts, roughly 30 percent, could be annotated by sequence similarity, yet that limited annotation was still enough to run meaningful enrichment analysis on differentially expressed genes [2]. Recent resources have improved the picture for some groups. For penaeid shrimp, the ShEnrich platform combined BLASTP against NCBI nr, InterProScan domain prediction, and eggNOG-mapper orthology to annotate 61,287 proteins representing 34,169 genes across five species, yielding 178,595 gene-pathway associations over 442 KEGG pathways and 138,077 proteins mapped to 12,897 GO terms [3]. Limited databases carry a real scientific risk beyond missing terms. In a survey of [protein kinase](/knowledge/molecular-biology/protein-kinase) families across ten invertebrate phyla, restricted database resources led to misidentification and mischaracterization of proteins, which in turn could mislead interpretation of experimental data [4]. Your enrichment results are only as good as the annotation layer underneath them, so building that layer carefully is most of the work. ## Four Annotation Strategies for Non-Model Species ![Flowchart from predicted proteins through BLAST, eggNOG-mapper, InterProScan and a custom package, ending in a TERM2GENE table for clusterProfiler.](/article-images/illustrations/enrichment-analysis-non-model-organism-r-annotation-branching.webp) ### Strategy 1: Homology transfer with BLAST Run blastx or blastp against a curated database such as UniProtKB/Swiss-Prot or NCBI nr, keep hits above an e-value and identity threshold you justify in your methods, and assign the top hit's GO terms or KEGG orthology to your gene. This is the fastest route and the one most people start with. Its weaknesses are propagated errors from automated annotation chains and complete blindness to species-specific genes [1]. It is also sensitive to what is in the database: for poorly sampled taxa, the closest relative in nr may be quite distant. ### Strategy 2: Orthology-based mapping with eggNOG-mapper eggNOG-mapper assigns your genes to orthologous groups and transfers GO terms, KEGG Orthology (KO) identifiers, and COG functional categories in one pass. Because it uses precomputed orthology across thousands of genomes rather than a single best BLAST hit, it tends to be more conservative and more reproducible. ShEnrich used exactly this combination of BLASTP, InterProScan, and eggNOG-mapper as its annotation backbone [3]. For a do-it-yourself workflow, eggNOG-mapper output is the single most useful file format to work from in R. ### Strategy 3: Domain-based annotation with InterProScan InterProScan, maintained at EMBL-EBI, scans your proteins against the InterPro member databases (Pfam, SMART, PROSITE, and others) and returns InterPro and GO terms per domain. Domain-level annotation helps when full-length homology is weak, which is common in divergent lineages. The tradeoff is that a domain is not a gene function, so domain-derived GO terms are broader and noisier. ### Strategy 4: Build a custom annotation package and let clusterProfiler do the rest A recent protocol describes building custom annotation R packages from NCBI genomic data that load as ordinary libraries and plug directly into clusterProfiler for GO and pathway enrichment. The authors demonstrate the approach with an annotation package for Mycobacterium tuberculosis H37Rv, then use it to analyze differentially expressed genes from an RNA-seq dataset of the bacterium's transcriptional response to rifampicin [1]. This is the cleanest long-term solution if you will analyze the same species repeatedly: you invest once in the package, then every downstream analysis becomes as easy as it is for mouse or human. Bioconductor's AnnotationForge package provides the scaffolding (`makeOrgPackage`) for this route. ## Step-by-Step: Enrichment Analysis for a Non-Model Organism in R This walkthrough assumes [RNA-seq differential expression](/knowledge/bioinformatics/rna-seq-differential-expression-deseq2-edger) results from DESeq2, but the same logic applies to proteomics, where published pipelines follow the same sequence of ORA and GSEA steps after differential abundance testing [5]. ### Step 1: Install the R environment ```r install.packages("BiocManager") BiocManager::install(c( "clusterProfiler", "enrichplot", "AnnotationForge", "GO.db", "KEGGREST", "DESeq2", "readr", "dplyr", "tidyr", "stringr", "ggplot2" )) ``` If you are new to the ecosystem, think of this step as the front door of any Bioconductor tutorial: `BiocManager` installs packages built and tested against your R version, which avoids the version mismatches that plague CRAN installs of Bioconductor software. ### Step 2: Generate the annotation files From the shell, with your predicted proteins in `proteins.fasta`: ```bash # Orthology, GO, and KO assignment in one pass (install eggnog-mapper via conda) emapper.py -i proteins.fasta -o my_species --cpu 8 # Optional domain-level support interproscan.sh -i proteins.fasta -f TSV -goterms -pa # Optional raw similarity hits for auditing makeblastdb -in uniprot_sprot.fasta -dbtype prot blastp -query proteins.fasta -db uniprot_sprot.fasta -evalue 1e-5 \ -max_target_seqs 1 -outfmt 6 -out swiss_hits.tsv ``` Keep the raw BLAST hits even if you use eggNOG-mapper as your primary assignment. When an enriched term looks suspicious later, you will want to audit which homolog supported it. ### Step 3: Build the TERM2GENE table in R clusterProfiler's `enricher()` and `GSEA()` functions accept any custom ontology through two data frames. `TERM2GENE` maps each term ID to the genes annotated to it (a many-to-many table with term in the first column, gene in the second). `TERM2NAME` optionally maps term IDs to human-readable names. ```r library(readr) library(dplyr) library(tidyr) library(stringr) raw <- readLines("my_species.emapper.annotations") hdr_i <- grep("^#query", raw)[1] cols <- sub("^#", "", strsplit(raw[hdr_i], "\t")[[1]]) eggnog <- read_tsv(paste(raw[(hdr_i + 1):length(raw)], collapse = "\n"), col_names = cols, na = c("", "NA", "-")) t2g_go <- eggnog %>% select(gene = query, GOs) %>% filter(!is.na(GOs)) %>% separate_rows(GOs, sep = ",") %>% filter(str_detect(GOs, "^GO:\\d{7}$")) %>% select(term = GOs, gene) t2g_ko <- eggnog %>% select(gene = query, KOs) %>% filter(!is.na(KOs)) %>% separate_rows(KOs, sep = ",") %>% mutate(gene = str_replace(gene, "^ko:", "")) %>% select(term = KOs, gene) ``` eggNOG versions differ in which delimiter separates terms within a cell, so open the file once and confirm the separator before trusting `separate_rows()`. A quick coverage check tells you how much of your transcriptome you actually annotated, and this number belongs in your paper: ```r n_annotated <- nrow(distinct(t2g_go, gene)) n_total <- nrow(distinct(eggnog, query)) round(100 * n_annotated / n_total, 1) # e.g. 62.4 means 62.4 percent of genes carry at least one GO term ``` For the Manila clam example cited above, that coverage figure was around 30 percent, and the authors still produced a usable enrichment analysis [2]. Low coverage narrows what you can detect. It does not invalidate the analysis if you handle the background correctly, which is covered in Step 5. ### Step 4: Prepare the gene list and the universe Two inputs matter more than any other parameter in this workflow. The first is the query gene set, for example genes with adjusted p value below 0.05 and absolute log2 fold change above 1. The second is the universe, the set of genes that could possibly have been detected in your experiment. For RNA-seq, that is every gene that passed your expression filter in DESeq2, not every gene in the annotation. ```r de <- read_tsv("de_results.tsv") universe <- de %>% filter(!is.na(padj), baseMean > 10) %>% # match your DESeq2 filtering exactly pull(gene_id) %>% unique() sig <- de %>% filter(padj < 0.05, abs(log2FoldChange) > 1) %>% pull(gene_id) %>% unique() ``` ### Step 5: Run over-representation analysis ```r library(clusterProfiler) library(GO.db) go_names <- data.frame(term = keys(GO.db, "GOID")) go_names$name <- mapIds(GO.db, keys = go_names$term, column = "TERM", keytype = "GOID") ego <- enricher( gene = sig, universe = universe, TERM2GENE = t2g_go, TERM2NAME = go_names, pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.2, minGSSize = 10, maxGSSize = 500 ) head(as.data.frame(ego)) ``` ORA asks a simple question: for each term, is the proportion of your query genes annotated to it higher than expected given the universe, tested with a hypergeometric or Fisher's exact model? The Benjamini-Hochberg method (`pAdjustMethod = "BH"`) controls the false discovery rate across all terms tested. Expected output looks like this (values are illustrative): | ID | Description | GeneRatio | BgRatio | pvalue | p.adjust | geneID | |--|--|--|--|--|--|--| | GO:0006955 | immune response | 18/240 | 95/9800 | 3.1e-06 | 0.0041 | GENE01, GENE07, ... | | GO:0006096 | glycolytic process | 12/240 | 61/9800 | 8.7e-05 | 0.0520 | GENE03, GENE11, ... | `GeneRatio` is query genes in the term over query genes with any annotation. `BgRatio` is all genes in the term over the universe. If `GeneRatio` looks high but `p.adjust` is not significant, the term is simply small and you lack power, which is normal in sparsely annotated species. ### Step 6: Run GSEA on the ranked list ORA discards effect size. GSEA uses the full ranked list, all genes ordered by a statistic such as log2 fold change, and asks whether genes of a term cluster at one end of the ranking. This is more sensitive when individual genes change modestly but a pathway shifts coherently, and it is the standard second analysis in published proteomics and transcriptomics workflows [5]. ```r ranked <- de %>% filter(!is.na(log2FoldChange), baseMean > 10) %>% arrange(desc(log2FoldChange)) %>% distinct(gene_id, .keep_all = TRUE) gene_vector <- ranked$log2FoldChange names(gene_vector) <- ranked$gene_id egsea <- GSEA( gene_vector, TERM2GENE = t2g_go, TERM2NAME = go_names, pAdjustMethod = "BH", pvalueCutoff = 0.05, minGSSize = 10, maxGSSize = 500, seed = TRUE ) ``` Three things to know about this call. The ranked vector must have unique names sorted in decreasing order, so `distinct()` and `arrange()` are not optional. `seed = TRUE` makes the permutation-based p values reproducible across runs. And unlike ORA, GSEA reports a normalized enrichment score whose sign tells you the direction: positive means the term's genes sit at the top of the list (upregulated), negative at the bottom (downregulated). ### Step 7: KEGG pathway enrichment via KO identifiers Most non-model species have no species code in KEGG, but KO identifiers are organism-agnostic. Assign KOs with eggNOG-mapper, map KOs to KEGG pathways with KEGGREST, and run enrichment at the KO level. ```r library(KEGGREST) link <- keggLink("pathway", "ko") ko2path <- data.frame(term = unname(link), gene = sub("^ko:", "", names(link))) path_names <- keggList("pathway", "ko") t2n_kegg <- data.frame(term = names(path_names), name = unname(path_names)) de_ko <- eggnog %>% select(gene = query, KOs) %>% filter(!is.na(KOs)) %>% separate_rows(KOs, sep = ",") %>% mutate(KOs = str_replace(KOs, "^ko:", "")) %>% filter(gene %in% sig) %>% pull(KOs) %>% unique() kegg_res <- enricher( gene = de_ko, TERM2GENE = ko2path, TERM2NAME = t2n_kegg, pAdjustMethod = "BH", pvalueCutoff = 0.05, minGSSize = 5 ) ``` Enriching KOs rather than gene IDs has a bonus: it partly sidesteps gene duplication and annotation-transfer noise, because ortholog groups are more stable than individual transferred labels. ### Step 8: Sanity checks before you believe anything Run these four checks every time: 1. **Coverage.** What fraction of your universe carries any annotation? Below roughly 30 percent, interpret enrichments as statements about the annotated subset only. 2. **Direction agreement.** For terms significant in both ORA and GSEA, do the member genes actually move in the reported direction? Pull them with `ego@geneInTerm` or via `clusterProfiler::gsingerene` equivalents and cross-check against the DE table. 3. **External validation.** A good example comes from deltamethrin-resistant Aedes aegypti transcriptomes, where GO and KEGG enrichment pointed at cytochrome P450 terms and xenobiotic metabolism pathways, and those findings were validated by qPCR with a correlation of R = 0.89 (p = 0.00049) between RNA-seq and qPCR measurements [6]. Your strongest enriched terms should suggest a similar testable prediction. 4. **Reproducibility.** Set the seed for GSEA, save the exact TERM2GENE table with a date stamp, and record package versions with `sessionInfo()`. Annotations change, and an analysis that cannot be regenerated later is an analysis you cannot defend at review. ## Parameter Reference | Parameter | Applies to | What it does | When to change it | |--|--|--|--| | `universe` | ORA | Defines the background set for the Fisher test | Always set it. Never let clusterProfiler default to all genes in TERM2GENE | | `pAdjustMethod` | Both | Multiple-testing correction, BH is the field standard | Keep BH unless you have a specific reason otherwise | | `pvalueCutoff` | Both | Raw p value threshold for reporting | Loosen to 0.1 when exploring, but filter on `p.adjust` in the paper | | `qvalueCutoff` | Both | FDR-based q value threshold on top of BH | Set to 1 if you want to filter on `p.adjust` alone | | `minGSSize` | Both | Minimum genes per term to test | Lower to 5 for sparse annotations, but expect instability | | `maxGSSize` | Both | Maximum genes per term to test | Lower to 500 to suppress terms so broad they are uninterpretable | | `seed` | GSEA | Fixes permutation randomness | Always TRUE for reproducible results | | `by` | GSEA | Algorithm backend, fgsea is default and fast | Leave default unless reproducing legacy results | ## Visualization: From Dot Plots to Cytoscape Networks ```r library(enrichplot) dotplot(ego, showCategory = 20, split = NULL) + ggplot2::ggtitle("ORA, control vs treated") ridgeplot(egsea, showCategory = 20) # fold-change distribution per term edo_pair <- pairwise_termsim(egsea) emapplot(edo_pair, showCategory = 30) # terms clustered by shared genes cnetplot(ego, showCategory = 5) # gene-concept network for top terms ``` Bar plots and dot plots summarize, but they hide relationships between terms. Enrichment maps connect terms that share genes, which turns a list of 200 significant terms into a handful of interpretable functional modules. A Current Protocols workflow extends this in three useful directions: calculating term similarity explicitly, exporting networks to Cytoscape with the `cyemapplot` function, and overlaying node-level properties such as pathway coverage (what fraction of a pathway's genes are in your query) and expression directionality, so you can distinguish complete pathway perturbation from partial drift [7]. If your term list is large, do the exploratory clustering in R and the final figure in Cytoscape, where manual curation is practical. ## Adapting Other Tools When clusterProfiler Is Not Enough DAVID's functional classification tool can work with custom annotations, but it requires an R script to import custom gene identifiers and novel annotation information, and the clustering behavior depends heavily on the kappa score. In rainbow trout datasets, the number of genes harbored in enriched modules generally decreased as kappa increased, with the peak number of genes in enriched modules occurring at a kappa of 0.5 for comparisons with 3,340 and 1,313 differentially expressed genes and at 0.3 for a comparison with 555 genes. Large non-enriched modules could be reanalyzed at a higher kappa to resolve smaller enriched modules [8]. The practical takeaway is that kappa should be tuned to the size of your gene list, roughly 0.3 for short lists and 0.5 for large ones, and reported in your methods. If your audience is wet-lab collaborators, web-based options lower the barrier. RNfuzzyApp provides an R Shiny interface for differential expression, time-series clustering with Mfuzz, and cluster enrichment, and it was built explicitly because most meta-analysis and enrichment tooling is tailored to model organisms and not easily adaptable to other species [9]. Similar species-specific platforms now exist for individual taxonomic groups, such as the shrimp-focused ShEnrich, which builds GMT libraries per species and runs enrichment through clusterProfiler behind a database front end [3]. Checking whether a platform exists for your taxon can save you the entire annotation build. The same custom-term-to-gene logic extends past transcriptomics. Quantitative proteomics pipelines perform ORA and GSEA on differentially abundant proteins with fully reproducible R code [5], and metabolomics workflows in R combine enrichment against prior knowledge networks with multi-omics integration [10]. In every case the statistical machinery is identical. Only the term-to-gene (or term-to-metabolite) mapping changes. ## Common Mistakes and Limitations **Leaving the universe unset.** The single most common error. Without `universe`, the background becomes every gene in your annotation table, which biases results toward well-annotated, highly expressed genes. The universe must be the set of genes that had a chance to be called significant in your experiment. **Identifier drift.** Gene IDs change between annotation runs, Ensembl releases, and database versions. If half your DE gene IDs silently fail to match TERM2GENE, clusterProfiler warns with `No gene can be mapped` only when the failure is total. Partial mismatches just quietly shrink your query. Before running enrichment, assert that a known fraction of your query matches: `mean(sig %in% t2g_go$gene)` should be close to your overall annotation coverage. **Treating transferred annotations as ground truth.** BLAST top-hit transfer mixes annotation errors into your term map, and distant homologs can carry the wrong function. Audit the genes behind any headline term before building a narrative on it. **Forgetting that your species biology may be genuinely absent.** Species-specific genes, recent duplications, and taxon-restricted pathways have no ortholog annotation at all. Enrichment will be silent about them. Complement enrichment with direct inspection of unannotated differentially expressed genes, for example by domain search against nr. **Running ORA on a ranked statistic you should have used GSEA for, and vice versa.** If most genes shift modestly in the same direction, ORA on a log2 fold change threshold will return nothing while GSEA detects coordinated movement. If you have a clean thresholded hit list from a targeted experiment, ORA is the right and simpler test. **Over-testing tiny terms.** With `minGSSize` below about 5 in sparse annotations, a single gene can make a term significant. Prefer 5 to 10 as a floor, and always report term sizes alongside p values. **Ignoring direction in ORA.** ORA pools up- and downregulated genes. A term can be enriched because its genes scatter in both directions. Either split the query by direction, or rely on GSEA for directional claims. **Discarding duplicates too late or too early.** GSEA aborts with `geneID should be unique` when the ranked vector has duplicated names. Collapse duplicate gene IDs by keeping the row with the largest absolute statistic, and document the collapse rule. **Untested kappa in DAVID-style clustering.** Picking kappa by habit rather than by sensitivity analysis can merge everything into one large uninformative module or shatter it into fragments. The trout analysis showed the optimum shifts with gene list size, so scan kappa and report the sensitivity [8]. Finally, remember that enrichment analysis describes your annotation layer as much as your biology. Every claim should be phrased as genes annotated to term X are overrepresented, not as pathway X is activated. Individual experimental follow-up, for example qPCR of candidate genes as in the Aedes aegypti study [6], remains the arbiter. ## Frequently Asked Questions ### Can I run clusterProfiler on a species with no official annotation package? Yes. Use `enricher()` and `GSEA()` with a custom `TERM2GENE` data frame built from BLAST, eggNOG-mapper, or InterProScan results, and optionally build a permanent OrgDb-style package with AnnotationForge [1]. ### What should I use as the background universe? Use every gene that passed your expression filters and had a chance to enter the differential expression test, matching the exact filters you applied in DESeq2 or edgeR. Do not use all genes in the annotation or all genes in the genome. ### Is BLAST-based annotation transfer reliable? It is the standard first pass, but it propagates errors and misses species-specific functions [1]. Use conservative thresholds, audit top hits for important terms, and prefer orthology-based assignment when you can get it. ### Should I use ORA or GSEA? Use ORA when you have a defined hit list from thresholded results and GSEA when you want to use the full ranked spectrum of effect sizes. Running both and reporting where they agree is the strongest presentation. ### How do I choose a kappa score for DAVID functional clustering in a non-model organism? Scan kappa and report sensitivity rather than fixing it by default. Published trout analyses found the peak number of genes in enriched modules at kappa 0.5 for large DEG lists and 0.3 for a small list of 555 genes [8]. ### Why do I get almost no enriched terms? Check annotation coverage first, then the universe, then identifier matching. Sparse annotations and partial ID mismatches shrink the effective query, and overly strict `maxGSSize` or direction-splitting can remove the signal you expect. ### Can I do KEGG pathway enrichment for a species not in KEGG? Yes, through KO identifiers. Assign KOs with eggNOG-mapper, map KOs to pathways with KEGGREST, and run `enricher()` at the KO level with a custom `TERM2GENE` table. ### How do I visualize enrichment results for a paper-quality figure? Start with dot plots and ridge plots in R, cluster terms by gene overlap with `emapplot` and `pairwise_termsim`, then export the network to Cytoscape for interactive curation, pathway coverage overlays, and expression directionality [7]. ## Sources 1. [Enhanced RNA-Seq Expression Profiling and Functional Enrichment in Non-model Organisms Using Custom Annotations.](https://pubmed.ncbi.nlm.nih.gov/42367215/) 2. [Transcriptome sequencing and microarray development for the Manila clam, Ruditapes philippinarum: genomic tools for environmental monitoring.](https://pubmed.ncbi.nlm.nih.gov/21569398/) 3. [ShEnrich: a functional annotation database and enrichment platform for penaeid shrimp.](https://pubmed.ncbi.nlm.nih.gov/42963765/) 4. [Crustacean Protein Kinases A and C: Bioinformatic Characterization in Decapods and Other Non-Model Organisms.](https://pubmed.ncbi.nlm.nih.gov/41226621/) 5. [Computational and Statistical Framework for Quantitative Proteomics Analysis.](https://pubmed.ncbi.nlm.nih.gov/42681023/) 6. [Transcriptome Analysis of Deltamethrin Resistance Mechanisms in Aedes aegypti Populations from Metro Manila, Philippines.](https://pubmed.ncbi.nlm.nih.gov/42792987/) 7. [Enhancing Interpretation of Functional Enrichment Analysis Results with Advanced Network Visualizations.](https://pubmed.ncbi.nlm.nih.gov/42905494/) 8. [Use of DAVID algorithms for clustering custom annotated gene lists in a non-model organism, rainbow trout.](https://pubmed.ncbi.nlm.nih.gov/29361978/) 9. [RNfuzzyApp: an R shiny RNA-seq data analysis app for visualisation, differential expression analysis, time-series clustering and enrichment analysis.](https://pubmed.ncbi.nlm.nih.gov/35186266/) 10. [Multi-Omics Guided Pathway and Network Analysis of Clinical Metabolomics and Proteomics Data.](https://pubmed.ncbi.nlm.nih.gov/42680979/)