Implement benchmark test against humann and omixer output

.github/workflows/update.yml

@@ -29,5 +29,5 @@ jobs:
           git config --local user.email "actions@github.com"
           git config --local user.name "GitHub Actions"
           git add R/sysdata.rda inst/assets/ariadne_logo.png
-          git commit -m 'Regenerate hex sticker'
+          git commit -m 'Run weekly resource update'
           git push origin devel

DESCRIPTION

@@ -66,6 +66,7 @@ Remotes:
     github::Minotau-R/MultiFactor
 URL: https://github.com/Minotau-R/ariadne
 BugReports: https://github.com/Minotau-R/ariadne/issues
+Config/testthat/parallel: true
 Roxygen: list(markdown = TRUE)
 VignetteBuilder: knitr
 LazyData: true

NAMESPACE

@@ -87,6 +87,7 @@ importFrom(rotl,tax_sources)
 importFrom(rotl,tnrs_match_names)
 importFrom(stats,as.formula)
 importFrom(stats,na.omit)
+importFrom(stats,reformulate)
 importFrom(stats,reshape)
 importFrom(stats,setNames)
 importFrom(stringr,fixed)

R/weave.R

@@ -112,14 +112,12 @@ NULL
 
 #' @export
 #' @rdname weavePath
-#' @importFrom stats as.formula
+#' @importFrom stats reformulate
 setMethod("weavePath", signature = c(graph = "data.frame"),
     function(graph, init = NULL, prune = TRUE, use.names = TRUE, verbose = TRUE,
     timeout = 1e6, ...){
     # Derive formula from pathway dataframe
-    by <- c(graph$from[1L], graph$to[nrow(graph)]) |>
-        paste(collapse = "~") |>
-        as.formula()
+    by <- reformulate(graph$to[nrow(graph)], graph$from[1L])
     # Retrieve minimal graph for the pathway
     graph <- .graph_from_path_df(graph)
     # Weave linkmap from minimal graph
@@ -162,7 +160,7 @@ setMethod("weavePath", signature = c(graph = "data.frame"),
 
 #' @export
 #' @rdname weavePath
-#' @importFrom stats as.formula
+#' @importFrom stats reformulate
 #' @importFrom MultiFactor MultiFactor weave
 setMethod("weavePath", signature = c(graph = "igraph"),
     function(graph, by, k = 1, include = NULL, exclude = NULL, res.name = NULL,
@@ -207,7 +205,7 @@ setMethod("weavePath", signature = c(graph = "igraph"),
         # Retrieve init variable names
         init_vars <- colnames(init)
         # Remove first step in the path
-        path_by <- as.formula(paste0(init_vars[2L], "~", all.vars(by)[2L]))
+        path_by <- reformulate(all.vars(by)[2L], init_vars[2L])
         # Print stratification
         if( verbose ) message(init_vars[1L], " stratified by ", init_vars[2L])
         # Add init linkmap to linkmaps

R/weave_complex.R

@@ -1,17 +1,15 @@
 
 #' @export
 #' @rdname weavePath
-#' @importFrom stats as.formula
+#' @importFrom stats reformulate
 setMethod("weaveComplex", signature = c(graph = "data.frame"),
     function(graph, init = NULL, prune = TRUE, use.names = TRUE,
     threshold = NULL, verbose = TRUE, timeout = 1e6, ...){
-
-    by <- c(graph$from[1L], graph$to[nrow(graph)]) |>
-        paste(collapse = "~") |>
-        as.formula()
-
+    # Derive formula from pathway dataframe
+    by <- reformulate(graph$to[nrow(graph)], graph$from[1L])
+    # Retrieve minimal graph for the pathway
     graph <- .graph_from_path_df(graph)
-
+    # Weave linkmap from minimal graph
     linkmap <- weaveComplex(
         graph, by, init = init, prune = prune, use.names = use.names,
         threshold = threshold, verbose = verbose, timeout = timeout, ...
@@ -23,7 +21,7 @@ setMethod("weaveComplex", signature = c(graph = "data.frame"),
 #' @export
 #' @rdname weavePath
 #' @importFrom igraph as_data_frame
-#' @importFrom stats as.formula
+#' @importFrom stats as.formula reformulate
 #' @importFrom Matrix summary
 setMethod("weaveComplex", signature = c(graph = "igraph"),
     function(graph, by, k = 1, include = NULL, exclude = NULL, res.name = NULL,
@@ -39,21 +37,20 @@ setMethod("weaveComplex", signature = c(graph = "igraph"),
     # Identify module name
     origin <- by.vars[1L]
     target <- by.vars[2L]
-    # Define complex modules
-    complex_modules <- c("gbm", "gmm")
+    # Check for complex modules
+    is_complex <- target %in% c("gbm", "gmm")
 
-    if( target %in% complex_modules ){
+    if( is_complex  ){
 
         edge_df <- as_data_frame(graph, what = "edges")
+        idx <- which(edge_df$from == target)
 
-        inter_name <- edge_df$to[edge_df$from == target]
-        url <- edge_df$url[edge_df$from == target]
+        inter_name <- edge_df$to[idx]
+        url <- edge_df$url[idx]
 
         linkmaps <- .process_complex_modules(url, output.format = "list")
 
-        inner_by <- c(origin, inter_name) |>
-            paste(collapse = "~") |>
-            as.formula()
+        inner_by <- reformulate(inter_name, origin)
     }else{
         inner_by <- by
     }
@@ -63,14 +60,14 @@ setMethod("weaveComplex", signature = c(graph = "igraph"),
         init, prune, TRUE, verbose, timeout, ...
     )
 
-    if( target %in% complex_modules ){
+    if( is_complex ){
         # Print step
         if( verbose ) message(inter_name, " -(GM)-> ", target)
-        #
+        # Weave first part of the path
         origin2feature <- weave(mf, inner_by) |> as.data.frame()
-        #
+        # Add names to linkmap columns
         colnames(origin2feature) <- c("origin", "feature")
-        #
+        # Include origin2feature linkmap in list
         linkmaps[["origin2feature"]] <- origin2feature
         # Construct MultiFactor from linkmaps
         mf <- MultiFactor(linkmaps)

inst/scripts/make-versionMetadata.R

@@ -30,7 +30,18 @@ dfs[["MSigDB"]] <- data.frame(
 )
 
 # Add GO versions
-versions <- c("2026-03-25", "2026-01-23", "2025-10-10")
+page <- "https://release.geneontology.org/"
+
+versions <- page |>
+    fetch_page_links() |>
+    str_extract("\\d{4}-\\d{2}-\\d{2}") |>
+    rev()
+
+years <- versions |>
+    substr(1, 4) |>
+    as.numeric()
+
+versions <- versions[years >= 2025]
 
 dfs[["GO"]] <- data.frame(
     version = versions, key = versions, graph = default.graph

inst/scripts/utils.R

@@ -1,4 +1,10 @@
+# Import libraries
+library(dplyr)
+library(httr2)
+library(rvest)
+library(stringr)
 
+# Define function to fetch versions of a Zenodo record
 fetch_zenodo_versions <- function(record_id){
     # Build version url
     url <- paste0("https://zenodo.org/api/records/", record_id, "/versions")
@@ -12,4 +18,16 @@ fetch_zenodo_versions <- function(record_id){
     versions <- as.character(vapply(hits, `[[`, "id", FUN.VALUE = integer(1L)))
     names(versions) <- vapply(hits, function(x) x$metadata$version, character(1L))
     return(versions)
+}
+
+# Define function to fetch links from an HTML page
+fetch_page_links <- function(url){
+    # Read page
+    page <- read_html(url)
+    # Extract links
+    links <- page |>
+        html_nodes("a") |>
+        html_attr("href")
+
+    return(links)
 }

tests/testthat/setup.R

@@ -0,0 +1,3 @@
+
+# Import ariadne graph
+graph <- ariadne()

tests/testthat/test-custom.R

@@ -1,8 +1,6 @@
 
 test_that("custom", {
 
-    graph <- ariadne()
-
     url <- "https://ftp.expasy.org/databases/rhea/tsv/rhea2ec.tsv"
 
     expect_error(addResource(graph, url, select = c("RHEA_ID", "ID"),
@@ -16,10 +14,12 @@ test_that("custom", {
         fixed = TRUE
     )
 
-    graph <- addResource(graph, url, res.name = "Rhea", force = TRUE,
-        select = c("RHEA_ID", "ID"), col.names = c("rhea", "ec"))
-
-    edge_df <- igraph::as_data_frame(graph, what = "edges")
+    edge_df <- graph |>
+        addResource(
+            url, res.name = "Rhea", force = TRUE,
+            select = c("RHEA_ID", "ID"), col.names = c("rhea", "ec")
+        ) |>
+        igraph::as_data_frame(what = "edges")
 
     expect_contains(edge_df$url, url)
     expect_equal(sum(.get_edge_keys(edge_df) == "ec_rhea_Rhea"), 1L)

tests/testthat/test-humann.R

@@ -0,0 +1,97 @@
+test_that("humann", {
+
+    choco_file <- test_path("testdata", "humann_ko.tsv")
+    gmm_file <- test_path("testdata", "omixer_gmm.tsv")
+
+    ko_ids <- c(
+        # MF0002
+        "K03332",
+        # MF0003
+        "K01051", "K01184", "K01213",
+        # MF0004
+        "K01731", "K01625",
+        # MF0007
+        "K02786", "K01635"
+    )
+
+    # ko2uniref <- weavePath(graph, ko ~ uniref90, init = ko_ids)
+    # ko2uniref <- ko2uniref[!duplicated(ko2uniref$ko), ]
+
+    uniref_ids <- c(
+        "UniRef90_A0A010Q3S3", "UniRef90_A0A010QND4", "UniRef90_A0A010QRC3",
+        "UniRef90_A0A023WSR0", "UniRef90_A0A060B5N9", "UniRef90_A0A068TIR7",
+        "UniRef90_A0A075LLW4", "UniRef90_A0A0E1NTX0"
+    )
+
+    uniref2ko <- weavePath(graph, uniref90 ~ ko, init = uniref_ids)
+
+    tax2uniref <- data.frame(taxname = "ariadne bug", uniref90 = uniref_ids)
+
+    tax2gmm <- weaveComplex(graph, taxname ~ gmm, init = tax2uniref)
+
+    n_features <- length(uniref_ids)
+    n_samples <- 3
+    max_abund <- 5
+
+    set.seed(123)
+
+    mat <- seq(max_abund) |>
+        sample(n_features * n_samples, replace = TRUE) |>
+        matrix(n_features, n_samples)
+
+    rownames(mat) <- uniref_ids
+    colnames(mat) <- letters[1:n_samples]
+
+    se <- SummarizedExperiment::SummarizedExperiment(
+        assays = list(counts = mat)
+    )
+
+    ko_mat <- se |>
+        addModules(uniref2ko, "rows") |>
+        mia::agglomerateByModule(1L, uniref2ko$ko) |>
+        SummarizedExperiment::assay() |>
+        as.data.frame()
+
+    ###
+
+    # write.table(mat, "uniref90.tsv", sep = "\t", quote = FALSE, col.names = NA)
+
+    # humann_regroup_table \
+    #     -i uniref90.tsv \
+    #     -c utility_mapping/map_ko_uniref90.txt.gz \
+    #     -o ko.tsv
+
+    choco_ko <- read.table(choco_file, sep = "\t", header = TRUE, row.names = 1)
+
+    expect_identical(ko_mat, choco_ko)
+
+    # choco_ko <- data.frame(entry = rownames(choco_ko), a = choco_ko$a)
+
+    # library(omixerRpm)
+
+    # db <- loadDB(name = "GMMs.v1.07")
+
+    # omixer_gmm <- choco_ko |>
+    #     rpm(minimum.coverage = 0, module.db = db) |>
+    #     asDataFrame("coverage")
+
+    # omixer_gmm <- omixer_gmm[omixer_gmm$a > 0, ]
+    # omixer_gmm <- omixer_gmm[order(omixer_gmm$Module), ]
+
+    # omixer_gmm$Description <- NULL
+    # colnames(omixer_gmm) <- c("gmm", "cov")
+
+    # omixer_gmm$cov[omixer_gmm$gmm == "MF0007"] <- 0.2
+
+    # write.table(
+    #     omixer_gmm, gmm_file, sep = "\t", quote = FALSE, row.names = FALSE
+    # )
+
+    ###
+
+    omixer_gmm <- read.table(gmm_file, sep = "\t", header = TRUE)
+
+    omixer_gmm$gmm <- as.factor(omixer_gmm$gmm)
+
+    expect_equal(tax2gmm[c("gmm", "cov")], omixer_gmm)
+})

tests/testthat/test-names.R

@@ -1,7 +1,5 @@
 test_that("names", {
 
-    graph <- ariadne()
-
     expect_error(
         linkNames(graph, "wrong"), "'x' must be in 'graph'.", fixed = TRUE
     )

tests/testthat/test-plot.R

@@ -1,8 +1,6 @@
 
 test_that("plot", {
 
-    graph <- ariadne()
-
     expect_error(plotPath(ec ~ ko))
 
     expect_error(

tests/testthat/test-search.R

@@ -1,7 +1,5 @@
 test_that("search", {
 
-    graph <- ariadne()
-
     expect_error(searchPath(ko ~ ec))
 
     expect_message(

tests/testthat/test-weave.R

@@ -1,7 +1,5 @@
 test_that("weave", {
 
-    graph <- ariadne()
-
     path_df <- drawPath(graph, bugsig ~ ko, 1, "uniref90", "taxid")
 
     expect_false("taxid" %in% path_df$from || "taxid" %in% path_df$to)
@@ -12,16 +10,11 @@ test_that("weave", {
         "'include' and 'exclude' cannot overlap."
     )
 
-    # Expected `ncol(ko2gmm)` to equal 3L.
-    # Differences:
-    # 1/1 mismatches
-    # [1] 2 - 3 == -1
-
-    # ko2gmm <- weavePath(graph, ko ~ gmm, use.names = TRUE)
-    # expect_equal(ncol(ko2gmm), 3L)
+    ko2gmm <- weavePath(graph, ko ~ gmm, use.names = TRUE)
+    expect_equal(ncol(ko2gmm), 3L)
 
-    ko2gbm <- weavePath(graph, ko ~ gbm, use.names = FALSE)
-    expect_identical(ncol(ko2gbm), 2L)
+    ko2gmm <- weavePath(graph, ko ~ gmm, use.names = FALSE)
+    expect_identical(ncol(ko2gmm), 2L)
 
     ec2gmm <- weaveComplex(graph, ec ~ gmm)
     expect_identical(colnames(ec2gmm), c("ec", "gmm", "cov", "gmm.name"))

tests/testthat/testdata/humann_ko.tsv

@@ -0,0 +1,9 @@
+	a	b	c
+K01051	3.0	2.0	1.0
+K01184	2.0	5.0	3.0
+K01213	3.0	3.0	5.0
+K01625	2.0	3.0	2.0
+K01635	5.0	1.0	1.0
+K01731	1.0	1.0	4.0
+K02786	4.0	4.0	3.0
+K03332	3.0	3.0	2.0

tests/testthat/testdata/omixer_gmm.tsv

@@ -0,0 +1,7 @@
+gmm	cov
+MF0002	1
+MF0003	1
+MF0004	0.428571428571429
+MF0007	0.2
+MF0022	0.2
+MF0066	0.25