GO and KEGG enrichment for a non-model organism

Project Overview

The client required an RNA-seq analysis for the haploid non-model organism Marchantia polymorpha. Because standard functional annotation resources were limited, the reference annotations had to be carefully curated and converted for downstream enrichment analysis. This project demonstrates my workflow for constructing reproducible GO annotations and KO-based pathway mappings, where KO identifiers provide the standardized link between non-model organism genes and KEGG pathways for over-representation analysis and GSEA, using the relationship “gene → KO ID → KEGG pathway”.

Key Takeaways

• In a non-model organism such as Marchantia polymorpha, enrichment analysis is not just a downstream step; it depends heavily on how the reference annotation is constructed.

• Because M. polymorpha is haploid and has sex-specific genomic resources, the RNA-seq pipeline requires extra care when choosing the genome, GTF annotation, and functional annotation reference.

• The gene IDs used in the genome/GTF annotation must be matched to the gene IDs in the GO and KO annotation tables. If the IDs are inconsistent, downstream enrichment analysis can fail silently or produce incomplete biological interpretation.

• For KO-based pathway analysis, the key relationship is:

  gene → KO ID → KEGG pathway

Reference genome and annotations

The genome was obtained from NCBI. Gene annotations were obtained from MarpolBase, using the MpTak2_v7.1 female Tak-2 reference genome.

• Index of /data/v7.1_reference_genomes/
• NCBI genome reference

Goal

• Process the annotation of MarpolBase to run enrichment analysis using R package, clusterprofiler, with “TERM2NAME” option.
• female Tak-2 reference

Pipeline - GO

1. Turn MpTak2_v7.1.GO.tsv from transcript-level to gene-level
2. Tidy into gene-level TERM2GENE format
3. Tidy TERM2NAME: mapping GO ID to GO term name
4. Save as RData with the following COLUMNs:
   • TERM2GENE
   • TERM2NAME
   • UNIVERSE
5. Results:
   • MpTak2_v7.1.GO.tsv original rows: 51353
   • gene-level rows: 10407
   • transcript-gene mapping rows: 13306
   • unmapped transcripts: 0
   • TERM2GENE rows: 39892
   • TERM2NAME rows: 4267
   • UNIVERSE size: 10407

Pipeline - KO

1. MarpolBase annotation is described as below:

   # This file contains annotation only from KEGG,KOG,Pfam,NCBIfam.
   # Annotations are joined with '; '.
   Mp1g00070.1	KEGG:K10683:BRCA1-associated RING domain protein 1; KOG:KOG0311:Predicted E3 ubiquitin ligase, C-term missing, [O]; KOG:KOG4362:Transcriptional regulator BRCA1, [LK]; Pfam:PF00533:BRCA1 C Terminus (BRCT) domain; Pfam:PF13771:PHD-like zinc-binding domain; Pfam:PF13923:Zinc finger, C3HC4 type (RING finger)
   Mp1g00080.1	KEGG:K01599:uroporphyrinogen decarboxylase [EC:4.1.1.37]; KOG:KOG2872:Uroporphyrinogen decarboxylase, [H]; NCBIfam:TIGR01464:uroporphyrinogen decarboxylase; Pfam:PF01208:Uroporphyrinogen decarboxylase (URO-D)
   Mp1g00090.1	KOG:KOG2896:UV radiation resistance associated protein, [R]; Pfam:PF10186:Vacuolar sorting 38 and autophagy-related subunit 14
   ...

   1. COLUMN 1 is transcript ID → Mp1g00070.1, Mp1g00080.1, Mp1g00090.1…etc → It has to match ID in genome reference
   2. KEGG:Kxxxx is the KO we want.

2. We need transcriptID ↔ geneID mapping information. See GTF:

   #gtf-version 2.2
   #!genome-build MpTak2_v7.1
   #!genome-build-accession NCBI_Assembly:GCA_037833965.1
   AP031359.1      DDBJ    gene    143608  197944  .       -       .       gene_id "MPTK2_Ug00010"; transcript_id ""; gbkey "Gene"; gene_biotype "protein_coding"; locus_tag "MPTK2_Ug00010"; 
   AP031359.1      DDBJ    transcript      143608  197944  .       -       .       gene_id "MPTK2_Ug00010"; transcript_id "unassigned_transcript_1"; gbkey "mRNA"; locus_tag "MPTK2_Ug00010"; note "submitter_gene_id: MpUg00010, submitter_transcript_id: MpUg00010.1"; transcript_biotype "mRNA"; 
   AP031359.1      DDBJ    exon    197518  197944  .       -       .       gene_id "MPTK2_Ug00010"; transcript_id "unassigned_transcript_1"; locus_tag "MPTK2_Ug00010"; note "submitter_gene_id: MpUg00010, submitter_transcript_id: MpUg00010.1"; transcript_biotype "mRNA"; exon_number "1"; 
   ...

   1. GTF submitter_gene_id in “note” is what I need

      note "submitter_gene_id: MpUg00010, submitter_transcript_id: MpUg00010.1"

