Refine networks

DESCRIPTION

@@ -1,13 +1,12 @@
 Package: immApex
 Title: Tools for Adaptive Immune Receptor Sequence-Based Machine and Deep Learning
-Version: 1.5.4
+Version: 1.7.1
 Authors@R: c(
     person(given = "Nick", family = "Borcherding", role = c("aut", "cre", "cph"), email = "ncborch@gmail.com"),
 	person(given = "Qile", family = "Yang", role = "ctb", email = "qile.yang@berkeley.edu", comment = c(ORCID = "0009-0005-0148-2499")))
 Description: A set of tools to for machine and deep learning in R from amino acid and nucleotide sequences focusing on adaptive immune receptors. The package includes pre-processing of sequences, unifying gene nomenclature usage, encoding sequences, and combining models. This package will serve as the basis of future immune receptor sequence functions/packages/models compatible with the scRepertoire ecosystem.
 License: MIT + file LICENSE 
 Encoding: UTF-8
-RoxygenNote: 7.3.3
 biocViews: Software, ImmunoOncology, SingleCell, Classification, Annotation, Sequencing, MotifAnnotation
 Depends: 
 	R (>= 4.3.0)
@@ -43,3 +42,4 @@ VignetteBuilder: knitr
 Language: en-US
 URL: https://github.com/BorchLab/immApex/
 BugReports: https://github.com/BorchLab/immApex/issues
+Config/roxygen2/version: 8.0.0

NEWS.md

@@ -1,3 +1,22 @@
+# immApex VERSION 1.7.1
+
+## PERFORMANCE
+`buildNetwork()` now scales to large, clonally expanded repertoires. The C++ engine was rewritten around candidate generation instead of all-pairs comparison, with no new dependencies and identical edge output for the default settings:
+
+* Identical sequences collapse to one representative before any distance is computed, so clonal expansion no longer drives quadratic work.
+* Length and V/J blocking act as an index. Only length-compatible representatives in the same block are compared.
+* Levenshtein and Damerau use deletion-neighborhood (SymSpell) hashing for small edit distances, with Damerau indexed at radius `2k`. Hamming uses pigeonhole segment indexing. Both fall back to length-blocked comparison outside their range.
+* Edges accumulate as integer node indices and gain their labels in R, which avoids building millions of barcode strings inside the parallel core.
+
+Speedups range from roughly 2x on mixed repertoires with many V genes to about 70x when one block holds many unique sequences at a tight threshold.
+
+## NEW FEATURES
+* Added the `expand` argument to `buildNetwork()`. `"clique"` (default) materializes every pairwise edge and reproduces the exact edge multiplicity that community-detection clustering expects. `"star"` links identical sequences through a single hub and connects related groups hub to hub. `"star"` produces far fewer edges and preserves connected components exactly, so it is a large memory win when the downstream step depends on connectivity.
+
+## BUG FIXES
+* Fixed a banded Levenshtein error that dropped pairs whose only optimal alignment rides the band edge, including identical sequences at a tight normalized threshold. The boundary cell read the wrong diagonal value. The fix only recovers true edges and never removes valid ones.
+* Made the normalized-threshold cutoff robust to floating-point rounding so a normalized distance exactly equal to the threshold is kept, matching the documented inclusive behavior.
+
  # immApex VERSION 1.5.4
 
 ## BUG FIXES

R/RcppExports.R

@@ -9,7 +9,7 @@ encodeSequences_cpp <- function(sequences, mode = "onehot", alphabet = as.charac
     .Call(`_immApex_encodeSequences_cpp`, sequences, mode, alphabet, prop_mat_, pad_token, summary, max_len, nthreads)
 }
 
