0.3.5

Author: traversc

ChangeLog

@@ -1,8 +1,8 @@
-Version 0.3.3 (2025-10-04)
-   * (Performance) Implement Ukkonen style banding to limit DP calculations to near diagonals
+Version 0.3.5 (2025-10-19)
+   * (Performance) Specialized fast `single_gap_search` algorithm, see `RadixTree$single_gap_search`
 
-Version 0.3.2 (2025-10-03)
-   * (Performance) Track DP column min to avoid redundant calculation
+Version 0.3.3 (2025-10-04)
+   * (Performance) Implement Ukkonen style banding to limit DP calculations to near diagonals, improves performance 30-40%
 
 Version 0.3.1 (2025-10-03)
    * (Performance) Reworked trie alignment searches to reuse dynamic-programming workspaces
@@ -24,3 +24,4 @@ Version 0.2.7 (2024-01-25)
 
 Version 0.2.6 (2024-01-25)
    * Minor fix for CRAN checks, std::basic_string<uint8_t> deprecated in Clang 18
+ 

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: seqtrie
 Title: Radix Tree and Trie-Based String Distances
-Version: 0.3.4
-Date: 2025-10-04
+Version: 0.3.5
+Date: 2025-10-19
 Authors@R: c(
     person("Travers", "Ching", email = "traversc@gmail.com", role = c("aut", "cre", "cph")),
     person("Martin", "Moene", role = c("ctb", "cph"), comment = "span-lite C++ library"),
@@ -20,7 +20,7 @@ LinkingTo:
 Imports: 
     Rcpp (>= 0.12.18.3), RcppParallel (>= 5.1.3), R6, rlang, dplyr, stringi
 Suggests:
-    knitr, rmarkdown, stringdist, qs, pwalign, igraph, ggplot2
+    knitr, rmarkdown, stringdist, pwalign, igraph, ggplot2
 VignetteBuilder: knitr
 RoxygenNote: 7.3.2
 Roxygen: list(markdown = TRUE)

Makefile

@@ -74,6 +74,7 @@ test:
 	IS_LOCAL=Yes Rscript tests/test_pairwise.R && unset IS_LOCAL
 	IS_LOCAL=Yes Rscript tests/test_RadixTree.R && unset IS_LOCAL
 	IS_LOCAL=Yes Rscript tests/test_RadixForest.R && unset IS_LOCAL
+	IS_LOCAL=Yes Rscript tests/test_single_gap_search.R && unset IS_LOCAL
 
 test-trie:
 	IS_LOCAL=Yes Rscript tests/test_RadixTree.R && unset IS_LOCAL

R/RadixForest.R

@@ -141,9 +141,7 @@ RadixForest <- R6::R6Class("RadixForest", list(
     }
 
     if(!is.null(max_distance)) {
-      if(length(max_distance) == 1) {
-        max_distance <- rep(max_distance, length(query))
-      }
+      max_distance <- recycle_arg(max_distance, query)
     } else if(!is.null(max_fraction)) {
       max_distance <- as.integer(nchar(query) * max_fraction)
     } else {

R/RadixTree.R

@@ -141,15 +141,13 @@ RadixTree <- R6::R6Class("RadixTree", public = list(
   #' @param show_progress `r rdoc("show_progress")`
   #' @return The output is a data.frame of all matches with columns "query" and "target".
   #' For anchored searches, the output also includes attributes "query_size" and "target_size" which are vectors containing the portion of the query and target sequences that are aligned.
-  search = function(query, max_distance = NULL, max_fraction = NULL, mode = "levenshtein", cost_matrix = NULL, gap_cost = NULL, gap_open_cost = NULL, nthreads = 1, show_progress = FALSE) {
+  search = function(query, max_distance = NULL, max_fraction = NULL, mode = "levenshtein", cost_matrix = NULL, gap_cost = NA_integer_, gap_open_cost = NA_integer_, nthreads = 1, show_progress = FALSE) {
     charset <- unique(c(CharCounter_keys(self$char_counter_pointer), get_charset(query)))
     check_alignment_params(mode, cost_matrix, gap_cost, gap_open_cost, charset, diag_must_be_zero = TRUE)
     mode <- normalize_mode_parameter(mode)
 
     if (!is.null(max_distance)) {
-      if (length(max_distance) == 1) {
-        max_distance <- rep(max_distance, length(query))
-      }
+      max_distance <- recycle_arg(max_distance, query)
     } else if (!is.null(max_fraction)) {
       max_distance <- as.integer(nchar(query) * max_fraction)
     } else {
@@ -159,14 +157,12 @@ RadixTree <- R6::R6Class("RadixTree", public = list(
       stop("max_distance/max_fraction must be non-negative")
     }
 
-    # defaults for C++ plain ints
-    if (is.null(gap_cost)) gap_cost <- 1L
-    if (is.null(gap_open_cost)) gap_open_cost <- 0L
-    # Align conventions with pwalign/Biostrings: first gap includes one extension
-    if (gap_open_cost > 0L) gap_open_cost <- gap_open_cost + gap_cost
-    result <- RadixTree_search(self$root_pointer, query, max_distance, mode, cost_matrix, as.integer(gap_cost), as.integer(gap_open_cost), nthreads, show_progress)
+    if (!is.na(gap_open_cost) && !is.na(gap_cost) && gap_open_cost > 0L) {
+      gap_open_cost <- gap_open_cost + gap_cost
+    }
+    result <- RadixTree_search(self$root_pointer, query, max_distance, mode, cost_matrix, gap_cost, gap_open_cost, nthreads, show_progress)
     if (mode == "anchored") { # Append query_size and target_size attributes
-      result2 <- c_dist_pairwise(result$query, result$target, mode, cost_matrix, as.integer(gap_cost), as.integer(gap_open_cost), nthreads, show_progress = FALSE)
+      result2 <- c_dist_pairwise(result$query, result$target, mode, cost_matrix, gap_cost, gap_open_cost, nthreads, show_progress = FALSE)
       if (any(result$distance != result2)) {
         stop("Internal error: anchored search results do not match pairwise results")
       }
@@ -175,6 +171,29 @@ RadixTree <- R6::R6Class("RadixTree", public = list(
     }
     result
   },
+  #' @description A specialized algorithm for searching for sequences allowing at most a single gap within the alignment itself. The mode is always "anchored" and does not penalize end gaps.
+  #' @param query `r rdoc("query")`
+  #' @param max_distance `r rdoc("max_distance")`
+  #' @param gap_cost `r rdoc("gap_cost")`
+  #' @param nthreads `r rdoc("nthreads")`
+  #' @param show_progress `r rdoc("show_progress")`
+  #' @return The output is a data.frame of matches with columns "query", "target" and "distance".
+  single_gap_search = function(query,
+                               max_distance,
+                               gap_cost = 1L,
+                               nthreads = 1,
+                               show_progress = FALSE) {
+    if (!is.null(max_distance)) {
+      max_distance <- recycle_arg(max_distance, query)
+    } else {
+      stop("Either max_distance must be non-null")
+    }
+    if (any(max_distance < 0)) {
+      stop("max_distance/max_fraction must be non-negative")
+    }
+
+    RadixTree_single_gap_search(self$root_pointer, query, max_distance, gap_cost, nthreads, show_progress)
+  },
   #' @description Validate the tree
   #' @return A logical indicating whether the tree is valid (TRUE) or not (FALSE). This is mostly an internal function for debugging purposes and should always return TRUE.
   validate = function() {

R/RadixTree_search_helpers.R

@@ -23,7 +23,7 @@
 #' dist_search(c("ACGT", "AAAA"), c("ACG", "ACGT"), max_distance = 1, mode = "levenshtein")
 #' @name dist_search
 dist_search <- function(query, target, max_distance = NULL, max_fraction = NULL, mode = "levenshtein",
-                        cost_matrix = NULL, gap_cost = NULL, gap_open_cost = NULL, tree_class = "RadixTree",
+                        cost_matrix = NULL, gap_cost = NA_integer_, gap_open_cost = NA_integer_, tree_class = "RadixTree",
                         nthreads = 1, show_progress = FALSE) {
   if (!tree_class %in% c("RadixTree", "RadixForest")) {
     stop("tree_class must be one of RadixTree or RadixForest")
@@ -33,7 +33,9 @@ dist_search <- function(query, target, max_distance = NULL, max_fraction = NULL,
     obj$insert(target)
     obj$search(query, max_distance, max_fraction, mode, cost_matrix, gap_cost, gap_open_cost, nthreads, show_progress)
   } else if(tree_class == "RadixForest") {
-    if(!is.null(cost_matrix) || !is.null(gap_cost) || !is.null(gap_open_cost)) {
+    gap_cost_provided <- !(length(gap_cost) == 1L && is.na(gap_cost))
+    gap_open_provided <- !(length(gap_open_cost) == 1L && is.na(gap_open_cost))
+    if(!is.null(cost_matrix) || gap_cost_provided || gap_open_provided) {
       stop("cost_matrix, gap_cost and gap_open_cost are not supported for RadixForest")
     }
     obj <- RadixForest$new()
@@ -117,4 +119,4 @@ split_search <- function(query, target, query_split, target_split, edge_trim = 0
   results <- dplyr::mutate(results, distance = .data$distance.left + .data$distance.right)
   results <- dplyr::filter(results, .data$distance <= max_distance)
   as.data.frame(results[c("query", "target", "distance")])
-}
+}

R/RcppExports.R

@@ -105,15 +105,19 @@ RadixTree_create <- function() {
     .Call(`_seqtrie_RadixTree_create`)
 }
 
-RadixTree_search <- function(xp, query, max_distance, mode = "global", cost_matrix = NULL, gap_cost = 1L, gap_open_cost = 0L, nthreads = 1L, show_progress = FALSE) {
+RadixTree_search <- function(xp, query, max_distance, mode = "global", cost_matrix = NULL, gap_cost = NA_integer_, gap_open_cost = NA_integer_, nthreads = 1L, show_progress = FALSE) {
     .Call(`_seqtrie_RadixTree_search`, xp, query, max_distance, mode, cost_matrix, gap_cost, gap_open_cost, nthreads, show_progress)
 }
 
-c_dist_matrix <- function(query, target, mode = "global", cost_matrix = NULL, gap_cost = 1L, gap_open_cost = 0L, nthreads = 1L, show_progress = FALSE) {
+RadixTree_single_gap_search <- function(xp, query, max_distance, gap_cost = 1L, nthreads = 1L, show_progress = FALSE) {
+    .Call(`_seqtrie_RadixTree_single_gap_search`, xp, query, max_distance, gap_cost, nthreads, show_progress)
+}
+
+c_dist_matrix <- function(query, target, mode = "global", cost_matrix = NULL, gap_cost = NA_integer_, gap_open_cost = NA_integer_, nthreads = 1L, show_progress = FALSE) {
     .Call(`_seqtrie_c_dist_matrix`, query, target, mode, cost_matrix, gap_cost, gap_open_cost, nthreads, show_progress)
 }
 
-c_dist_pairwise <- function(query, target, mode = "levenshtein", cost_matrix = NULL, gap_cost = 1L, gap_open_cost = 0L, nthreads = 1L, show_progress = FALSE) {
+c_dist_pairwise <- function(query, target, mode = "levenshtein", cost_matrix = NULL, gap_cost = NA_integer_, gap_open_cost = NA_integer_, nthreads = 1L, show_progress = FALSE) {
     .Call(`_seqtrie_c_dist_pairwise`, query, target, mode, cost_matrix, gap_cost, gap_open_cost, nthreads, show_progress)
 }
 

R/pairwise.R

@@ -18,16 +18,14 @@
 #' @examples
 #' dist_matrix(c("ACGT", "AAAA"), c("ACG", "ACGT"), mode = "global")
 #' @name dist_matrix
-dist_matrix <- function(query, target, mode, cost_matrix = NULL, gap_cost = NULL, gap_open_cost = NULL, nthreads = 1, show_progress = FALSE) {
+dist_matrix <- function(query, target, mode, cost_matrix = NULL, gap_cost = NA_integer_, gap_open_cost = NA_integer_, nthreads = 1, show_progress = FALSE) {
   charset <- unique(c(get_charset(query), get_charset(target)))
   check_alignment_params(mode, cost_matrix, gap_cost, gap_open_cost, charset, diag_must_be_zero = FALSE)
   mode <- normalize_mode_parameter(mode)
-  # defaults for C++ plain ints
-  if (is.null(gap_cost)) gap_cost <- 1L
-  if (is.null(gap_open_cost)) gap_open_cost <- 0L
-  # Align conventions with pwalign/Biostrings: first gap includes one extension
-  if (gap_open_cost > 0L) gap_open_cost <- gap_open_cost + gap_cost
-  c_dist_matrix(query, target, mode, cost_matrix, as.integer(gap_cost), as.integer(gap_open_cost), nthreads, show_progress)
+  if (!is.na(gap_open_cost) && !is.na(gap_cost) && gap_open_cost > 0L) {
+    gap_open_cost <- gap_open_cost + gap_cost
+  }
+  c_dist_matrix(query, target, mode, cost_matrix, gap_cost, gap_open_cost, nthreads, show_progress)
 }
 
 #' @title Pairwise distance between two sets of sequences
@@ -50,17 +48,15 @@ dist_matrix <- function(query, target, mode, cost_matrix = NULL, gap_cost = NULL
 #' @examples
 #' dist_pairwise(c("ACGT", "AAAA"), c("ACG", "ACGT"), mode = "global")
 #' @name dist_pairwise
-dist_pairwise <- function(query, target, mode, cost_matrix = NULL, gap_cost = NULL, gap_open_cost = NULL, nthreads = 1, show_progress = FALSE) {
+dist_pairwise <- function(query, target, mode, cost_matrix = NULL, gap_cost = NA_integer_, gap_open_cost = NA_integer_, nthreads = 1, show_progress = FALSE) {
   if (length(query) != length(target)) {
     stop("query and target must be the same length")
   }
   charset <- unique(c(get_charset(query), get_charset(target)))
   check_alignment_params(mode, cost_matrix, gap_cost, gap_open_cost, charset, diag_must_be_zero = FALSE)
   mode <- normalize_mode_parameter(mode)
-  # defaults for C++ plain ints
-  if (is.null(gap_cost)) gap_cost <- 1L
-  if (is.null(gap_open_cost)) gap_open_cost <- 0L
-  # Align conventions with pwalign/Biostrings: first gap includes one extension
-  if (gap_open_cost > 0L) gap_open_cost <- gap_open_cost + gap_cost
-  c_dist_pairwise(query, target, mode, cost_matrix, as.integer(gap_cost), as.integer(gap_open_cost), nthreads, show_progress)
+  if (!is.na(gap_open_cost) && !is.na(gap_cost) && gap_open_cost > 0L) {
+    gap_open_cost <- gap_open_cost + gap_cost
+  }
+  c_dist_pairwise(query, target, mode, cost_matrix, gap_cost, gap_open_cost, nthreads, show_progress)
 }

R/utils.R

@@ -58,6 +58,14 @@ is_integerlike <- function(x) {
   }
 }
 
+is_missing_arg <- function(x) {
+  is.null(x) || (length(x) == 1L && is.na(x))
+}
+
+recycle_arg <- function(recycle_vector, target_vector) {
+  rep(recycle_vector, length.out = length(target_vector))
+}
+
 # A cost matrix must
 # 1) Be a square matrix _OR_ NULL
 # 2) If gap_cost is NULL, cost_matrix must have a "gap" entry instead
@@ -103,16 +111,16 @@ check_cost_matrix <- function(cost_matrix, gap_cost, gap_open_cost, charset, dia
 
   # gap parameters are provided separately in the R APIs now
   # validate gap_cost / gap_open_cost here, independent of matrix contents
-  if (!is.null(gap_cost)) {
+  if (!is_missing_arg(gap_cost)) {
     if (!(is_integerlike(gap_cost) && length(gap_cost) == 1 && gap_cost > 0)) {
       stop("gap_cost must be a single positive integer when provided")
     }
   }
-  if (!is.null(gap_open_cost)) {
+  if (!is_missing_arg(gap_open_cost)) {
     if (!(is_integerlike(gap_open_cost) && length(gap_open_cost) == 1 && gap_open_cost >= 0)) {
       stop("gap_open_cost must be a single non-negative integer when provided")
     }
-    if (is.null(gap_cost)) {
+    if (is_missing_arg(gap_cost)) {
       stop("If gap_open_cost is defined, gap_cost must also be defined")
     }
   }
@@ -133,7 +141,7 @@ check_alignment_params <- function(mode, cost_matrix, gap_cost, gap_open_cost, c
     stop("mode must be one of hamming (hm), global (gb, lv, levenshtein) or anchored (an, en, extension)")
   }
   check_cost_matrix(cost_matrix, gap_cost, gap_open_cost, charset, diag_must_be_zero)
-  if ((mode %in% c("hamming", "hm")) && (!is.null(cost_matrix) || !is.null(gap_cost) || !is.null(gap_open_cost))) {
+  if ((mode %in% c("hamming", "hm")) && (!is.null(cost_matrix) || !is_missing_arg(gap_cost) || !is_missing_arg(gap_open_cost))) {
     warning("cost_matrix and gap parameters are ignored when mode is 'hamming'")
   }
 }

README.md

@@ -336,17 +336,6 @@ tree$prefix_search("car")
     ## 3   car carburetor
     ## 4   car  carbuncle
 
-### Why not just use Bowtie2, BWA or other fast alignment software?
-
-There are no apples-to-apples comparisons. With NGS alignment software,
-you are looking for alignments of reads (queries) *within* a genome
-reference (target). Here, we’re looking for alignments from the query to
-the *full* target. However, many NGS aligners do use Tries and similar
-data structures.
-
-Compared to pairwise alignment packages, calculating all alignment pairs
-takes much longer, but on the other hand gives you more information.
-
 ### References and literature
 
 - “Fast string correction with Levenshtein automata” (2002)