3. TIDY annotation reference（MpTak2_v7.1.func_annotation.1_line.tsv） turn transcript-level into gene-level, we produced following files:

• Gene-level annotation

  Mp1g00070	KEGG:K10683:BRCA1-associated RING domain protein 1; KOG:KOG0311:Predicted E3 ubiquitin ligase, C-term missing, [O]; KOG:KOG4362:Transcriptional regulator BRCA1, [LK]; Pfam:PF00533:BRCA1 C Terminus (BRCT) domain; Pfam:PF13771:PHD-like zinc-binding domain; Pfam:PF13923:Zinc finger, C3HC4 type (RING finger)
  Mp1g00080	KEGG:K01599:uroporphyrinogen decarboxylase [EC:4.1.1.37]; KOG:KOG2872:Uroporphyrinogen decarboxylase, [H]; NCBIfam:TIGR01464:uroporphyrinogen decarboxylase; Pfam:PF01208:Uroporphyrinogen decarboxylase (URO-D)
  Mp1g00090	KOG:KOG2896:UV radiation resistance associated protein, N-term missing, [R]; Pfam:PF10186:Vacuolar sorting 38 and autophagy-related subunit 14 || KOG:KOG2896:UV radiation resistance associated protein, [R]; Pfam:PF10186:Vacuolar sorting 38 and autophagy-related subunit 14

• Mapping information

  transcript_id	gene_id
  Mp1g00070.1	Mp1g00070
  Mp1g00080.1	Mp1g00080
  Mp1g00090.1	Mp1g00090
  Mp1g00090.2	Mp1g00090
  ...

• Unmapped information

1. Convert gene-level annotation into enricher() / GSEA() RData → reproducible

   • download KEGG annotations

   curl -L https://rest.kegg.jp/link/pathway/ko -o ko2path.tsv
   curl -L https://rest.kegg.jp/list/pathway -o pathway_names.tsv

   • data process

Code (build_custom_kegg_cache.R)

#!/usr/bin/env Rscript

suppressPackageStartupMessages({
  library(data.table)
})

options(timeout = 100000)

script_dir <- function() {
  cmd_args <- commandArgs(trailingOnly = FALSE)
  file_arg <- grep("^--file=", cmd_args, value = TRUE)
  if (length(file_arg) > 0L) {
    return(dirname(normalizePath(sub("^--file=", "", file_arg[[1]]), winslash = "/", mustWork = TRUE)))
  }
  getwd()
}

get_arg <- function(flag, default = NULL) {
  args <- commandArgs(trailingOnly = TRUE)

  idx <- match(flag, args)
  if (!is.na(idx)) {
    if (idx == length(args)) {
      stop("Missing value for ", flag, call. = FALSE)
    }
    return(args[[idx + 1]])
  }

  prefix <- paste0(flag, "=")
  hit <- args[startsWith(args, prefix)]
  if (length(hit) > 0L) {
    return(sub(prefix, "", hit[[1]], fixed = TRUE))
  }

  default
}

has_flag <- function(flag) {
  args <- commandArgs(trailingOnly = TRUE)
  flag %in% args || any(startsWith(args, paste0(flag, "=")))
}

normalize_missing <- function(x) {
  x <- trimws(as.character(x))
  x[x %in% c("", "\"\"", "NA")] <- NA_character_
  x
}

split_ko_terms <- function(x) {
  x <- normalize_missing(x)
  if (is.na(x)) {
    return(character(0))
  }

  tokens <- unlist(strsplit(x, "\\s*;\\s*", perl = TRUE), use.names = FALSE)
  tokens <- trimws(tokens)
  tokens <- tokens[nzchar(tokens)]
  tokens <- sub("^KEGG:", "", tokens)
  tokens <- sub("^ko:", "", tokens, ignore.case = TRUE)
  unique(tokens[nzchar(tokens)])
}

