Fix AA-identity/property-value mismatch in sequenceEncoder()/sequenceDecoder() property mode

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: immApex
 Title: Tools for Adaptive Immune Receptor Sequence-Based Machine and Deep Learning
-Version: 1.7.1
+Version: 1.7.2
 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")))

NEWS.md

@@ -1,3 +1,16 @@
+# immApex VERSION 1.7.2
+
+## BUG FIXES
+* Fixed `sequenceEncoder()` and `sequenceDecoder()` property mode silently mismatching amino-acid identity to property values for `crucianiProperties`, `MSWHIM`, and `ProtFP`. These property sets' column order in `Peptides::AAdata` is not already canonical (MSWHIM/ProtFP are stored alphabetically; crucianiProperties has an isolated E/Q swap), and `.aa.property.matrix()` was returning that raw storage order instead of reordering to match the canonical amino-acid order its callers assume positionally. For example, encoding "R" under MSWHIM silently returned Cysteine's values. `.aa.property.matrix()` now always returns canonical-AA-ordered columns, matching `calculateProperty()`'s existing contract and the documented `sequenceEncoder()`/`sequenceDecoder()` interface. Verified against `Peptides`' own independent public scoring functions (`mswhimScores()`, `crucianiProperties()`, `protFP()`).
+* Fixed `sequenceEncoder()` and `sequenceDecoder()` pairing residues with the wrong property values when `property.set` was combined with a non-default `sequence.dictionary`. Both aligned the property matrix to the canonical 20 amino acids regardless of the dictionary actually in use, and neither checked the two agreed. `.aa.property.matrix()` now takes the `sequence.dictionary` and aligns columns to it by name, which also supports restricted alphabets. Requesting a residue a scale has no values for is now an error naming that residue, rather than a silently wrong matrix (this also covers `"pK"`, which exists in `Peptides::AAdata` but only defines 9 residues).
+
+## DOCUMENTATION
+* Corrected `sequenceDecoder()`'s `property.matrix` documentation, which described the matrix as `20 x P`. It is `P x 20`, the transpose of `sequenceEncoder()`'s `property.matrix` argument.
+
+## UNDERLYING CHANGES
+* `amino.acids` was defined twice, in `R/utils.R` and again in `R/calculateProperty.R`. Consolidated to a single definition in `R/aaa-constants.R`, which collates before both.
+* Fixed `test-buildNetwork-star.R` comparing igraph's arbitrary component ID labels rather than the component partition itself. The IDs depend on edge row order, which is not stable between `expand = "clique"` and `expand = "star"`, so the test failed intermittently depending on merge order.
+
 # immApex VERSION 1.7.1
 
 ## PERFORMANCE

R/aaa-constants.R

@@ -0,0 +1,13 @@
+# Package-level constants.
+#
+# This file is named `aaa-constants.R` so it collates first. The built-in
+# property scales in `calculateProperty.R` reference `amino.acids` at load
+# time (in their `dimnames`), and the package has no `Collate:` field, so R
+# sources files in alphabetical order -- `calculateProperty.R` before
+# `utils.R`. Keeping the definition here lets there be exactly one copy.
+
+#' Standard 20 amino acids
+#'
+#' Vector of one-letter codes for the 20 standard amino acids.
+#' @export
+amino.acids <- c("A", "R", "N", "D", "C", "Q", "E", "G", "H", "I", "L", "K", "M", "F", "P", "S", "T", "W", "Y", "V")

R/calculateProperty.R

@@ -68,6 +68,8 @@ calculateProperty <- function(input.sequences,
     stop("'property.set' must be a recognised name or a numeric matrix.")
   }
 
+  # `.aa.property.matrix()` already returns canonical order; this only has to
+  # re-order a user-supplied `property.set` matrix.
   S <- S[ , amino.acids, drop = FALSE]                 # enforce AA order
   k <- nrow(S)
 
@@ -119,18 +121,49 @@ calculateProperty <- function(input.sequences,
   Summ
 }
 
