Fix crash from factor sample.class and sparse tmb tables (#11)

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: SelectSim
 Title: Selected Events Linked by Evolutionary Conditions in Cancer
-Version: 0.1.6
+Version: 0.1.7
 Authors@R: c(
     person("Arvind", "Iyer", , "ayalurarvind@gmail.com", role = c("aut", "cre", "cph"),
            comment = c(ORCID = "0000-0002-8247-700X")),

NEWS.md

@@ -1,3 +1,10 @@
+# SelectSim 0.1.7
+
+* Fixed `replacement has length zero` crash in `retrieveOutliers()` caused by a factor `sample.class`/`alteration.class` silently losing sample names in `get.blocks()` (`which()` drops names on a named factor but not a named character vector); factor input is now auto-coerced to character (#11).
+* Fixed silent mis-alignment of per-sample TMB totals in `new.AL.general()` when `am$tmb` tables only listed a subset of samples (e.g. only samples with a nonzero mutation count for that type); TMB tables are now matched and zero-filled by sample identity instead of summed by raw vector position (#11).
+* `new.AL.general()` now errors clearly if an `am$tmb` table references a sample id absent from the corresponding `am$M` matrix, or contains duplicate sample entries.
+* Added regression tests covering factor `sample.class`, partial/sparse `tmb` tables, and unknown sample ids in `tmb`.
+
 # SelectSim 0.1.6
 
 * Removed `reshape2` dependency; replaced with base R equivalents (`tapply`, `as.data.frame.table`).

R/gam_utils.r

@@ -2,7 +2,7 @@
 # Author  : Marco Mina , Arvind Iyer
 # Project : SelectSim
 # Desc    : Functions to process the maf to gam
-# Version : 0.1.6
+# Version : 0.1.7
 ###
 
 

R/selectX_create.r

@@ -2,7 +2,7 @@
 # Author  : Arvind Iyer, Miljan Petrovic
 # Project : SelectSim
 # Desc    : Implementation of SelectSim algorithm
-# Version : 0.1.6
+# Version : 0.1.7
 ###
 
 #' Create an Alteration Landscape (AL) object
@@ -14,13 +14,17 @@
 #'   alteration matrices, each genes x samples) and \code{tmb} (a named list of data
 #'   frames, one per matrix in \code{M}, each with columns \code{sample} and
 #'   \code{mutation}). The names of \code{M} and \code{tmb} must match. All matrices
-#'   in \code{M} must have identical row and column names.
-#' @param feat.covariates Named character vector of alteration-type annotations, one
-#'   entry per feature (gene). Names must match rownames of the matrices in \code{M}.
-#'   If \code{NULL}, all features are labelled \code{"MUT"}.
-#' @param sample.covariates Named character vector of sample-type annotations, one
-#'   entry per sample. Names must match colnames of the matrices in \code{M}. If
-#'   \code{NULL}, all samples are labelled \code{"sample"}.
+#'   in \code{M} must have identical row and column names. Each \code{tmb} table may
+#'   be partial (list only the samples with a nonzero \code{mutation} count for that
+#'   type) — any sample present in \code{M} but absent from a \code{tmb} table is
+#'   treated as having \code{mutation = 0} for that type.
+#' @param feat.covariates Named character (or factor, auto-coerced) vector of
+#'   alteration-type annotations, one entry per feature (gene). Names must match
+#'   rownames of the matrices in \code{M}. If \code{NULL}, all features are labelled
+#'   \code{"MUT"}.
+#' @param sample.covariates Named character (or factor, auto-coerced) vector of
+#'   sample-type annotations, one entry per sample. Names must match colnames of the
+#'   matrices in \code{M}. If \code{NULL}, all samples are labelled \code{"sample"}.
 #' @param min.freq Minimum number of samples a gene must be mutated in (strictly
 #'   greater than) to be retained. Features with \code{rowSums <= min.freq} are
 #'   dropped.
@@ -96,24 +100,53 @@ new.AL.general <- function(am,
     al$alterations$alteration.class <- rep("MUT", nrow(al$am[[1]]))
     names(al$alterations$alteration.class) <- rownames(al$am[[1]])
   } else {
+    if (is.factor(feat.covariates)) {
+      if (verbose) message("feat.covariates was a factor; converting to character.")
+      feat.covariates <- setNames(as.character(feat.covariates), names(feat.covariates))
+    }
     al$alterations$alteration.class <- feat.covariates
   }
   if (is.null(sample.covariates)) {
     al$samples$sample.class <- rep("sample", ncol(al$am[[1]]))
     names(al$samples$sample.class) <- colnames(al$am[[1]])
   } else {
+    if (is.factor(sample.covariates)) {
+      if (verbose) message("sample.covariates was a factor; converting to character.")
+      sample.covariates <- setNames(as.character(sample.covariates), names(sample.covariates))
+    }
     al$samples$sample.class <- sample.covariates
   }
-  # set the tumor mutation burden vector
+  # set the tumor mutation burden vector, aligning every tmb table to col.order
+  # by sample identity (rather than assuming row order/length already match) and
+  # zero-filling any sample missing from a given tmb table.
   al$tmb <- list()
   for (i in names(am$M)) {
-    al$tmb[[i]] <- am$tmb[[i]]
-  }
-  al$tmb[["total"]] <- c(rep(0, ncol(am$M[[1]])))
-  for (i in names(am$tmb)) {
-    al$tmb[["total"]] <- al$tmb[["total"]] + am$tmb[[i]][, c("mutation")]
+    tb <- am$tmb[[i]]
+    tb_sample <- as.character(tb$sample)
+    unknown <- setdiff(tb_sample, col.order)
+    if (length(unknown) > 0) {
+      stop(sprintf(
+        "am$tmb[['%s']] contains sample id(s) not present in am$M[['%s']]: %s",
+        i, i, paste(unknown[seq_len(min(5, length(unknown)))], collapse = ", ")
+      ))
+    }
+    if (anyDuplicated(tb_sample) > 0) {
+      stop(sprintf("am$tmb[['%s']] has duplicate 'sample' entries.", i))
+    }
+    missing <- setdiff(col.order, tb_sample)
+    if (verbose && length(missing) > 0) {
+      message(sprintf(
+        "am$tmb[['%s']]: %d of %d samples had no entry - filled with mutation = 0",
+        i, length(missing), length(col.order)
+      ))
+    }
+    mutation <- setNames(rep(0, length(col.order)), col.order)
+    mutation[tb_sample] <- tb$mutation
+    al$tmb[[i]] <- data.frame(sample = col.order, mutation = mutation,
+                              row.names = col.order, stringsAsFactors = FALSE)
   }
-  names(al$tmb$total) <- am$tmb[[1]]$sample
+  al$tmb[["total"]] <- Reduce(`+`, lapply(al$tmb, function(t) t[col.order, "mutation"]))
+  names(al$tmb$total) <- col.order
 
   class(al) <- "AL"
   return(al)
@@ -144,8 +177,13 @@ get.blocks <- function(al) {
     al$samples$sample.class <- rep("sample", ncol(al$am[[1]]))
     names(al$samples$sample.class) <- colnames(al$am[[1]])
   }
-  alteration.class <- al$alterations$alteration.class[rownames(al$am$full)]
-  sample.class <- al$samples$sample.class[colnames(al$am$full)]
+  # Coerce defensively to character: which() silently drops names when comparing
+  # a named factor (but not a named character vector), which would otherwise
+  # produce unnamed sample/feature blocks downstream.
+  alteration.class <- as.character(al$alterations$alteration.class[rownames(al$am$full)])
+  names(alteration.class) <- rownames(al$am$full)
+  sample.class <- as.character(al$samples$sample.class[colnames(al$am$full)])
+  names(sample.class) <- colnames(al$am$full)
   feature.blocks <- lapply(unique(alteration.class), function(x) which(alteration.class == x))
   names(feature.blocks) <- unique(alteration.class)
   sample.blocks <- lapply(unique(sample.class), function(x) which(sample.class == x))

R/selectX_plot.R

@@ -2,7 +2,7 @@
 # Author  : Arvind Iyer, Miljan Petrovic
 # Project : SelectSim
 # Desc    : The file contains plot related functions
-# Version : 0.1.6
+# Version : 0.1.7
 ###
 
 # Suppress R CMD check notes for ggplot2/ggridges column name variables

R/selectX_run.R

@@ -2,7 +2,7 @@
 # Author  : Arvind Iyer, Miljan Petrovic
 # Project : SelectSim
 # Desc    : Main file which to run the SelecSim algoritm via calling selectX function to create alteration object with background model and funtion to generate the table.
-# Version : 0.1.6
+# Version : 0.1.7
 # Notes:
 # - Better Error message and running text
 # - Edge case: When sample size in less than 2 there is error in computation (need to fix a number to do this analysis)

R/selectX_stats.r

@@ -2,7 +2,7 @@
 # Author  : Arvind Iyer
 # Project : SelectSim
 # Desc    : The file which contains the function to generate the stats and table
-# Version : 0.1.6
+# Version : 0.1.7
 ###
 
 #' Initialize an Alteration Landscape Stats (ALS) container

tests/testthat/test-new.AL.general-robustness.R

@@ -0,0 +1,96 @@
+make_synthetic_am <- function(n_genes = 20, n_samples = 30, seed = 1) {
+  set.seed(seed)
+  samples <- paste0("S", seq_len(n_samples))
+  genes <- paste0("G", seq_len(n_genes))
+  mk_gam <- function(p) {
+    matrix(rbinom(n_genes * n_samples, 1, p), nrow = n_genes,
+           dimnames = list(genes, samples))
+  }
+  M <- list(Nonsense = mk_gam(0.2), Missense = mk_gam(0.1))
+  mk_full_tmb <- function(gam) {
+    csum <- colSums(gam)
+    data.frame(sample = names(csum), mutation = csum, row.names = names(csum))
+  }
+  list(M = M, tmb = list(Nonsense = mk_full_tmb(M$Nonsense), Missense = mk_full_tmb(M$Missense)),
+       samples = samples, genes = genes)
+}
+
+test_that("factor sample.class/alteration.class no longer breaks get.blocks() naming", {
+  fx <- make_synthetic_am()
+  sample.class.chr <- setNames(rep(c("A", "B"), length.out = length(fx$samples)), fx$samples)
+  sample.class.factor <- factor(sample.class.chr)
+  alteration.class <- setNames(rep("MUT", length(fx$genes)), fx$genes)
+
+  al_chr <- new.AL.general(fx, feat.covariates = alteration.class,
+                            sample.covariates = sample.class.chr, min.freq = 1)
+  al_fac <- new.AL.general(fx, feat.covariates = alteration.class,
+                            sample.covariates = sample.class.factor, min.freq = 1)
+
+  blocks_chr <- get.blocks(al_chr)
+  blocks_fac <- get.blocks(al_fac)
+
+  # Previously: factor input made which() drop names, leaving blocks unnamed
+  # (length(names(...)) == 0), which cascaded into empty template matrices.
+  expect_true(all(vapply(blocks_fac$sample.blocks, function(b) length(names(b)) > 0, logical(1))))
+  expect_equal(blocks_chr$sample.blocks, blocks_fac$sample.blocks)
+  expect_equal(al_chr$tmb$total, al_fac$tmb$total)
+})
+
+test_that("selectX() runs to completion with a factor sample.class (regression for #11)", {
+  fx <- make_synthetic_am()
+  sample.class <- factor(setNames(rep(c("BM", "PCT"), length.out = length(fx$samples)), fx$samples))
+  alteration.class <- setNames(rep("MUT", length(fx$genes)), fx$genes)
+
+  expect_no_error({
+    result <- selectX(M = fx, sample.class = sample.class,
+                       alteration.class = alteration.class,
+                       min.freq = 1, n.permut = 10, n.cores = 1, verbose = FALSE)
+  })
+  expect_true(is.list(result))
+})
+
+test_that("partial (sparse) tmb tables are zero-filled by sample identity, not recycled", {
+  fx <- make_synthetic_am()
+  # Drop most Missense rows, keeping only a handful of samples - mirrors the
+  # reporter's tmb tables, which only listed samples with >=1 mutation of that type.
+  sparse_tmb <- fx$tmb$Missense[1:5, ]
+  fx_sparse <- fx
+  fx_sparse$tmb$Missense <- sparse_tmb
+
+  sample.class <- setNames(rep("sample", length(fx$samples)), fx$samples)
+  alteration.class <- setNames(rep("MUT", length(fx$genes)), fx$genes)
+
+  expect_no_warning({
+    al <- new.AL.general(fx_sparse, feat.covariates = alteration.class,
+                          sample.covariates = sample.class, min.freq = 1)
+  })
+
+  # Samples absent from the sparse Missense table should contribute mutation = 0,
+  # not a positionally-recycled value from an unrelated sample.
+  dropped_samples <- setdiff(fx$samples, sparse_tmb$sample)
+  expect_true(all(al$tmb$Missense[dropped_samples, "mutation"] == 0))
+
+  # Total should equal Nonsense (full) + Missense (zero-filled for dropped
+  # samples, actual value otherwise), matched by sample identity.
+  missense_zero_filled <- setNames(rep(0, length(fx$samples)), fx$samples)
+  missense_zero_filled[sparse_tmb$sample] <- sparse_tmb$mutation
+  expected_total <- fx$tmb$Nonsense[fx$samples, "mutation"] + missense_zero_filled[fx$samples]
+  expect_equal(unname(al$tmb$total[fx$samples]), unname(expected_total))
+})
+
+test_that("am$tmb with an unknown sample id errors clearly", {
+  fx <- make_synthetic_am()
+  bad_tmb <- fx$tmb$Nonsense
+  bad_tmb$sample[1] <- "NOT_A_REAL_SAMPLE"
+  rownames(bad_tmb)[1] <- "NOT_A_REAL_SAMPLE"
+  fx$tmb$Nonsense <- bad_tmb
+
+  sample.class <- setNames(rep("sample", length(fx$samples)), fx$samples)
+  alteration.class <- setNames(rep("MUT", length(fx$genes)), fx$genes)
+
+  expect_error(
+    new.AL.general(fx, feat.covariates = alteration.class,
+                    sample.covariates = sample.class, min.freq = 1),
+    "not present in am\\$M"
+  )
+})