extract_ko_terms_from_annotation <- function(x, ko_regex) {
  x <- normalize_missing(x)
  if (is.na(x)) {
    return(character(0))
  }

  hits <- regmatches(x, gregexpr(ko_regex, x, perl = TRUE))[[1]]
  if (length(hits) == 0L) {
    return(character(0))
  }

  hits <- sub("^KEGG:", "", hits)
  hits <- sub("^ko:", "", hits, ignore.case = TRUE)
  unique(hits[nzchar(hits)])
}

collect_row_kos <- function(ko_text, annotation_text, ko_regex) {
  unique(c(
    split_ko_terms(ko_text),
    extract_ko_terms_from_annotation(annotation_text, ko_regex)
  ))
}

read_gene_level_table <- function(path) {
  if (!file.exists(path)) {
    stop("Gene-level annotation file not found: ", path, call. = FALSE)
  }

  dt <- data.table::fread(
    path,
    header = TRUE,
    sep = "\t",
    quote = "",
    fill = TRUE,
    data.table = FALSE
  )

  if (!("gene_id" %in% names(dt))) {
    dt <- data.table::fread(
      path,
      header = FALSE,
      sep = "\t",
      quote = "",
      fill = TRUE,
      data.table = FALSE
    )
    if (ncol(dt) < 2L) {
      stop("Gene-level annotation file must contain at least two columns.", call. = FALSE)
    }
    dt <- as.data.frame(dt[, 1:2], stringsAsFactors = FALSE)
    names(dt) <- c("gene_id", "annotation")
  } else {
    if (!("annotation" %in% names(dt))) {
      stop("Gene-level annotation file must contain an `annotation` column.", call. = FALSE)
    }
    dt <- as.data.frame(dt[, c("gene_id", "annotation")], stringsAsFactors = FALSE)
  }

  if (nrow(dt) == 0L) {
    stop("Gene-level annotation file is empty: ", path, call. = FALSE)
  }

  dt
}

read_kegg_table <- function(source, label) {
  if (is.null(source) || !nzchar(source)) {
    stop("Missing source for ", label, call. = FALSE)
  }

  tbl <- if (file.exists(source)) {
    data.table::fread(
      source,
      header = FALSE,
      sep = "\t",
      quote = "",
      fill = TRUE,
      data.table = FALSE
    )
  } else {
    tryCatch(
      data.table::fread(
        source,
        header = FALSE,
        sep = "\t",
        quote = "",
        fill = TRUE,
        data.table = FALSE
      ),
      error = function(e) {
        stop(
          "Failed to read ", label, " from source `", source, "`: ",
          conditionMessage(e),
          call. = FALSE
        )
      }
    )
  }

  if (nrow(tbl) == 0L) {
    stop("No rows returned for ", label, " from source `", source, "`.", call. = FALSE)
  }

  tbl
}

normalize_ko2path_table <- function(tbl, source_label) {
  if (ncol(tbl) < 2L) {
    stop("Unexpected KO-to-pathway table shape from source: ", source_label, call. = FALSE)
  }

  col1 <- as.character(tbl[[1]])
  col2 <- as.character(tbl[[2]])

  col1_ko <- sum(grepl("^ko:", col1))
  col1_path <- sum(grepl("^path:", col1))
  col2_ko <- sum(grepl("^ko:", col2))
  col2_path <- sum(grepl("^path:", col2))

  if (col1_ko > 0L && col2_path > 0L && col1_ko >= col1_path && col2_path >= col2_ko) {
    ko_raw <- col1
    pathway_raw <- col2
  } else if (col2_ko > 0L && col1_path > 0L && col2_ko >= col2_path && col1_path >= col1_ko) {
    ko_raw <- col2
    pathway_raw <- col1
  } else if (col1_ko + col2_path >= col2_ko + col1_path) {
    ko_raw <- col1
    pathway_raw <- col2
  } else {
    ko_raw <- col2
    pathway_raw <- col1
  }

  data.frame(
    pathway = sub("^path:", "", pathway_raw),
    ko = sub("^ko:", "", ko_raw),
    stringsAsFactors = FALSE
  )
}

