Fixes a bug in binomial testing

as well as adding a test for it, and then fixing a test in gocats importing

DESCRIPTION

@@ -1,5 +1,5 @@
 Package: categoryCompare2
-Version: 0.100.15
+Version: 0.100.16
 Title: Meta-Analysis of High-Throughput Experiments Using Feature
         Annotations
 Author: Robert M. Flight <[email protected]>

R/binomial_enrichment.R

@@ -49,6 +49,11 @@ binomial_feature_enrichment = function(binomial_features,
   n_feature = sapply(tmp_annot_feature, length)
   keep_annot = n_feature >= min_features
 
+  if (sum(keep_annot) == 0) {
+    warning("No annotations had more than min_features, no tests ran!")
+    return(NULL)
+  }
+
   tmp_annot_feature = tmp_annot_feature[keep_annot]
 
   binomial_features@annotation@annotation_features = tmp_annot_feature

R/gocats.R

@@ -44,7 +44,7 @@ gocats_to_annotation = function(ancestors_file = "ancestors.json",
 
   if (!is.null(namespace_file)) {
     if (!file.exists(namespace_file)) {
-      message(paste0(namespace_file, " does not exist. GO namespace will not be updated."))
+      warning(paste0(namespace_file, " does not exist. GO namespace will not be updated."))
       namespaces_short = character(0)
     } else {
       namespaces = jsonlite::fromJSON(namespace_file) |> unlist()

inst/extdata/generate_test_data.R

@@ -20,6 +20,10 @@ estrogen_raw <- ReadAffy(filenames = file.path(data_dir, rownames(pData(pheno_da
 
 estrogen_rma <- rma(estrogen_raw)
 
+affy_ids = featureNames(estrogen_rma)
+affy_other = AnnotationDbi::select(hgu95av2.db, keys = affy_ids, columns = c("ENTREZID", "SYMBOL", "UNIPROT")) |>
+  dplyr::mutate(affy = PROBEID)
+
 e10 <- estrogen_rma[, estrogen_rma$time.h == 10]
 e10 <- nsFilter(e10, remove.dupEntrez=TRUE, var.filter=FALSE, 
                 feature.exclude="^AFFX")$eset
@@ -33,8 +37,7 @@ fit10_2 <- contrasts.fit(fit10, c10)
 eB10 <- eBayes(fit10_2)
 table10 <- topTable(eB10, number=nrow(e10), p.value=1, adjust.method="BH")
 table10$affy <- rownames(table10)
-table10$Entrez <- AnnotationDbi::select(hgu95av2.db, keys = table10$affy, columns = "ENTREZID")[, "ENTREZID"]
-table10$Symbol <- AnnotationDbi::select(hgu95av2.db, keys = table10$affy, columns = "SYMBOL")[, "SYMBOL"]
+table10 = dplyr::left_join(table10, affy_other, by = "affy")
 
 e48 <- estrogen_rma[, estrogen_rma$time.h == 48]
 e48 <- nsFilter(e48, remove.dupEntrez=TRUE, var.filter=FALSE, 
@@ -49,24 +52,31 @@ fit48_2 <- contrasts.fit(fit48, c48)
 eB48 <- eBayes(fit48_2)
 table48 <- topTable(eB48, number=nrow(e48), p.value=1, adjust.method="BH")
 table48$affy <- rownames(table48)
-table48$Entrez <- AnnotationDbi::select(hgu95av2.db, keys = table48$affy, columns = "ENTREZID")[, "ENTREZID"]
-table48$Symbol <- AnnotationDbi::select(hgu95av2.db, keys = table48$affy, columns = "SYMBOL")[, "SYMBOL"]
 
+table48 = dplyr::left_join(table48, affy_other, by = "affy")
 adj_cut <- 0.05
 
-t48_entrez <- dplyr::filter(table48, adj.P.Val <= adj_cut) %>% extract2("Entrez")
-t10_entrez <- dplyr::filter(table10, adj.P.Val <= adj_cut) %>% extract2("Entrez")
-universe_entrez <- table10$Entrez
+t48_entrez <- dplyr::filter(table48, adj.P.Val <= adj_cut) %>% dplyr::pull("ENTREZID") %>% unique()
+t10_entrez <- dplyr::filter(table10, adj.P.Val <= adj_cut) %>% dplyr::pull("ENTREZID") %>% unique()
+universe_entrez <- unique(table10$ENTREZID)
 
 cat(t48_entrez, sep = "\n", file = "executable_related/test_data/48_entrez.txt")
 cat(t10_entrez, sep = "\n", file = "executable_related/test_data/10_entrez.txt")
 cat(universe_entrez, sep = "\n", file = "executable_related/test_data/universe_entrez.txt")
 
-t48_symbol <- dplyr::filter(table48, adj.P.Val <= adj_cut) %>% extract2("Symbol")
-t10_symbol <- dplyr::filter(table10, adj.P.Val <= adj_cut) %>% extract2("Symbol")
-universe_symbol <- table10$Symbol
+t48_symbol <- dplyr::filter(table48, adj.P.Val <= adj_cut) %>% extract2("SYMBOL") %>% unique()
+t10_symbol <- dplyr::filter(table10, adj.P.Val <= adj_cut) %>% extract2("SYMBOL") %>% unique()
+universe_symbol <- unique(table10$SYMBOL)
 
 cat(t48_symbol, sep = "\n", file = "executable_related/test_data/48_symbol.txt")
 cat(t10_symbol, sep = "\n", file = "executable_related/test_data/10_symbol.txt")
 cat(universe_symbol, sep = "\n", file = "executable_related/test_data/universe_symbol.txt")
 
+t48_uniprot <- dplyr::filter(table48, adj.P.Val <= adj_cut) %>% extract2("UNIPROT") %>% unique()
+t10_uniprot <- dplyr::filter(table10, adj.P.Val <= adj_cut) %>% extract2("UNIPROT") %>% unique()
+universe_uniprot <- unique(table10$UNIPROT)
+
+cat(t48_uniprot, sep = "\n", file = "executable_related/test_data/48_uniprot.txt")
+cat(t10_uniprot, sep = "\n", file = "executable_related/test_data/10_uniprot.txt")
+cat(universe_uniprot, sep = "\n", file = "executable_related/test_data/universe_uniprot.txt")
+

tests/testthat/test-binomial_enrichment.R

@@ -69,3 +69,11 @@ test_that("binomial greater works", {
 
   expect_snapshot_value(binomial_results, style = "serialize")
 })
+
+
+test_that("binomial warning works", {
+  lipid_features = readRDS("lipid_binomial_testing.rds")
+  expect_warning(binomial_feature_enrichment(lipid_features, min_features = 6), "No annotations had more than")
+
+  expect_no_warning(binomial_feature_enrichment(lipid_features, min_features = 1))
+})

tests/testthat/test-gocats.R

@@ -12,7 +12,7 @@ test_that("gocats annotation importing works", {
   list_translation = as.list(ensembl_uniprot)
 
   expect_error(gocats_to_annotation("random_file"))
-  expect_warning(gocats_to_annotation(ancestors_file, "random_file"), "does not exist")
+  expect_warning(gocats_to_annotation(ancestors_file, "random_file"), "random_file does not exist")
   expect_error(gocats_to_annotation(ancestors_file, feature_translation = list_translation), "must be a data.frame")
   expect_error(gocats_to_annotation(ancestors_file, feature_translation = bad_translation), "must contain the columns")