diff --git a/DESCRIPTION b/DESCRIPTION index 28c95bc..a91108d 100644 --- a/DESCRIPTION +++ b/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"))) diff --git a/NEWS.md b/NEWS.md index a631ee0..0fccbc8 100644 --- a/NEWS.md +++ b/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 diff --git a/R/aaa-constants.R b/R/aaa-constants.R new file mode 100644 index 0000000..79f5ed6 --- /dev/null +++ b/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") diff --git a/R/calculateProperty.R b/R/calculateProperty.R index 801baeb..23c9d90 100644 --- a/R/calculateProperty.R +++ b/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) diff --git a/R/sequenceDecoder.R b/R/sequenceDecoder.R index 1fb4f85..ae94425 100644 --- a/R/sequenceDecoder.R +++ b/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) } diff --git a/R/sequenceEncoder.R b/R/sequenceEncoder.R index 1a9dd0c..d1b6b56 100644 --- a/R/sequenceEncoder.R +++ b/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)) { diff --git a/R/utils.R b/R/utils.R index 5cca9e5..8e19819 100644 --- a/R/utils.R +++ b/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 diff --git a/man/amino.acids.Rd b/man/amino.acids.Rd index 05f04ec..0deecdb 100644 --- a/man/amino.acids.Rd +++ b/man/amino.acids.Rd @@ -1,5 +1,5 @@ % Generated by roxygen2: do not edit by hand -% Please edit documentation in R/utils.R +% Please edit documentation in R/aaa-constants.R \name{amino.acids} \alias{amino.acids} \title{Standard 20 amino acids} diff --git a/man/sequenceDecoder.Rd b/man/sequenceDecoder.Rd index 629bfce..a6aae84 100644 --- a/man/sequenceDecoder.Rd +++ b/man/sequenceDecoder.Rd @@ -26,8 +26,11 @@ This is typically inferred if `encoded.object` is a list from `sequenceEncoder`. names (e.g., `"atchleyFactors"`) that were used for the original encoding. See `?sequenceEncoder`. This is ignored if `property.matrix` is supplied.} -\item{property.matrix}{For `mode = "property"`, the exact numeric matrix -(with dimensions `20 x P`) that was used for encoding. This overrides +\item{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`.} \item{call.threshold}{A numeric confidence threshold for making a call. diff --git a/tests/testthat/test-buildNetwork-star.R b/tests/testthat/test-buildNetwork-star.R index b458177..a2a255d 100644 --- a/tests/testthat/test-buildNetwork-star.R +++ b/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")) diff --git a/tests/testthat/test-calculateProperty.R b/tests/testthat/test-calculateProperty.R index 452ef42..e0347cc 100644 --- a/tests/testthat/test-calculateProperty.R +++ b/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", { diff --git a/tests/testthat/test-sequenceDecoder.R b/tests/testthat/test-sequenceDecoder.R index 96db487..1386e81 100644 --- a/tests/testthat/test-sequenceDecoder.R +++ b/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, diff --git a/tests/testthat/test-sequenceEncoder.R b/tests/testthat/test-sequenceEncoder.R index d5ddb21..1ba1d63 100644 --- a/tests/testthat/test-sequenceEncoder.R +++ b/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")