diff --git a/DESCRIPTION b/DESCRIPTION index 91ecd43..28c95bc 100644 --- a/DESCRIPTION +++ b/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 diff --git a/NEWS.md b/NEWS.md index 0f35831..a631ee0 100644 --- a/NEWS.md +++ b/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 diff --git a/R/RcppExports.R b/R/RcppExports.R index 7188a6c..77d0d71 100644 --- a/R/RcppExports.R +++ b/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) } diff --git a/R/buildNetwork.R b/R/buildNetwork.R index 09d6bdb..69aba55 100644 --- a/R/buildNetwork.R +++ b/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") diff --git a/data-raw/make-buildNetwork-golden.R b/data-raw/make-buildNetwork-golden.R new file mode 100644 index 0000000..89ea189 --- /dev/null +++ b/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") diff --git a/man/amino.acids.Rd b/man/amino.acids.Rd index 8ccde83..05f04ec 100644 --- a/man/amino.acids.Rd +++ b/man/amino.acids.Rd @@ -1,16 +1,11 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/utils.R -\docType{data} \name{amino.acids} \alias{amino.acids} \title{Standard 20 amino acids} -\format{ -An object of class \code{character} of length 20. -} \usage{ amino.acids } \description{ Vector of one-letter codes for the 20 standard amino acids. } -\keyword{datasets} diff --git a/man/buildNetwork.Rd b/man/buildNetwork.Rd index 1e045cd..3b42ca0 100644 --- a/man/buildNetwork.Rd +++ b/man/buildNetwork.Rd @@ -20,7 +20,8 @@ buildNetwork( filter.j = FALSE, ids = NULL, output = c("edges", "sparse"), - weight = c("dist", "binary") + weight = c("dist", "binary"), + expand = c("clique", "star") ) } \arguments{ @@ -86,6 +87,15 @@ Only used when `metric` is "nw" or "sw".} \item{weight}{`"dist"` (store the edit distance) **or** `"binary"` (all edges get weight 1). Ignored when `output = "edges"`.} + +\item{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.} } \value{ edge-list `data.frame` **or** sparse adjacency `dgCMatrix` of distances diff --git a/man/immApex-package.Rd b/man/immApex-package.Rd index 69248e2..83fd26b 100644 --- a/man/immApex-package.Rd +++ b/man/immApex-package.Rd @@ -19,6 +19,11 @@ Useful links: \author{ \strong{Maintainer}: Nick Borcherding \email{ncborch@gmail.com} [copyright holder] +Authors: +\itemize{ + \item Nick Borcherding \email{ncborch@gmail.com} [copyright holder] +} + Other contributors: \itemize{ \item Qile Yang \email{qile.yang@berkeley.edu} (\href{https://orcid.org/0009-0005-0148-2499}{ORCID}) [contributor] diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 5d2318d..d0412e3 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -44,8 +44,8 @@ BEGIN_RCPP END_RCPP } // fast_edge_list -DataFrame fast_edge_list(CharacterVector seqs, double thresh, Nullable v_gene, Nullable j_gene, bool match_v, bool match_j, Nullable ids, std::string metric, std::string normalize, Nullable 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 v_gene, Nullable j_gene, bool match_v, bool match_j, std::string metric, std::string normalize, Nullable 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 >::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 >::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 >::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 } diff --git a/src/fastEditEdges.cpp b/src/fastEditEdges.cpp index c856b25..d4d450c 100644 --- a/src/fastEditEdges.cpp +++ b/src/fastEditEdges.cpp @@ -4,6 +4,8 @@ #include #include #include +#include +#include #include #ifdef _OPENMP #include @@ -11,6 +13,49 @@ using namespace Rcpp; +// Lift floating-point products back over an integer they should have reached: +// thresh * norm_len can land at e.g. 2.9999999996 when the true value is 3, and +// a bare (int) cast would drop a pair whose normalized distance exactly equals +// the threshold (the inclusive `fd <= thresh` boundary). The epsilon dwarfs the +// FP rounding error (~norm_len * 2e-16) while never reaching the next integer. +static const double THRESH_EPS = 1e-9; + +// ============================================================================ +// Deletion-neighborhood (SymSpell) helpers +// ============================================================================ +// All strings reachable by deleting up to `k` characters. For Levenshtein, +// edit(a,b) <= k implies their deletion neighborhoods intersect, so colliding +// reps are exact candidates (no false negatives) verified later by banded DP. +// A variant of length L-d is always reached with exactly d deletions, so the +// "insert-then-recurse only on first sight" prune is loss-free. +static void gen_deletes_rec(const std::string& cur, int rem, + std::unordered_set& out) { + if (rem == 0) return; + for (size_t i = 0; i < cur.size(); ++i) { + std::string nxt; + nxt.reserve(cur.size() - 1); + nxt.append(cur, 0, i); + nxt.append(cur, i + 1, std::string::npos); + if (out.insert(nxt).second) gen_deletes_rec(nxt, rem - 1, out); + } +} +static void gen_deletes(const std::string& s, int k, + std::unordered_set& out) { + out.clear(); + out.insert(s); + if (k > 0) gen_deletes_rec(s, k, out); +} +// Rough upper bound on variants per sequence: sum_{j=0}^{k} C(L, j), capped. +static long deletion_variant_bound(int L, int k) { + long tot = 1, c = 1; + for (int j = 1; j <= k; ++j) { + c = c * (L - j + 1) / j; + tot += c; + if (tot > 1000000L) return tot; + } + return tot; +} + // ============================================================================ // Data Structures & Buffers (For Speed) // ============================================================================ @@ -58,13 +103,18 @@ static inline int levenshtein_opt(const std::string &a, const std::string &b, in std::iota(row.begin(), row.begin() + m + 1, 0); for (int i = 1; i <= n; ++i) { - int prev = row[0]; - row[0] = i; - int row_min = i; - int start_j = std::max(1, i - thr); int end_j = std::min(m, i + thr); - + + // prev must hold the diagonal D[i-1][start_j-1]. When the band starts past + // column 1 this is row[start_j-1] BEFORE it is overwritten by the boundary + // guard - seeding it with D[i-1][0] (the old behaviour) corrupts the first + // banded cell and makes the routine miss pairs whose only optimal alignment + // rides the band edge (e.g. distance == |len diff| == thr). + int prev = row[start_j - 1]; + row[0] = i; + int row_min = i; + // Fill valid band if (start_j > 1) row[start_j - 1] = thr + 2; // Guard boundary @@ -204,13 +254,14 @@ DataFrame fast_edge_list(CharacterVector seqs, Nullable j_gene = R_NilValue, bool match_v = false, bool match_j = false, - Nullable ids = R_NilValue, std::string metric = "levenshtein", std::string normalize = "none", Nullable subst_matrix = R_NilValue, int gap_open = -10, - int gap_extend = -1) + int gap_extend = -1, + std::string expand = "clique") { + const bool star = (expand == "star"); int n = seqs.size(); if (n < 2) stop("Need at least 2 sequences."); @@ -223,16 +274,42 @@ DataFrame fast_edge_list(CharacterVector seqs, if(match_v && v_gene.isNotNull()) v_vec = as>(v_gene); if(match_j && j_gene.isNotNull()) j_vec = as>(j_gene); - std::vector lbl(n); - if(ids.isNotNull()) lbl = as>(ids); - else for(int i=0; i, std::vector> group_map; + // 2. Deduplicate to representatives. + // A representative is a unique (sequence, V-used, J-used) tuple. V/J are + // only part of the key when the corresponding filter is on, so a rep never + // spans more than one (V,J) block. All distance work runs over reps; the + // member node indices are expanded back into edges afterwards. + const char SEP = '\x1f'; // unit separator: cannot occur in seq/gene + std::unordered_map rep_of; + std::vector> rep_members; // rep id -> node indices (ascending) + std::vector rep_seq; + std::vector rep_len; + std::map> block_reps; // block key -> rep ids + + rep_of.reserve(n * 2); for(int i=0; i{i}); + block_reps[block_key].push_back(r); + } else { + rep_members[it->second].push_back(i); + } } - + // 3. Prepare Substitution Matrix SubstMatrix smat; if (metric == "nw" || metric == "sw") { @@ -254,96 +331,287 @@ DataFrame fast_edge_list(CharacterVector seqs, } } - // 4. Output Storage - std::vector out_from, out_to; + // 4. Output Storage (0-based node indices; labels are attached in R) + std::vector out_i, out_k; std::vector out_d; - // 5. Parallel Processing -#ifdef _OPENMP - int nthreads = omp_get_max_threads(); -#else - int nthreads = 1; -#endif - - std::vector> all_groups; - for(auto &kv : group_map) if(kv.second.size() > 1) all_groups.push_back(kv.second); - - for (const auto& grp : all_groups) { - int g_size = grp.size(); - -#pragma omp parallel num_threads(nthreads) + // 5. Per-pair decision (the exact, brute-force-faithful verifier). + // Encodes metric / normalize as ints so the hot path avoids string + // comparisons. Returns true and sets `fd_out` when an edge is kept, + // reproducing the original maxd / length-gate / normalize / post-gate + // logic byte-for-byte. + enum Metric { LEV, HAM, DAM, MNW, MSW }; + enum Norm { N_NONE, N_MAXLEN, N_MEANLEN }; + Metric mcode = (metric == "hamming") ? HAM : + (metric == "damerau") ? DAM : + (metric == "nw") ? MNW : + (metric == "sw") ? MSW : LEV; + Norm ncode = (normalize == "maxlen") ? N_MAXLEN : + (normalize == "length") ? N_MEANLEN : N_NONE; + + auto compute_pair = [&](const std::string& A, int la, + const std::string& B, int lb, + DPWorkspace& ws, double& fd_out) -> bool { + int max_l = la > lb ? la : lb; + if (max_l == 0) return false; // matches original `max_l==0 continue` + int maxd; + if (thresh >= 1.0) { + maxd = (int)thresh; + } else { + double norm_len = (ncode == N_MEANLEN) ? (la + lb) / 2.0 : (double)max_l; + maxd = (int)(thresh * norm_len + THRESH_EPS); + } + if (mcode != MNW && mcode != MSW && std::abs(la - lb) > maxd) return false; + + int d = 0; + switch (mcode) { + case LEV: d = levenshtein_opt(A, B, maxd, ws.v1); break; + case HAM: + if (la != lb) d = maxd + 1; + else { for (int p = 0; p < la; ++p) if (A[p] != B[p]) { if (++d > maxd) break; } } + break; + case DAM: d = damerau_opt(A, B, maxd, ws); break; + case MNW: d = nw_opt(A, B, smat, ws); break; + case MSW: d = sw_opt(A, B, smat, ws); break; + } + if (d > maxd) return false; + + double fd = (double)d; + if (ncode == N_MAXLEN) fd /= max_l; + else if (ncode == N_MEANLEN) fd /= ((la + lb) / 2.0); + if (thresh < 1.0 && ncode != N_NONE && fd > thresh) return false; + + fd_out = fd; + return true; + }; + + // 6. Per-block index: reps sorted by id, length buckets, and a block-level + // integer distance cap k_block. Because every per-pair maxd <= k_block + // (norm_len <= Lmax_block), restricting candidates to length-compatible + // buckets never drops a real edge - the exact per-pair verifier decides. + int R = (int)rep_seq.size(); + struct BlockIdx { + std::vector rep_ids; // ascending rep id + std::map> buckets; // length -> ascending rep ids + int kblock = 0; + int Lmax = 0; + }; + std::vector bidx; + bidx.reserve(block_reps.size()); + std::vector rep_block(R); + { + int b = 0; + for (auto& kv : block_reps) { + BlockIdx bi; + bi.rep_ids = std::move(kv.second); + int Lmax = 0; + for (int r : bi.rep_ids) { + rep_block[r] = b; + bi.buckets[rep_len[r]].push_back(r); + if (rep_len[r] > Lmax) Lmax = rep_len[r]; + } + bi.kblock = (thresh >= 1.0) ? (int)thresh : (int)(thresh * (double)Lmax + THRESH_EPS); + bi.Lmax = Lmax; + bidx.push_back(std::move(bi)); + ++b; + } + } + + // 6b. Hamming pigeonhole index (one per equal-length bucket). Two equal-length + // sequences within Hamming distance k must share at least one of k+1 + // positioned segments, so candidate partners are the reps that collide on + // a segment. Buckets where k >= L (every pair trivially passes) fall back + // to a plain bucket scan. + struct SegBucket { + bool brute = false; + int nseg = 0; + std::vector seg_start, seg_len; + std::unordered_map> index; // pos-tag+seg -> ascending rep ids + }; + std::vector> ham(mcode == HAM ? bidx.size() : 0); + if (mcode == HAM) { + for (size_t b = 0; b < bidx.size(); ++b) { + for (auto& lb : bidx[b].buckets) { + int L = lb.first; + const std::vector& reps = lb.second; + if (L == 0 || reps.size() < 2) continue; // no possible edges + int kseg = (thresh >= 1.0) ? (int)thresh : (int)(thresh * (double)L + THRESH_EPS); + SegBucket sb; + if (kseg >= L) { sb.brute = true; ham[b].emplace(L, std::move(sb)); continue; } + sb.nseg = kseg + 1; + int base = L / sb.nseg, rem = L % sb.nseg, pos = 0; + for (int j = 0; j < sb.nseg; ++j) { + int len = base + (j < rem ? 1 : 0); + sb.seg_start.push_back(pos); + sb.seg_len.push_back(len); + pos += len; + } + sb.index.reserve(reps.size() * sb.nseg); + for (int r : reps) { + const std::string& sq = rep_seq[r]; + for (int j = 0; j < sb.nseg; ++j) { + std::string key; + key.reserve(sb.seg_len[j] + 1); + key += (char)(j + 1); // position tag (letters are >= 'A', no clash) + key.append(sq, sb.seg_start[j], sb.seg_len[j]); + sb.index[key].push_back(r); // reps ascending -> postings ascending + } + } + ham[b].emplace(L, std::move(sb)); + } + } + } + + // 6c. Deletion-neighborhood (SymSpell) index per block for edit-distance + // metrics with a small radius. Levenshtein uses radius = k_block; the + // guard keeps the variant explosion bounded, otherwise a block falls back + // to the length-blocked bucket scan. (Damerau, radius 2*k_block, is added + // alongside in step 8.) + const int K_CAP = 3; + const long MAXVAR = 4096; + struct DelIdx { + bool use = false; + int radius = 0; + std::unordered_map> index; // variant -> ascending rep ids + }; + std::vector del(bidx.size()); + if (mcode == LEV || mcode == DAM) { + std::unordered_set vars; + for (size_t b = 0; b < bidx.size(); ++b) { + int k = bidx[b].kblock; + if (k < 0 || k > K_CAP) continue; + int radius = (mcode == DAM) ? 2 * k : k; // transposition costs 2 in Levenshtein + if (deletion_variant_bound(bidx[b].Lmax, radius) > MAXVAR) continue; + DelIdx& di = del[b]; + di.use = true; + di.radius = radius; + di.index.reserve(bidx[b].rep_ids.size() * 4); + for (int r : bidx[b].rep_ids) { + gen_deletes(rep_seq[r], radius, vars); + for (const std::string& v : vars) di.index[v].push_back(r); + } + } + } + + // 7. Parallel over reps. Each rep emits its identical-member edges (clique or + // star) plus the edges to length-compatible reps with a higher id (so each + // unordered pair is generated exactly once). Expansion maps rep-pairs back + // to member node indices. +#pragma omp parallel { - // Thread-Local Storage - DPWorkspace ws; - std::vector loc_f, loc_t; - std::vector loc_d_vec; - -#pragma omp for schedule(dynamic, 64) - for (int i_idx = 0; i_idx < g_size; ++i_idx) { - int i = grp[i_idx]; - for (int k_idx = i_idx + 1; k_idx < g_size; ++k_idx) { - int k = grp[k_idx]; - - int max_l = std::max(lens[i], lens[k]); - if (max_l == 0) continue; - - int maxd; - if (thresh >= 1.0) { - // Absolute distance threshold - maxd = (int)thresh; - } else { - // Normalized distance threshold - use same length as normalization - double norm_len; - if (normalize == "maxlen") { - norm_len = (double)max_l; - } else if (normalize == "length") { - norm_len = (lens[i] + lens[k]) / 2.0; - } else { - // normalize == "none": interpret threshold as fraction of max_l - norm_len = (double)max_l; + DPWorkspace ws; + std::vector loc_f, loc_t; + std::vector loc_d; + std::vector cand; // reusable candidate-partner buffer + std::unordered_set vars; // reusable deletion-variant buffer + + // verify a rep-pair and, if kept, expand to member node pairs + auto verify_expand = [&](int r, int sv) { + double fd; + if (!compute_pair(rep_seq[r], rep_len[r], rep_seq[sv], rep_len[sv], ws, fd)) return; + const std::vector& mr = rep_members[r]; + const std::vector& ms = rep_members[sv]; + if (star) { + loc_f.push_back(mr[0]); loc_t.push_back(ms[0]); loc_d.push_back(fd); + } else { + for (int a : mr) for (int b : ms) { + loc_f.push_back(a); loc_t.push_back(b); loc_d.push_back(fd); + } + } + }; + // consider partners in an ascending bucket whose rep id exceeds r + auto consider = [&](int r, const std::vector& bucket) { + for (auto it = std::upper_bound(bucket.begin(), bucket.end(), r); + it != bucket.end(); ++it) + verify_expand(r, *it); + }; + +#pragma omp for schedule(dynamic, 32) + for (int r = 0; r < R; ++r) { + // intra-rep edges (identical sequences, mutual distance d_self). + const std::vector& mr = rep_members[r]; + double fd; + if (mr.size() >= 2 && + compute_pair(rep_seq[r], rep_len[r], rep_seq[r], rep_len[r], ws, fd)) { + if (star) { + int hub = mr[0]; + for (size_t a = 1; a < mr.size(); ++a) { + loc_f.push_back(hub); loc_t.push_back(mr[a]); loc_d.push_back(fd); } - maxd = (int)(thresh * norm_len); + } else { + for (size_t a = 0; a < mr.size(); ++a) + for (size_t b = a + 1; b < mr.size(); ++b) { + loc_f.push_back(mr[a]); loc_t.push_back(mr[b]); loc_d.push_back(fd); + } } - // ===================================================================== - - if (metric != "nw" && metric != "sw" && std::abs(lens[i] - lens[k]) > maxd) continue; - - int d = 0; - if (metric == "levenshtein") d = levenshtein_opt(s[i], s[k], maxd, ws.v1); - else if (metric == "hamming") { - if (lens[i] != lens[k]) d = maxd + 1; - else { - for(int p=0; pmaxd) break; } + } + + // inter-rep candidates, restricted by the (V,J,length) index. + const BlockIdx& bi = bidx[rep_block[r]]; + if (mcode == MNW || mcode == MSW) { + consider(r, bi.rep_ids); // alignment scores: no length pruning + } else if (mcode == HAM) { + // Hamming: equal length only, candidates via pigeonhole segments. + const std::unordered_map& hmap = ham[rep_block[r]]; + auto hit = hmap.find(rep_len[r]); + if (hit == hmap.end()) { + // bucket had < 2 reps or empty strings -> no edges + } else if (hit->second.brute) { + auto bit = bi.buckets.find(rep_len[r]); + if (bit != bi.buckets.end()) consider(r, bit->second); + } else { + const SegBucket& sb = hit->second; + const std::string& sq = rep_seq[r]; + cand.clear(); + for (int j = 0; j < sb.nseg; ++j) { + std::string key; + key.reserve(sb.seg_len[j] + 1); + key += (char)(j + 1); + key.append(sq, sb.seg_start[j], sb.seg_len[j]); + auto pit = sb.index.find(key); + if (pit != sb.index.end()) + for (auto p = std::upper_bound(pit->second.begin(), pit->second.end(), r); + p != pit->second.end(); ++p) + cand.push_back(*p); } + std::sort(cand.begin(), cand.end()); + cand.erase(std::unique(cand.begin(), cand.end()), cand.end()); + for (int s : cand) verify_expand(r, s); } - else if (metric == "damerau") d = damerau_opt(s[i], s[k], maxd, ws); - else if (metric == "nw") d = nw_opt(s[i], s[k], smat, ws); - else if (metric == "sw") d = sw_opt(s[i], s[k], smat, ws); - - if (d <= maxd) { - double fd = (double)d; - if (normalize == "maxlen") fd /= max_l; - else if (normalize == "length") fd /= ((lens[i]+lens[k])/2.0); - - if (thresh < 1.0 && normalize != "none" && fd > thresh) { - continue; // Skip edges that exceed normalized threshold + } else { + // edit distance (Levenshtein / Damerau) + const DelIdx& di = del[rep_block[r]]; + if (di.use) { + // SymSpell: candidates are reps sharing a deletion variant + gen_deletes(rep_seq[r], di.radius, vars); + cand.clear(); + for (const std::string& v : vars) { + auto it = di.index.find(v); + if (it != di.index.end()) + for (auto p = std::upper_bound(it->second.begin(), it->second.end(), r); + p != it->second.end(); ++p) + cand.push_back(*p); + } + std::sort(cand.begin(), cand.end()); + cand.erase(std::unique(cand.begin(), cand.end()), cand.end()); + for (int s : cand) verify_expand(r, s); + } else { + // fallback: length-blocked bucket scan, |dL| <= k_block + int L = rep_len[r], k = bi.kblock; + for (int Lp = L - k; Lp <= L + k; ++Lp) { + auto it = bi.buckets.find(Lp); + if (it != bi.buckets.end()) consider(r, it->second); } - // ===================================================================== - - loc_f.push_back(lbl[i]); - loc_t.push_back(lbl[k]); - loc_d_vec.push_back(fd); } } } #pragma omp critical { - out_from.insert(out_from.end(), loc_f.begin(), loc_f.end()); - out_to.insert(out_to.end(), loc_t.begin(), loc_t.end()); - out_d.insert(out_d.end(), loc_d_vec.begin(), loc_d_vec.end()); + out_i.insert(out_i.end(), loc_f.begin(), loc_f.end()); + out_k.insert(out_k.end(), loc_t.begin(), loc_t.end()); + out_d.insert(out_d.end(), loc_d.begin(), loc_d.end()); } } - } - - return DataFrame::create(_["from"] = out_from, _["to"] = out_to, _["dist"] = out_d, _["stringsAsFactors"] = false); + + return DataFrame::create(_["i"] = out_i, _["k"] = out_k, _["dist"] = out_d, _["stringsAsFactors"] = false); } diff --git a/tests/testthat/fixtures/buildNetwork_golden.rds b/tests/testthat/fixtures/buildNetwork_golden.rds new file mode 100644 index 0000000..90f0f38 Binary files /dev/null and b/tests/testthat/fixtures/buildNetwork_golden.rds differ diff --git a/tests/testthat/helper-buildNetwork-regression.R b/tests/testthat/helper-buildNetwork-regression.R new file mode 100644 index 0000000..a962d1a --- /dev/null +++ b/tests/testthat/helper-buildNetwork-regression.R @@ -0,0 +1,193 @@ +# --------------------------------------------------------------------------- +# Regression harness for buildNetwork() +# +# Goal: freeze the *set* of edges (and their distances) that the current +# implementation produces, so that the candidate-indexing rewrite can be +# proven bit-identical in `expand = "clique"` mode across the full parameter +# matrix. Row order is intentionally ignored (OpenMP merge order is not +# deterministic); we compare canonical sets. +# +# Both the golden-snapshot generator (data-raw style script) and the +# regression test source these helpers, so the scenario grid lives in one +# place. +# --------------------------------------------------------------------------- + +# Canonical, order-independent representation of a buildNetwork() result. +# Handles both the edge-list data.frame and the sparse dgCMatrix output. +.bn_canon <- function(res) { + if (inherits(res, "data.frame")) { + if (nrow(res) == 0L) return(character(0)) + a <- pmin(res$from, res$to) + b <- pmax(res$from, res$to) + # %.12g round-trips the integer and normalized doubles exactly while + # staying robust to trivial formatting noise. method = "radix" sorts in + # byte order, so the canonical form does not depend on LC_COLLATE (the + # golden was generated under one locale and R CMD check runs under C). + key <- paste(a, b, sprintf("%.12g", res$dist), sep = "\031") + sort(key, method = "radix") + } else if (inherits(res, "dgCMatrix")) { + res <- methods::as(res, "TsparseMatrix") + nm <- rownames(res) + fr <- nm[res@i + 1L] + to <- nm[res@j + 1L] + a <- pmin(fr, to) + b <- pmax(fr, to) + key <- paste(a, b, sprintf("%.12g", res@x), sep = "\031") + sort(unique(key), method = "radix") + } else { + stop("Unexpected buildNetwork result class: ", class(res)[1]) + } +} + +# Deterministic data generators (fixed seeds) covering the tricky cases. +.bn_data <- function() { + aa <- c("A","C","D","E","F","G","H","I","K","L", + "M","N","P","Q","R","S","T","V","W","Y") + rseq <- function(len) paste(sample(aa, len, replace = TRUE), collapse = "") + + out <- list() + + # 1. Duplicate-heavy: many identical (seq,V,J) tuples across distinct ids, + # plus near-neighbors. This is the dedup/expansion stress case. + set.seed(101) + base <- replicate(15, rseq(sample(12:16, 1))) + seqs <- sample(base, 120, replace = TRUE) # heavy repetition + # inject a handful of 1-2 edit neighbors of base[1] + mut <- base[1] + near <- vapply(1:6, function(k) { + s <- strsplit(mut, "")[[1]] + p <- sample(seq_along(s), 1); s[p] <- sample(aa, 1) + paste(s, collapse = "") + }, character(1)) + seqs <- c(seqs, near) + n <- length(seqs) + out$dup_heavy <- data.frame( + sequence_id = paste0("dup_", seq_len(n)), + junction_aa = seqs, + v_call = paste0("IGHV", sample(1:3, n, replace = TRUE)), + j_call = paste0("IGHJ", sample(1:2, n, replace = TRUE)), + stringsAsFactors = FALSE) + + # 2. Mostly unique sequences, mixed lengths (length-blocking stress). + set.seed(202) + m <- 90 + out$unique_mixed <- data.frame( + sequence_id = paste0("u_", seq_len(m)), + junction_aa = vapply(seq_len(m), function(i) rseq(sample(10:20, 1)), character(1)), + v_call = paste0("IGHV", sample(1:4, m, replace = TRUE)), + j_call = paste0("IGHJ", sample(1:3, m, replace = TRUE)), + stringsAsFactors = FALSE) + + # 3. Empty strings + single-member groups + short sequences. + out$edge_cases <- data.frame( + sequence_id = paste0("e_", 1:12), + junction_aa = c("", "", "AC", "ACD", "ACDE", "ACDF", + "ACDEFGHIK", "ACDEFGHIL", "ACDEFGHIK", + "LONELYSEQ1", "WXYV", "WXYW"), + v_call = c("V1","V1","V1","V1","V1","V1", + "V2","V2","V2","V9","V3","V3"), + j_call = c("J1","J1","J1","J1","J1","J1", + "J1","J1","J1","J9","J1","J1"), + stringsAsFactors = FALSE) + + # 4. Equal-length block for Hamming (with duplicates + 1/2 substitutions). + set.seed(303) + L <- 14L; q <- 60 + eq <- replicate(8, rseq(L)) + hs <- sample(eq, q, replace = TRUE) + out$equal_len <- data.frame( + sequence_id = paste0("h_", seq_len(q)), + junction_aa = hs, + v_call = paste0("IGHV", sample(1:2, q, replace = TRUE)), + j_call = "IGHJ1", + stringsAsFactors = FALSE) + + # 5. Pure transposition pairs (Damerau radius-2k trap). seq vs adjacent swap. + out$transpose <- data.frame( + sequence_id = paste0("t_", 1:8), + junction_aa = c("ABCDEFGHIK", "BACDEFGHIK", # swap pos 1-2 (Dam 1, Lev 2) + "ABCDEFGHIK", "ABDCEFGHIK", # swap pos 3-4 + "MNPQRSTVWY", "MNPQRSTVYW", # swap pos 9-10 + "MNPQRSTVWY", "NMPQRSTVWY"), # swap pos 1-2 + v_call = "IGHV1", + j_call = "IGHJ1", + stringsAsFactors = FALSE) + + out +} + +# The parameter grid. Each scenario is list(data=, args=). +# `args` are passed to buildNetwork() along with seq_col/ids. +.bn_scenarios <- function() { + datasets <- .bn_data() + grid <- list() + add <- function(data_name, ...) { + grid[[length(grid) + 1L]] <<- list(data = data_name, args = list(...)) + } + + filters <- list(c(FALSE, FALSE), c(TRUE, FALSE), + c(FALSE, TRUE), c(TRUE, TRUE)) + + # Levenshtein / Damerau: absolute + normalized thresholds, all filters. + for (dn in c("dup_heavy", "unique_mixed", "edge_cases")) { + for (mt in c("levenshtein", "damerau")) { + for (f in filters) { + add(dn, dist_type = mt, threshold = 1, filter.v = f[1], filter.j = f[2]) + add(dn, dist_type = mt, threshold = 2, filter.v = f[1], filter.j = f[2]) + add(dn, dist_type = mt, threshold = 3, filter.v = f[1], filter.j = f[2]) + add(dn, dist_type = mt, threshold = 0.1, normalize = "maxlen", + filter.v = f[1], filter.j = f[2]) + add(dn, dist_type = mt, threshold = 0.15, normalize = "length", + filter.v = f[1], filter.j = f[2]) + # the loose-then-refilter pattern scRepertoire uses today + add(dn, dist_type = mt, threshold = 0.85, normalize = "length", + filter.v = f[1], filter.j = f[2]) + } + } + } + + # Hamming: equal-length and edge-case data, absolute + normalized. + for (dn in c("equal_len", "edge_cases", "dup_heavy")) { + for (f in filters) { + add(dn, dist_type = "hamming", threshold = 1, filter.v = f[1], filter.j = f[2]) + add(dn, dist_type = "hamming", threshold = 2, filter.v = f[1], filter.j = f[2]) + add(dn, dist_type = "hamming", threshold = 3, filter.v = f[1], filter.j = f[2]) + add(dn, dist_type = "hamming", threshold = 0.15, normalize = "maxlen", + filter.v = f[1], filter.j = f[2]) + } + } + + # Damerau transposition stress. + add("transpose", dist_type = "damerau", threshold = 1) + add("transpose", dist_type = "damerau", threshold = 2) + add("transpose", dist_type = "levenshtein", threshold = 1) + add("transpose", dist_type = "levenshtein", threshold = 2) + + # NW / SW (no pruning path): a few thresholds. + for (dn in c("unique_mixed", "equal_len")) { + add(dn, dist_type = "nw", dist_mat = "BLOSUM62", threshold = 5, filter.v = TRUE) + add(dn, dist_type = "sw", dist_mat = "BLOSUM62", threshold = 5, filter.v = TRUE) + add(dn, dist_type = "nw", dist_mat = "PAM30", threshold = 8, filter.v = TRUE) + } + + # A couple of sparse-output scenarios to lock the matrix path. + add("dup_heavy", dist_type = "levenshtein", threshold = 2, + filter.v = TRUE, output = "sparse", weight = "dist") + add("dup_heavy", dist_type = "levenshtein", threshold = 2, + filter.v = TRUE, output = "sparse", weight = "binary") + + list(datasets = datasets, grid = grid) +} + +# Run one scenario through buildNetwork() and return its canonical set. +.bn_run <- function(datasets, scen, extra = list()) { + d <- datasets[[scen$data]] + args <- modifyList(scen$args, extra) + args$input.data <- d + args$seq_col <- "junction_aa" + args$v_col <- "v_call" + args$j_col <- "j_call" + args$ids <- d$sequence_id + res <- do.call(buildNetwork, args) + .bn_canon(res) +} diff --git a/tests/testthat/test-buildNetwork-damerau.R b/tests/testthat/test-buildNetwork-damerau.R new file mode 100644 index 0000000..2ac9e5c --- /dev/null +++ b/tests/testthat/test-buildNetwork-damerau.R @@ -0,0 +1,21 @@ +# Guards the SymSpell radius trap: a pure adjacent transposition is Damerau +# distance 1 but Levenshtein distance 2, so the deletion index for Damerau must +# use radius 2*k. A radius-k index would silently miss these pairs. + +test_that("Damerau finds adjacent transpositions that Levenshtein does not", { + d <- data.frame( + sequence_id = c("a", "b", "c", "d"), + junction_aa = c("ABCDEFGHIK", "BACDEFGHIK", # swap 1-2 + "MNPQRSTVWY", "MNPQRSTVYW"), # swap 9-10 + v_call = "V1", j_call = "J1", stringsAsFactors = FALSE) + + e_dam <- suppressMessages(buildNetwork(d, seq_col = "junction_aa", ids = d$sequence_id, + dist_type = "damerau", threshold = 1)) + e_lev <- suppressMessages(buildNetwork(d, seq_col = "junction_aa", ids = d$sequence_id, + dist_type = "levenshtein", threshold = 1)) + + key <- function(df) if (nrow(df)) sort(paste(pmin(df$from, df$to), pmax(df$from, df$to))) else character(0) + expect_identical(key(e_dam), c("a b", "c d")) + expect_true(all(e_dam$dist == 1)) + expect_identical(key(e_lev), character(0)) # swaps are Levenshtein distance 2 +}) diff --git a/tests/testthat/test-buildNetwork-levenshtein-exact.R b/tests/testthat/test-buildNetwork-levenshtein-exact.R new file mode 100644 index 0000000..f811473 --- /dev/null +++ b/tests/testthat/test-buildNetwork-levenshtein-exact.R @@ -0,0 +1,46 @@ +# Direct correctness check for the banded Levenshtein engine, independent of the +# golden snapshot. The banded routine previously missed pairs whose only optimal +# alignment rides the band edge (true distance == |length difference| == thr), +# including identical sequences at a tight normalized threshold. This compares +# buildNetwork's edge set against brute-force utils::adist for several +# thresholds and a length-diverse sequence set (which guarantees band edges are +# exercised). + +test_that("Levenshtein edges match brute-force adist (incl. band edges)", { + set.seed(2024) + aa <- c("A","C","D","E","F","G","H","I","K","L","M","N","P","Q","R","S","T","V","W","Y") + # wide length spread so |len diff| frequently equals the threshold + seqs <- vapply(1:120, function(i) paste(sample(aa, sample(8:22, 1), replace = TRUE), collapse = ""), + character(1)) + ids <- paste0("s", seq_along(seqs)) + D <- as.matrix(utils::adist(seqs)) + + for (thr in 0:4) { + # brute-force truth: all i maxd = floor(0.05 * 9) = 0 + e <- buildNetwork(d, seq_col = "junction_aa", ids = d$sequence_id, + threshold = 0.05, normalize = "maxlen", filter.v = TRUE) + expect_equal(nrow(e), 3L) # all three identical pairs + expect_true(all(e$dist == 0)) +}) diff --git a/tests/testthat/test-buildNetwork-regression.R b/tests/testthat/test-buildNetwork-regression.R new file mode 100644 index 0000000..2015467 --- /dev/null +++ b/tests/testthat/test-buildNetwork-regression.R @@ -0,0 +1,39 @@ +# Regression: the candidate-indexing rewrite must reproduce, as a SET, the +# exact edges + distances of the original brute-force implementation across the +# full parameter matrix. `expand = "clique"` is the bit-identical mode (star +# changes edge multiplicity by design and is checked separately). +# +# The golden snapshot in fixtures/buildNetwork_golden.rds was generated from the +# pre-rewrite implementation. Regenerate ONLY with intent via +# data-raw/make-buildNetwork-golden.R (never blindly). + +test_that("buildNetwork reproduces frozen golden edge sets (clique mode)", { + golden_path <- test_path("fixtures", "buildNetwork_golden.rds") + skip_if_not(file.exists(golden_path), "golden snapshot missing") + golden <- readRDS(golden_path) + sc <- .bn_scenarios() + + # Pass expand = "clique" only once the argument exists (post-rewrite). + extra <- if ("expand" %in% names(formals(buildNetwork))) list(expand = "clique") else list() + + fails <- character(0) + for (i in seq_along(golden)) { + g <- golden[[i]] + scen <- sc$grid[[i]] + got <- tryCatch( + suppressMessages(.bn_run(sc$datasets, scen, extra)), + error = function(e) paste0("ERROR:", conditionMessage(e))) + # Compare as SETS. The contract is the edge set, not row order, so a + # difference in collation order (e.g. C locale under R CMD check vs the + # locale the golden was built in) must not register as a failure. + only_new <- setdiff(got, g$set) + only_old <- setdiff(g$set, got) + if (length(only_new) || length(only_old)) { + lab <- sprintf("[%d] %s {%s}", i, g$data, + paste(names(scen$args), unlist(scen$args), sep = "=", collapse = ", ")) + fails <- c(fails, sprintf("%s +%d/-%d", lab, length(only_new), length(only_old))) + } + } + expect_identical(fails, character(0), + info = paste0("\n", paste(utils::head(fails, 20), collapse = "\n"))) +}) diff --git a/tests/testthat/test-buildNetwork-star.R b/tests/testthat/test-buildNetwork-star.R new file mode 100644 index 0000000..b458177 --- /dev/null +++ b/tests/testthat/test-buildNetwork-star.R @@ -0,0 +1,52 @@ +# `expand = "star"` is NOT bit-identical to clique (it changes edge +# multiplicity by design). Its contract is: identical connected-components +# membership and identical vertex set. We verify that against clique mode on the +# duplicate-heavy fixtures, which is exactly what scRepertoire's default +# `cluster.method = "components"` relies on. + +test_that("star expansion preserves connected components and vertex set", { + skip_if_not_installed("igraph") + sc <- .bn_scenarios() + datasets <- sc$datasets + + # Each connected component as a sorted vertex-set signature; the set of these + # signatures is the partition, compared order-independently. + component_sets <- function(edge_df) { + if (nrow(edge_df) == 0) return(character(0)) + g <- igraph::graph_from_data_frame(edge_df[, c("from", "to")], directed = FALSE) + memb <- igraph::membership(igraph::components(g, mode = "weak")) + groups <- split(names(memb), as.integer(memb)) + sort(vapply(groups, function(v) paste(sort(v), collapse = "\031"), character(1)), + method = "radix") + } + + cfgs <- list( + list(data = "dup_heavy", args = list(dist_type = "levenshtein", threshold = 2, filter.v = TRUE)), + list(data = "dup_heavy", args = list(dist_type = "levenshtein", threshold = 0.15, normalize = "length", filter.v = TRUE)), + list(data = "equal_len", args = list(dist_type = "hamming", threshold = 2, filter.v = TRUE)), + list(data = "dup_heavy", args = list(dist_type = "damerau", threshold = 2, filter.v = TRUE, filter.j = TRUE)), + list(data = "unique_mixed", args = list(dist_type = "nw", dist_mat = "BLOSUM62", threshold = 5, filter.v = TRUE)) + ) + + for (cfg in cfgs) { + d <- datasets[[cfg$data]] + base_args <- c(cfg$args, list(input.data = d, seq_col = "junction_aa", + v_col = "v_call", j_col = "j_call", + ids = d$sequence_id)) + e_clique <- suppressMessages(do.call(buildNetwork, c(base_args, list(expand = "clique")))) + e_star <- suppressMessages(do.call(buildNetwork, c(base_args, list(expand = "star")))) + + # vertex sets identical + v_clique <- sort(unique(c(e_clique$from, e_clique$to))) + v_star <- sort(unique(c(e_star$from, e_star$to))) + expect_identical(v_star, v_clique, + info = paste(cfg$data, "| vertex set")) + + # connected-component partition identical + expect_identical(component_sets(e_star), component_sets(e_clique), + info = paste(cfg$data, "| component partition")) + + # star is never larger than clique + expect_lte(nrow(e_star), nrow(e_clique)) + } +})