tests/testthat/test-pathway_gsea.R

In ggpicrust2: Make 'PICRUSt2' Output Analysis and Visualization Easier

# Tests for pathway_gsea function

# Helper: create standard GSEA test data
create_gsea_test_data <- function(n_features = 50, n_samples = 10) {
  set.seed(123)
  abundance <- matrix(rnorm(n_features * n_samples), nrow = n_features, ncol = n_samples)
  colnames(abundance) <- paste0("Sample", 1:n_samples)
  rownames(abundance) <- paste0("K", sprintf("%05d", 1:n_features))

  metadata <- data.frame(
    sample_name = colnames(abundance),
    group = factor(rep(c("Control", "Treatment"), each = n_samples / 2)),
    stringsAsFactors = FALSE
  )
  rownames(metadata) <- metadata$sample_name

  list(abundance = abundance, metadata = metadata)
}

test_that("pathway_gsea works with fgsea method", {
  skip_if_not_installed("fgsea")

  test_data <- create_gsea_test_data()

  local_mocked_bindings(
    prepare_gene_sets = function(...) {
      list("path:ko00001" = c("K00001", "K00002"), "path:ko00002" = c("K00002", "K00003"))
    },
    run_fgsea = function(...) {
      data.frame(
        pathway_id = c("path:ko00001", "path:ko00002"),
        pathway_name = c("path:ko00001", "path:ko00002"),
        size = c(2, 2), ES = c(0.5, -0.3), NES = c(1.2, -0.8),
        pvalue = c(0.01, 0.05), p.adjust = c(0.02, 0.1),
        leading_edge = c("K00001;K00002", "K00003"),
        stringsAsFactors = FALSE
      )
    }
  )

  result <- pathway_gsea(
    abundance = test_data$abundance,
    metadata = test_data$metadata,
    group = "group",
    method = "fgsea"
  )

  expect_s3_class(result, "data.frame")
  expect_true(all(c("pathway_id", "NES", "pvalue", "method") %in% colnames(result)))
})

test_that("pathway_gsea works with camera method", {
  skip_if_not_installed("limma")

  test_data <- create_gsea_test_data()
  abundance <- abs(test_data$abundance) * 100

  local_mocked_bindings(
    prepare_gene_sets = function(...) {
      list(
        "ko00010" = paste0("K", sprintf("%05d", 1:10)),
        "ko00020" = paste0("K", sprintf("%05d", 11:20))
      )
    }
  )

  result <- pathway_gsea(
    abundance = abundance,
    metadata = test_data$metadata,
    group = "group",
    method = "camera",
    pathway_type = "KEGG"
  )

  expect_s3_class(result, "data.frame")
  expect_true(all(c("pathway_id", "direction", "pvalue") %in% colnames(result)))
  expect_equal(unique(result$method), "camera")
})

test_that("pathway_gsea validates inputs correctly", {
  test_data <- create_gsea_test_data()

  expect_error(pathway_gsea(abundance = "invalid", metadata = test_data$metadata, group = "group"),
               "'abundance' must be a data frame or matrix")
  expect_error(pathway_gsea(abundance = test_data$abundance, metadata = "invalid", group = "group"),
               "'metadata' must be a data frame")
  expect_error(pathway_gsea(abundance = test_data$abundance, metadata = test_data$metadata, group = "invalid_group"),
               "not found in metadata")
  expect_error(pathway_gsea(abundance = test_data$abundance, metadata = test_data$metadata, group = "group", method = "invalid"),
               "method must be one of")
})

test_that("prepare_gene_sets works for KEGG and MetaCyc pathway types", {
  gene_sets_kegg <- prepare_gene_sets("KEGG")
  expect_type(gene_sets_kegg, "list")
  expect_true(length(gene_sets_kegg) > 0)
  expect_false("ko01001" %in% names(gene_sets_kegg))
  expect_false("ko99980" %in% names(gene_sets_kegg))

  gene_sets_metacyc <- prepare_gene_sets("MetaCyc")
  expect_type(gene_sets_metacyc, "list")
})

test_that("calculate_rank_metric aligns by sample column and handles constant t-test rows", {
  calc <- getFromNamespace("calculate_rank_metric", "ggpicrust2")
  abundance <- matrix(
    c(
      1, 1, 1, 1,
      1, 2, 8, 9
    ),
    nrow = 2,
    byrow = TRUE,
    dimnames = list(c("K00001", "K00002"), paste0("S", 1:4))
  )
  metadata <- data.frame(
    sample_name = paste0("S", 1:4),
    group = c("A", "A", "B", "B"),
    stringsAsFactors = FALSE
  )

  metric <- calc(abundance, metadata, "group", method = "t_test")
  expect_named(metric, rownames(abundance))
  expect_equal(unname(metric["K00001"]), 0)
  expect_true(is.finite(metric["K00002"]))
})

test_that("pathway_gsea validates PICRUSt2 #NAME input after stripping the ID column", {
  skip_if_not_installed("fgsea")

  abundance <- data.frame(
    `#NAME` = paste0("K", sprintf("%05d", 1:8)),
    S1 = c(1:8),
    S2 = c(2:9),
    S3 = c(10:17),
    S4 = c(11:18),
    check.names = FALSE
  )
  metadata <- data.frame(
    sample = paste0("S", 1:4),
    group = c("A", "A", "B", "B"),
    stringsAsFactors = FALSE
  )

  local_mocked_bindings(
    prepare_gene_sets = function(...) {
      list("ko00010" = abundance$`#NAME`)
    },
    run_fgsea = function(...) {
      data.frame(
        pathway_id = "ko00010",
        pathway_name = "ko00010",
        size = 8,
        ES = 0.5,
        NES = 1.2,
        pvalue = 0.01,
        p.adjust = 0.02,
        leading_edge = "K00001",
        stringsAsFactors = FALSE
      )
    }
  )

  result <- pathway_gsea(
    abundance = abundance,
    metadata = metadata,
    group = "group",
    method = "fgsea"
  )
  expect_s3_class(result, "data.frame")
  expect_equal(result$pathway_id, "ko00010")
})

test_that("prepare_gene_sets warns when `organism` is set to a non-default value but still returns the same KO gene sets", {
  # Regression: `organism` was documented on both `prepare_gene_sets()`
  # and `pathway_gsea()` but the KEGG and GO branches both read
  # KO-based reference tables (`ko_to_kegg_reference`,
  # `ko_to_go_reference`) that are organism-independent by
  # construction. A user passing `organism = "hsa"` silently got the
  # same gene sets as the default -- a promise/implementation gap --
  # which we now surface as a deprecation warning. The gene-set
  # content must remain identical across organism values (the
  # argument truly has no effect), both for KEGG and GO.
  expect_warning(
    gs_hsa <- prepare_gene_sets("KEGG", organism = "hsa"),
    regexp = "organism.*deprecated.*no effect"
  )
  gs_default <- suppressWarnings(prepare_gene_sets("KEGG", organism = "ko"))
  expect_identical(gs_hsa, gs_default)

  # Default value does not warn.
  expect_warning(
    prepare_gene_sets("KEGG", organism = "ko"),
    regexp = NA
  )

  # GO branch: same contract.
  expect_warning(
    go_hsa <- prepare_gene_sets("GO", organism = "hsa"),
    regexp = "organism.*deprecated.*no effect"
  )
  go_default <- suppressWarnings(prepare_gene_sets("GO", organism = "ko"))
  expect_identical(go_hsa, go_default)
})