build_custom_kegg_cache <- function(
  gene_level_path,
  cache_path,
  species,
  keytype = "custom-ko",
  gene_column = "gene_id",
  ko_column = "ko_terms",
  annotation_column = "annotation",
  ko_regex = "KEGG:K\\d{5}",
  readable_keytype = "SYMBOL",
  pathview_species = "ko",
  pathview_node_basis = "KO",
  pathview_node_sum = "max.abs",
  pathview_file_tag = "ko_maxabs",
  ko_link_source = "https://rest.kegg.jp/link/pathway/ko",
  path_name_source = "https://rest.kegg.jp/list/pathway",
  overwrite = FALSE
) {
  if (file.exists(cache_path) && !isTRUE(overwrite)) {
    stop(
      "Cache already exists: ", cache_path, "\n",
      "Use --overwrite to replace it.",
      call. = FALSE
    )
  }

  gene_level <- read_gene_level_table(gene_level_path)

  if (!gene_column %in% names(gene_level)) {
    stop("Missing gene column `", gene_column, "` in: ", gene_level_path, call. = FALSE)
  }
  if (!ko_column %in% names(gene_level) && !annotation_column %in% names(gene_level)) {
    stop(
      "The input file must contain either `", ko_column, "` or `", annotation_column, "` for KO extraction.",
      call. = FALSE
    )
  }

  gene_ids <- as.character(gene_level[[gene_column]])
  ko_source <- if (ko_column %in% names(gene_level)) as.character(gene_level[[ko_column]]) else rep(NA_character_, nrow(gene_level))
  annotation_source <- if (annotation_column %in% names(gene_level)) as.character(gene_level[[annotation_column]]) else rep(NA_character_, nrow(gene_level))

  gene_kos <- mapply(
    collect_row_kos,
    ko_text = ko_source,
    annotation_text = annotation_source,
    MoreArgs = list(ko_regex = ko_regex),
    SIMPLIFY = FALSE,
    USE.NAMES = FALSE
  )

  gene2ko <- data.frame(
    gene_id = rep(gene_ids, lengths(gene_kos)),
    ko = unlist(gene_kos, use.names = FALSE),
    stringsAsFactors = FALSE
  )
  gene2ko <- gene2ko[!is.na(gene2ko$gene_id) & nzchar(gene2ko$gene_id), , drop = FALSE]
  gene2ko <- gene2ko[!is.na(gene2ko$ko) & nzchar(gene2ko$ko), , drop = FALSE]
  gene2ko <- unique(gene2ko)

  if (nrow(gene2ko) == 0L) {
    stop("No KO terms could be extracted from: ", gene_level_path, call. = FALSE)
  }

  ko2path <- normalize_ko2path_table(
    read_kegg_table(ko_link_source, "KO-to-pathway mapping"),
    ko_link_source
  )
  ko2path <- unique(ko2path[!is.na(ko2path$pathway) & nzchar(ko2path$pathway) & !is.na(ko2path$ko) & nzchar(ko2path$ko), , drop = FALSE])

  matched_kos <- intersect(unique(gene2ko$ko), unique(ko2path$ko))
  if (length(matched_kos) == 0L) {
    stop(
      "No overlapping KO terms were found between the gene-level annotation and the KO-to-pathway table.\n",
      "Unique KO terms in annotation: ", length(unique(gene2ko$ko)), "\n",
      "Rows in KO-to-pathway table: ", nrow(ko2path), "\n",
      "If the KEGG REST source is unavailable, supply local files with ",
      "--ko2path-file and --path-name-file.",
      call. = FALSE
    )
  }
  gene_pathway <- merge(gene2ko, ko2path, by = "ko", all = FALSE)
  if (nrow(gene_pathway) == 0L) {
    stop("No gene-to-pathway mappings were produced from the provided KO terms.", call. = FALSE)
  }
  gene_pathway <- unique(data.frame(
    term = gene_pathway$pathway,
    gene = gene_pathway$gene_id,
    stringsAsFactors = FALSE
  ))

  path_names <- read_kegg_table(path_name_source, "pathway names")
  if (ncol(path_names) < 2L) {
    stop("Unexpected pathway name table shape from source: ", path_name_source, call. = FALSE)
  }
  names(path_names)[1:2] <- c("term", "name")
  path_names <- as.data.frame(path_names, stringsAsFactors = FALSE)
  path_names$term <- sub("^path:", "", path_names$term)
  path_names <- unique(path_names[!is.na(path_names$term) & nzchar(path_names$term), , drop = FALSE])

  custom_kegg_resources <- list(
    keytype = keytype,
    species = species,
    annotation_path = gene_level_path,
    gene_level_path = gene_level_path,
    cache_path = cache_path,
    ko_link_source = ko_link_source,
    path_name_source = path_name_source,
    ko_regex = ko_regex,
    TERM2GENE = gene_pathway,
    TERM2NAME = path_names,
    UNIVERSE = unique(gene_pathway$gene),
    gene2ko = gene2ko,
    ko2path = ko2path,
    readable_keytype = readable_keytype,
    pathview_species = pathview_species,
    pathview_node_basis = pathview_node_basis,
    pathview_node_sum = pathview_node_sum,
    pathview_file_tag = pathview_file_tag,
    created_at = Sys.time()
  )

  dir.create(dirname(cache_path), recursive = TRUE, showWarnings = FALSE)
  save(custom_kegg_resources, file = cache_path)

  message("Input rows: ", nrow(gene_level))
  message("Unique genes: ", length(unique(gene_ids)))
  message("Gene->KO pairs: ", nrow(gene2ko))
  message("Gene->pathway pairs: ", nrow(gene_pathway))
  message("Unique KO terms: ", length(unique(gene2ko$ko)))
  message("Unique pathways: ", length(unique(gene_pathway$term)))
  invisible(custom_kegg_resources)
}