-.aa.property.matrix <- function(key) {
-
+# Align a property scale's columns to `sequence.dictionary`, by name.
+#
+# Every consumer of a property matrix (calculateProperty()'s matrix product,
+# sequenceEncoder()/sequenceDecoder()'s C++ backend) indexes its columns
+# *positionally* against the alphabet, so column i must BE alphabet residue i.
+# Doing that lookup by name here is what makes that contract hold, rather than
+# trusting whatever order the scale happened to be stored in.
+.align.property.columns <- function(v, sequence.dictionary, key) {
+  cn <- colnames(v)
+  if (is.null(cn))
+    stop("Property set '", key, "' has no column names, so it cannot be ",
+         "aligned to 'sequence.dictionary'.", call. = FALSE)
+
+  unknown <- setdiff(sequence.dictionary, cn)
+  if (length(unknown))
+    stop("Property set '", key, "' has no values for the following ",
+         "'sequence.dictionary' entries: ",
+         paste(unknown, collapse = ", "),
+         ". Built-in property scales cover only the 20 canonical amino acids; ",
+         "supply 'property.matrix' directly to use a different alphabet.",
+         call. = FALSE)
+
+  v[, sequence.dictionary, drop = FALSE]
+}
+
+.aa.property.matrix <- function(key, sequence.dictionary = amino.acids) {
+
   if (exists(key, envir = .builtin_scales, inherits = FALSE))
-    return(.builtin_scales[[key]])
-  
+    return(.align.property.columns(.builtin_scales[[key]], sequence.dictionary, key))
+
   has_peptides <- requireNamespace("Peptides", quietly = TRUE)
-  
+
   if (has_peptides) {
-    acc <- utils::getFromNamespace("AAdata", "Peptides")  
+    acc <- utils::getFromNamespace("AAdata", "Peptides")
     if (key %in% names(acc)) {
       v <- do.call(rbind, acc[[key]])
-      return(v)
+      # Columns of `v` come from Peptides::AAdata in whatever order that
+      # scale's author originally stored it (e.g. alphabetical for MSWHIM/
+      # ProtFP), NOT necessarily the canonical `amino.acids` order every
+      # caller (sequenceEncoder/sequenceDecoder's C++ backend, indexed
+      # purely positionally) assumes. Reorder here, once, so every caller
+      # of this helper gets canonical order by construction.
+      return(.align.property.columns(v, sequence.dictionary, key))
     }
   }
 
@@ -145,7 +178,6 @@ calculateProperty <- function(input.sequences,
 }
 
 .builtin_scales <- new.env(parent = emptyenv())
-amino.acids <- c("A", "R", "N", "D", "C", "Q", "E", "G", "H", "I", "L", "K", "M", "F", "P", "S", "T", "W", "Y", "V")
 
 .builtin_scales$atchleyFactors <- t(matrix(c(
   # A (Alanine)

R/sequenceDecoder.R

@@ -11,8 +11,11 @@
 #' @param property.set For `mode = "property"`, a character vector of property
 #'   names (e.g., `"atchleyFactors"`) that were used for the original encoding.
 #'   See `?sequenceEncoder`. This is ignored if `property.matrix` is supplied.
-#' @param property.matrix For `mode = "property"`, the exact numeric matrix
-#'   (with dimensions `20 x P`) that was used for encoding. This overrides
+#' @param property.matrix For `mode = "property"`, the numeric matrix of the
+#'   `P` property scales used for encoding, with dimensions `P x 20`: one row
+#'   per property and one column per entry of `sequence.dictionary`, in that
+#'   order. Note this is the transpose of `sequenceEncoder()`'s
+#'   `property.matrix` argument, which is `20 x P`. This overrides
 #'   `property.set`.
 #' @param call.threshold A numeric confidence threshold for making a call.
 #'   - In `"onehot"` mode, this is the minimum required value in the vector (e.g., `0.9`).
@@ -80,7 +83,7 @@ sequenceDecoder <- function(encoded.object,
       if (!is.null(property.matrix)) {
         depth <- nrow(property.matrix)
       } else if (!is.null(property.set)) {
-        depth <- nrow(.aa.property.matrix(property.set))
+        depth <- nrow(.aa.property.matrix(property.set, sequence.dictionary))
       } else {
         stop("For flattened matrix input in 'property' mode, supply 'property.set' or 'property.matrix'.")
       }
@@ -107,10 +110,14 @@ sequenceDecoder <- function(encoded.object,
       if (is.null(property.set)) {
         stop("In 'property' mode, you must supply either 'property.set' or 'property.matrix'.")
       }
-      property.matrix <- .aa.property.matrix(property.set)
+      # Aligned to `sequence.dictionary` for the same reason as in
+      # sequenceEncoder(): .propertyDecoder() calls the residue by taking
+      # `sequence.dictionary[which.min(distances)]`, so column i of this
+      # matrix must be dictionary residue i.
+      property.matrix <- .aa.property.matrix(property.set, sequence.dictionary)
     }
     if (ncol(property.matrix) != length(sequence.dictionary)) {
-      stop("Rows in 'property.matrix' must match the length of 'sequence.dictionary'.")
+      stop("Columns in 'property.matrix' must match the length of 'sequence.dictionary'.")
     }
     decoded_sequences <- .propertyDecoder(cube, property.matrix, sequence.dictionary, padding.symbol, call.threshold)
   }

R/sequenceEncoder.R

@@ -129,7 +129,11 @@ sequenceEncoder <- function(input.sequences,
   prop_mat <- NULL
   if (mode == "property") {
     if (!is.null(property.set)) {
-      prop_mat <- t(.aa.property.matrix(property.set))
+      # Align to `sequence.dictionary`, not to the canonical 20: the C++
+      # backend indexes `prop_mat` rows positionally against `alphabet`, so a
+      # non-default dictionary would otherwise silently pair each residue with
+      # another residue's property values.
+      prop_mat <- t(.aa.property.matrix(property.set, sequence.dictionary))
     } else if (!is.null(property.matrix)) {
       if (!is.matrix(property.matrix) || !is.numeric(property.matrix) ||
           nrow(property.matrix) != length(sequence.dictionary)) {

R/utils.R

@@ -1,11 +1,5 @@
 `%||%` <- function(x, y) if (is.null(x)) y else x   # tiny helper
 
-#' Standard 20 amino acids
-#'
-#' Vector of one-letter codes for the 20 standard amino acids.
-#' @export
-amino.acids <- c("A", "R", "N", "D", "C", "Q", "E", "G", "H", "I", "L", "K", "M", "F", "P", "S", "T", "W", "Y", "V")
-
 "%!in%" <- Negate("%in%")
 
 # Hoist helper function from the original geometricEncoder.R

tests/testthat/test-buildNetwork-star.R

@@ -11,13 +11,23 @@ test_that("star expansion preserves connected components and vertex set", {
 
   # Each connected component as a sorted vertex-set signature; the set of these
   # signatures is the partition, compared order-independently.
+  # NOTE: the integer component IDs igraph assigns are arbitrary labels that
+  # depend on edge row order, and buildNetwork()'s row order is not stable
+  # between expand modes (OpenMP merge order). vapply() carries those IDs
+  # through as names, so they must be dropped -- otherwise expect_identical()
+  # compares the labels, not the partition, and fails whenever clique and star
+  # happen to number the same components differently. Vertex names are sorted
+  # with method = "radix" for the same reason .bn_canon() does: byte order, so
+  # the signature does not depend on LC_COLLATE.
   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")
+    sigs <- vapply(groups,
+                   function(v) paste(sort(v, method = "radix"), collapse = "\031"),
+                   character(1), USE.NAMES = FALSE)
+    sort(sigs, method = "radix")
   }
 
   cfgs <- list(
@@ -37,8 +47,8 @@ test_that("star expansion preserves connected components and vertex set", {
     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)))
+    v_clique <- sort(unique(c(e_clique$from, e_clique$to)), method = "radix")
+    v_star   <- sort(unique(c(e_star$from, e_star$to)), method = "radix")
     expect_identical(v_star, v_clique,
                      info = paste(cfg$data, "| vertex set"))
 

tests/testthat/test-calculateProperty.R

@@ -76,6 +76,52 @@ test_that("custom property matrix input works", {
   expect_equal(rownames(res), c("S1", "S2"))
 })
 
+test_that(".aa.property.matrix always returns canonical amino-acid column order", {
+  skip_if_not_installed("Peptides")
+
+  # Regression test: sequenceEncoder()/sequenceDecoder() assume column i of
+  # this helper's output IS canonical-order amino acid i, positionally (see
+  # sequenceEncoder.R's roxygen contract). Peptides::AAdata stores some
+  # scales in a different order (MSWHIM/ProtFP: alphabetical; crucianiProperties:
+  # E/Q swapped), so this helper must always reorder before returning -- not
+  # rely on each caller to do it (calculateProperty() used to be the only one
+  # that did).
+  peptides_sets <- c("crucianiProperties", "FASGAI", "kideraFactors", "MSWHIM",
+                     "ProtFP", "stScales", "tScales", "VHSE", "zScales")
+  for (key in peptides_sets) {
+    expect_equal(colnames(.aa.property.matrix(key)), amino.acids, info = key)
+  }
+  expect_equal(colnames(.aa.property.matrix("atchleyFactors")), amino.acids)
+})
+
+test_that(".aa.property.matrix aligns columns to a custom sequence.dictionary", {
+  skip_if_not_installed("Peptides")
+
+  # The 20 canonical residues in a non-canonical order: alignment must follow
+  # the dictionary, not silently fall back to `amino.acids`.
+  dict <- rev(amino.acids)
+  m <- .aa.property.matrix("MSWHIM", dict)
+  expect_identical(colnames(m), dict)
+  expect_equal(unname(m[, "R"]), unname(Peptides::mswhimScores("R")[[1]]))
+
+  # A restricted alphabet is a legitimate subset, not an error.
+  expect_identical(colnames(.aa.property.matrix("MSWHIM", c("A", "C", "R"))),
+                   c("A", "C", "R"))
+
+  # A residue the scale has no values for must fail loudly, naming the residue.
+  expect_error(.aa.property.matrix("MSWHIM", c("A", "X", "Z")), "X, Z")
+
+  # `pK` lives in Peptides::AAdata but covers only 9 residues; it used to fall
+  # through as a silently wrong 9-column matrix.
+  expect_error(.aa.property.matrix("pK"), "no values for")
+})
+
+test_that("built-in scales are stored in canonical amino-acid order", {
+  # Guards the single-definition/collation arrangement of `amino.acids`:
+  # `.builtin_scales` builds its dimnames from it at load time.
+  expect_identical(colnames(.builtin_scales$atchleyFactors), amino.acids)
+})
+
 #  ── Error handling ───────────────────────────────────────────────
 
 test_that("invalid inputs trigger errors", {

tests/testthat/test-sequenceDecoder.R

@@ -87,6 +87,39 @@ test_that(".propertyDecoder decodes standard sequences correctly", {
   )
 })
 
+test_that(".propertyDecoder assigns amino-acid identity correctly (regression)", {
+  skip_if_not_installed("Peptides")
+
+  # Built directly from Peptides' own ground-truth values, NOT via
+  # sequenceEncoder() -- a pure encode+decode round trip through immApex's
+  # own (formerly buggy) functions self-cancels this exact bug, since both
+  # sides would apply the same wrong order. This is the decoder counterpart
+  # to test-sequenceEncoder.R's "property mode assigns amino-acid identity
+  # correctly" regression test: .propertyDecoder used
+  # sequence.dictionary[which.min(distances)] against .aa.property.matrix()'s
+  # unreordered columns, so it could silently call the wrong amino acid.
+  true_R <- unname(Peptides::mswhimScores("R")[[1]])
+  cube <- array(true_R, dim = c(3, 1, 1))
+
+  decoded <- sequenceDecoder(cube, mode = "property", property.set = "MSWHIM",
+                             call.threshold = 0.01)
+  expect_equal(decoded, "R")
+})
+
+test_that(".propertyDecoder honours a non-default sequence.dictionary", {
+  skip_if_not_installed("Peptides")
+
+  # Counterpart to the encoder test: .propertyDecoder() calls the residue via
+  # `sequence.dictionary[which.min(distances)]`, so the property matrix has to
+  # be aligned to that dictionary rather than to canonical `amino.acids`.
+  dict <- rev(amino.acids)
+  cube <- array(unname(Peptides::mswhimScores("R")[[1]]), dim = c(3, 1, 1))
+  expect_equal(
+    sequenceDecoder(cube, mode = "property", property.set = "MSWHIM",
+                    sequence.dictionary = dict, call.threshold = 0.01),
+    "R")
+})
+
 test_that(".propertyDecoder handles padding and thresholds", {
   sequences <- c("CA", "R")
   encoded <- sequenceEncoder(sequences, 

tests/testthat/test-sequenceEncoder.R

@@ -57,6 +57,53 @@ test_that("property encoding with explicit matrix works and returns summary", {
   expect_equal(res$flattened, res.w$flattened)
 })
 
+test_that("property mode assigns amino-acid identity correctly (regression)", {
+  skip_if_not_installed("Peptides")
+
+  # sequenceEncoder() must assign each amino acid ITS OWN property values --
+  # not some other amino acid's, mismatched via positional row order. Checked
+  # against Peptides' own independent public scoring functions (not
+  # immApex/calculateProperty() internals), so this doesn't just re-check
+  # immApex against itself. Regression test for a bug where
+  # crucianiProperties/MSWHIM/ProtFP (whose column order in Peptides::AAdata
+  # is not already canonical) were passed through unreordered -- e.g.
+  # querying "R" under MSWHIM silently returned Cysteine's values.
+  check_single_residue <- function(aa, property.set, truth.fun) {
+    enc <- sequenceEncoder(aa, mode = "property", property.set = property.set,
+                           max.length = 1, verbose = FALSE)$flattened
+    truth <- truth.fun(aa)[[1]]
+    expect_equal(as.numeric(enc[1, paste0(names(truth), "_1")]), unname(truth),
+                tolerance = 1e-8, info = paste(property.set, aa))
+  }
+
+  check_single_residue("R", "MSWHIM", Peptides::mswhimScores)
+  check_single_residue("R", "ProtFP", Peptides::protFP)
+  # crucianiProperties' bug was an isolated E/Q swap -- check both explicitly
+  check_single_residue("E", "crucianiProperties", Peptides::crucianiProperties)
+  check_single_residue("Q", "crucianiProperties", Peptides::crucianiProperties)
+})
+
+test_that("property mode honours a non-default sequence.dictionary", {
+  skip_if_not_installed("Peptides")
+
+  # The C++ backend indexes property-matrix rows positionally against
+  # `alphabet`. Before this was fixed, `property.set` was always aligned to
+  # canonical `amino.acids`, so any other dictionary paired each residue with
+  # a different residue's values -- with no error raised.
+  dict <- rev(amino.acids)
+  enc <- sequenceEncoder("R", mode = "property", property.set = "MSWHIM",
+                         max.length = 1, sequence.dictionary = dict,
+                         verbose = FALSE)$flattened
+  expect_equal(as.numeric(enc[1, ]),
+               unname(Peptides::mswhimScores("R")[[1]]), tolerance = 1e-8)
+
+  expect_error(
+    sequenceEncoder("R", mode = "property", property.set = "MSWHIM",
+                    max.length = 1, sequence.dictionary = c(amino.acids, "X"),
+                    verbose = FALSE),
+    "no values for")
+})
+
 test_that("geometric encoding returns correct shape and content", {
 
   seqs <- c("CARDRST", "YYYGMD", "ACACACAC")