-fast_edge_list <- function(seqs, thresh = 1.0, v_gene = NULL, j_gene = NULL, match_v = FALSE, match_j = FALSE, ids = NULL, metric = "levenshtein", normalize = "none", subst_matrix = NULL, gap_open = -10L, gap_extend = -1L) {
-    .Call(`_immApex_fast_edge_list`, seqs, thresh, v_gene, j_gene, match_v, match_j, ids, metric, normalize, subst_matrix, gap_open, gap_extend)
+fast_edge_list <- function(seqs, thresh = 1.0, v_gene = NULL, j_gene = NULL, match_v = FALSE, match_j = FALSE, metric = "levenshtein", normalize = "none", subst_matrix = NULL, gap_open = -10L, gap_extend = -1L, expand = "clique") {
+    .Call(`_immApex_fast_edge_list`, seqs, thresh, v_gene, j_gene, match_v, match_j, metric, normalize, subst_matrix, gap_open, gap_extend, expand)
 }
 

R/buildNetwork.R

@@ -21,6 +21,14 @@
 #' `data.frame` **or** a symmetric `Matrix::dgCMatrix` adjacency matrix.
 #' @param weight `"dist"` (store the edit distance) **or** `"binary"`
 #' (all edges get weight 1). Ignored when `output = "edges"`.
+#' @param expand How identical sequences (and groups of related sequences) are
+#' expanded into edges. `"clique"` (default) materializes every pairwise edge,
+#' reproducing the exact edge multiplicity expected by community-detection
+#' clustering. `"star"` links each set of identical sequences through a single
+#' hub and connects related groups hub-to-hub. `"star"` produces far fewer edges
+#' (a large memory win on clonally expanded repertoires) and preserves connected
+#' components exactly, but changes edge multiplicity, so use it only when the
+#' downstream step depends on connectivity rather than weighted degree.
 #' @param dist_type Character string specifying the distance metric to use:
 #'  \itemize{
 #'    \item{`"levenshtein"`}  - Standard edit distance (default, backward compatible)
@@ -114,11 +122,13 @@ buildNetwork <- function(input.data        = NULL,
                          filter.j          = FALSE,
                          ids               = NULL,
                          output            = c("edges", "sparse"),
-                         weight            = c("dist", "binary")) {
-
+                         weight            = c("dist", "binary"),
+                         expand            = c("clique", "star")) {
+
   output    <- match.arg(output)
   weight    <- match.arg(weight)
   normalize <- match.arg(normalize) # Validate input
+  expand    <- match.arg(expand)
 
   ## 1. Decide where sequences come from 
   if (is.null(input.data)) {
@@ -192,35 +202,42 @@ buildNetwork <- function(input.data        = NULL,
     numeric_mat <- .fetch.matrix(dist_mat)
   }
 
-  ## 5. Call the C++ engine 
-  edge_df  <- fast_edge_list(
+  # Node labels: explicit `ids`, else row-names / positional (matches the
+  # former C++ default of "1".."n" for the bare-sequence path).
+  if (is.null(ids)) ids <- as.character(seq_len(n))
+
+  ## 5. Call the C++ engine
+  raw  <- fast_edge_list(
     seqs         = seq_vec,
     thresh       = threshold,
     v_gene       = v_vec,
     j_gene       = j_vec,
     match_v      = filter.v,
     match_j      = filter.j,
-    ids          = ids,
     metric       = dist_type,
     subst_matrix = numeric_mat,
     gap_open     = gap_open,
     gap_extend   = gap_extend,
-    normalize    = normalize  
+    normalize    = normalize,
+    expand       = expand
   )
-
-  # Remove Duplicate Edges
-  if (nrow(edge_df) > 0) {
-    edge_key <- paste(
-      pmin(edge_df$from, edge_df$to),
-      pmax(edge_df$from, edge_df$to),
-      sep = "\t"
-    )
-
-    # Check if there are duplicates
-    if (anyDuplicated(edge_key)) {
-      agg_idx <- !duplicated(edge_key)
-      edge_df <- edge_df[agg_idx, , drop = FALSE]
-    }
+
+  # The engine returns 0-based node indices; attach the labels here so the
+  # millions of barcode strings are never materialized inside the parallel core.
+  edge_df <- data.frame(
+    from = ids[raw$i + 1L],
+    to   = ids[raw$k + 1L],
+    dist = raw$dist,
+    stringsAsFactors = FALSE
+  )
+
+  # The engine emits each unordered node pair exactly once, so duplicate edges
+  # can only appear when non-unique `ids` collapse distinct nodes onto the same
+  # label. Skip the (otherwise wasteful) key build unless that is the case.
+  if (nrow(edge_df) > 0 && anyDuplicated(ids)) {
+    edge_key <- paste(pmin(edge_df$from, edge_df$to),
+                      pmax(edge_df$from, edge_df$to), sep = "\t")
+    edge_df <- edge_df[!duplicated(edge_key), , drop = FALSE]
   }
 
   if (output == "edges")

data-raw/make-buildNetwork-golden.R

@@ -0,0 +1,14 @@
+# Regenerate the buildNetwork regression golden snapshot.
+# Run deliberately (NOT as part of the rewrite) only when an intended change to
+# buildNetwork's output set is being accepted. From the package root:
+#   Rscript data-raw/make-buildNetwork-golden.R
+suppressMessages(devtools::load_all(".", quiet = TRUE))
+source("tests/testthat/helper-buildNetwork-regression.R")
+sc <- .bn_scenarios()
+extra <- if ("expand" %in% names(formals(buildNetwork))) list(expand = "clique") else list()
+golden <- lapply(sc$grid, function(scen) {
+  set <- suppressMessages(.bn_run(sc$datasets, scen, extra))
+  list(data = scen$data, args = scen$args, set = set)
+})
+saveRDS(golden, "tests/testthat/fixtures/buildNetwork_golden.rds")
+cat("Wrote", length(golden), "scenarios.\n")

src/RcppExports.cpp

@@ -44,8 +44,8 @@ BEGIN_RCPP
 END_RCPP
 }
 // fast_edge_list
-DataFrame fast_edge_list(CharacterVector seqs, double thresh, Nullable<CharacterVector> v_gene, Nullable<CharacterVector> j_gene, bool match_v, bool match_j, Nullable<CharacterVector> ids, std::string metric, std::string normalize, Nullable<NumericMatrix> subst_matrix, int gap_open, int gap_extend);
-RcppExport SEXP _immApex_fast_edge_list(SEXP seqsSEXP, SEXP threshSEXP, SEXP v_geneSEXP, SEXP j_geneSEXP, SEXP match_vSEXP, SEXP match_jSEXP, SEXP idsSEXP, SEXP metricSEXP, SEXP normalizeSEXP, SEXP subst_matrixSEXP, SEXP gap_openSEXP, SEXP gap_extendSEXP) {
+DataFrame fast_edge_list(CharacterVector seqs, double thresh, Nullable<CharacterVector> v_gene, Nullable<CharacterVector> j_gene, bool match_v, bool match_j, std::string metric, std::string normalize, Nullable<NumericMatrix> subst_matrix, int gap_open, int gap_extend, std::string expand);
+RcppExport SEXP _immApex_fast_edge_list(SEXP seqsSEXP, SEXP threshSEXP, SEXP v_geneSEXP, SEXP j_geneSEXP, SEXP match_vSEXP, SEXP match_jSEXP, SEXP metricSEXP, SEXP normalizeSEXP, SEXP subst_matrixSEXP, SEXP gap_openSEXP, SEXP gap_extendSEXP, SEXP expandSEXP) {
 BEGIN_RCPP
     Rcpp::RObject rcpp_result_gen;
     Rcpp::RNGScope rcpp_rngScope_gen;
@@ -55,13 +55,13 @@ BEGIN_RCPP
     Rcpp::traits::input_parameter< Nullable<CharacterVector> >::type j_gene(j_geneSEXP);
     Rcpp::traits::input_parameter< bool >::type match_v(match_vSEXP);
     Rcpp::traits::input_parameter< bool >::type match_j(match_jSEXP);
-    Rcpp::traits::input_parameter< Nullable<CharacterVector> >::type ids(idsSEXP);
     Rcpp::traits::input_parameter< std::string >::type metric(metricSEXP);
     Rcpp::traits::input_parameter< std::string >::type normalize(normalizeSEXP);
     Rcpp::traits::input_parameter< Nullable<NumericMatrix> >::type subst_matrix(subst_matrixSEXP);
     Rcpp::traits::input_parameter< int >::type gap_open(gap_openSEXP);
     Rcpp::traits::input_parameter< int >::type gap_extend(gap_extendSEXP);
-    rcpp_result_gen = Rcpp::wrap(fast_edge_list(seqs, thresh, v_gene, j_gene, match_v, match_j, ids, metric, normalize, subst_matrix, gap_open, gap_extend));
+    Rcpp::traits::input_parameter< std::string >::type expand(expandSEXP);
+    rcpp_result_gen = Rcpp::wrap(fast_edge_list(seqs, thresh, v_gene, j_gene, match_v, match_j, metric, normalize, subst_matrix, gap_open, gap_extend, expand));
     return rcpp_result_gen;
 END_RCPP
 }