main <- function() {
  default_gene_level <- "/mnt/space4/RNAseq/Ref/Marchantia_polymorpha/MpTak_v7.1.func_annotation.gene_level.tsv"
  default_species <- "common_liverwort"

  gene_level_path <- get_arg("--gene-level", default_gene_level)
  species <- get_arg("--species", default_species)
  default_cache_path <- file.path(script_dir(), "KEGG_database", "KO", paste0(species, ".RData"))
  cache_path <- get_arg("--cache-path", default_cache_path)

  if (!identical(species, default_species) && identical(gene_level_path, default_gene_level)) {
    stop(
      "When using --species ", species, ", you must also pass --gene-level explicitly.",
      call. = FALSE
    )
  }

  gene_column <- get_arg("--gene-column", "gene_id")
  ko_column <- get_arg("--ko-column", "ko_terms")
  annotation_column <- get_arg("--annotation-column", "annotation")
  ko_regex <- get_arg("--ko-regex", "KEGG:K\\d{5}")
  readable_keytype <- get_arg("--readable-keytype", "SYMBOL")
  pathview_species <- get_arg("--pathview-species", "ko")
  pathview_node_basis <- get_arg("--pathview-node-basis", "KO")
  pathview_node_sum <- get_arg("--pathview-node-sum", "max.abs")
  pathview_file_tag <- get_arg("--pathview-file-tag", "ko_maxabs")
  ko_link_source <- get_arg("--ko-link-source", "https://rest.kegg.jp/link/pathway/ko")
  path_name_source <- get_arg("--path-name-source", "https://rest.kegg.jp/list/pathway")
  ko2path_file <- get_arg("--ko2path-file", NULL)
  if (!is.null(ko2path_file) && nzchar(ko2path_file)) {
    ko_link_source <- ko2path_file
  }
  path_name_file <- get_arg("--path-name-file", NULL)
  if (!is.null(path_name_file) && nzchar(path_name_file)) {
    path_name_source <- path_name_file
  }
  keytype <- get_arg("--keytype", "custom-ko")

  overwrite <- has_flag("--overwrite")

  build_custom_kegg_cache(
    gene_level_path = gene_level_path,
    cache_path = cache_path,
    species = species,
    keytype = keytype,
    gene_column = gene_column,
    ko_column = ko_column,
    annotation_column = annotation_column,
    ko_regex = ko_regex,
    readable_keytype = readable_keytype,
    pathview_species = pathview_species,
    pathview_node_basis = pathview_node_basis,
    pathview_node_sum = pathview_node_sum,
    pathview_file_tag = pathview_file_tag,
    ko_link_source = ko_link_source,
    path_name_source = path_name_source,
    overwrite = overwrite
  )

  message("Gene-level input: ", gene_level_path)
  message("Species: ", species)
  message("Cache path: ", cache_path)
  message("Keytype: ", keytype)
  message("Overwrite: ", overwrite)
  message("KO regex: ", ko_regex)
  message("KO-to-pathway source: ", ko_link_source)
  message("Pathway-name source: ", path_name_source)
  message("Cache built successfully.")
}

main()

Rscript build_custom_kegg_cache.R \
  --gene-level /path/to/your_gene_level.tsv \
  --species your_species \
  --cache-path your_species.RData \
  --ko2path-file /path/to/ko2path.tsv \
  --path-name-file /path/to/path_names.tsv

• “gene ↔ KO” format: gene2ko
• import ko2path.tsv → ”KO ↔ pathway” mapping format
• concat gene2ko and ko2path → “gene ↔ pathway” mapping information → this is what TERM2GENE need for enricher()
• import path_names.tsv → convert to TERM2NAME → enrichment results show the name of pathways

1. Result
   • gene_level input: 12504 rows
   • gene2ko: 6565 rows
   • TERM2GENE: 29544 rows
   • TERM2NAME: 585 rows
   • UNIVERSE: 3866 genes