add new outlier method

DESCRIPTION

@@ -69,13 +69,13 @@ Imports:
     GSEABase,
     GSVA,
     harmony,
-    HiClimR,
     htmlwidgets,
     igraph,
     iheatmapr,
     IntNMF,
     IRanges,
     irlba,
+    isotree,
     isva,
     jsonlite,
     karyoploteR,

R/pgx-outlier.R

@@ -4,59 +4,112 @@
 ##
 
 #' @export
-detectOutlierSamples <- function(X, plot = TRUE, par = NULL) {
+detectOutlierSamples <- function(X,
+                                 methods = c("z.correlation", "z.distance",
+                                   "z.features", "z.isoforest")[1:3],
+                                 col = "grey70", plot = TRUE, par = TRUE) {
+
+  all.methods <- c("z.correlation", "z.distance", "z.features", "z.isoforest")
+  if (is.null(methods)) methods <- all.methods
+  ## errors on a typo. without this an unmatched name silently yields a
+  ## zero-column Z and NaN z-scores, which read as "no outliers".
+  methods <- match.arg(methods, all.methods, several.ok = TRUE)
+
   ## correlation and distance
   X <- head(X[order(-matrixStats::rowSds(X, na.rm = TRUE)), ], 1000)
   X <- X - median(X, na.rm = TRUE)
-  corX <- HiClimR::fastCor(X, optBLAS = TRUE)
+  #corX <- HiClimR::fastCor(X, optBLAS = TRUE)
+  corX <- cor(X, use="pairwise.complete")  
   distX <- as.matrix(dist(t(X)))
 
+  z1=z2=z3=z4=NULL
+
+  ## robust z-score. mad() is 0 as soon as over half the values tie (e.g. a
+  ## sample uploaded twice), which used to give 0/0 = NaN for every sample and
+  ## crash the caller. Fall back to the mean absolute deviation, and to zeros
+  ## when every value really is identical (nothing deviates, so nothing is an
+  ## outlier). Unchanged whenever mad() > 0.
+  zscore <- function(x) {
+    dx <- abs(x - median(x, na.rm = TRUE))
+    s <- mad(x, na.rm = TRUE)
+    if (s == 0) s <- mean(dx, na.rm = TRUE) / 0.7979 ## consistency constant
+    if (s == 0) return(dx * 0)
+    dx / s
+  }
+
   ## z-score based on correlation
-  ## cor.min
-  ## cor.max <- apply(abs(corX), 1, max, na.rm = TRUE)
   cor.median <- apply(abs(corX), 1, median, na.rm = TRUE)
-  ## cor.q10 <- apply(abs(corX), 1, quantile, probs = 0.1, na.rm = TRUE)
   x1 <- (cor.median - mean(cor.median, na.rm = TRUE))
-  z1 <- abs(x1 - median(x1, na.rm = TRUE)) / mad(x1, na.rm = TRUE)
+  z1 <- zscore(x1)
 
   ## z-score based on euclidean distance
   dist.max <- apply(distX, 1, max, na.rm = TRUE)
-  ## dist.median <- apply(distX, 1, median, na.rm = TRUE)
-  ## dist.q90 <- apply(distX, 1, quantile, probs = 0.9, na.rm = TRUE)
   dist.q10 <- apply(distX, 1, quantile, probs = 0.1, na.rm = TRUE)
   dist.r <- dist.q10 / dist.max
-  z2 <- abs(dist.r - median(dist.r, na.rm = TRUE)) / mad(dist.r, na.rm = TRUE)
+  z2 <- zscore(dist.r)
 
   ## gene-wise z-score
   xz <- abs(X - rowMeans(X, na.rm = TRUE)) / matrixStats::rowSds(X, na.rm = TRUE)
   xz <- colMeans(xz, na.rm = TRUE)
-  z3 <- abs(xz - median(xz, na.rm = TRUE)) / mad(xz, na.rm = TRUE)
+  z3 <- zscore(xz)
+
+  ## isoforest z-score. only on request: fitting 10k trees is expensive
+  if ("z.isoforest" %in% methods) {
+    z4 <- outlier.isoforest_zscore(X, ndim=2, ntrees=10000)
+  }
+
+  ## NULL columns are dropped by cbind(), so unrequested methods vanish here
+  Z <- cbind(z.correlation = z1, z.distance = z2, z.features = z3,
+    z.isoforest = z4)
+  Z <- Z[, which(colnames(Z) %in% methods), drop = FALSE]
 
-  Z <- cbind(z1, z2, z3)
-  colnames(Z) <- c("z.correlation", "z.distance", "z.features")
   zz <- rowMeans(Z, na.rm = TRUE)
   z0 <- 0.1 * mean(Z, na.rm = TRUE)
   zz2 <- exp(rowMeans(log(Z + z0), na.rm = TRUE)) - z0
 
   res <- list(z.outlier = zz, z.outlier2 = zz2, Z = Z)
-  if (plot) plotOutlierScores(res, par = par)
+  if (plot) plotOutlierScores(res, par = par, col=col)
   return(res)
 }
 
 #' @export
-plotOutlierScores <- function(res.outliers, z.threshold = c(3, 6, 9), par = TRUE) {
-  if (par) par(mfrow = c(2, 3), mar = c(8, 4, 2, 2))
+plotOutlierScores <- function(res.outliers, z.threshold = c(3, 6, 9),
+                              col = "grey70", par = TRUE) {
+  if (par) {
+    ## restore the caller's layout, we are a guest in their device
+    opar <- graphics::par(mfrow = c(2, 3), mar = c(8, 4, 2, 2))
+    on.exit(graphics::par(opar))
+  }
   Z <- res.outliers$Z
   zz <- res.outliers$z.outlier
   zz2 <- res.outliers$z.outlier2
   barplot2 <- function(x, ...) {
-    barplot(x, ylim = c(0, max(10, max(Z))), ylab = "z-score", ...)
+    barplot(x, col = col, las = 3,
+      ## finite only: an all-NA/Inf column gave "need finite 'ylim' values".
+      ## max(10, numeric(0)) is 10, so a fully non-finite Z still plots.
+      ylim = c(0, max(10, Z[is.finite(Z)])),
+      ylab = "z-score", ...)
     abline(h = z.threshold, lty = 3, col = "red")
   }
-  barplot2(zz, main = "z.outlier (mean)", las = 3)
-  barplot2(zz2, main = "z.outlier (geom.mean)", las = 3)
-  for (i in 1:ncol(Z)) {
+  barplot2(zz, main = "z.outlier (mean)")
+  barplot2(zz2, main = "z.outlier (geom.mean)")
+  for (i in seq_len(ncol(Z))) {
     z1 <- Z[, i]
-    barplot2(z1, main = colnames(Z)[i], las = 3)
+    barplot2(z1, main = colnames(Z)[i])
   }
 }
+
+#' 
+outlier.isoforest_zscore <- function(X, ndim=2, ntrees=10000) {
+  cX <- X - rowMeans(X, na.rm=TRUE)
+  cX <- cX[complete.cases(cX),,drop=FALSE]
+  ## ponytail: plain svd. callers cap cX at 1000 rows so a truncated
+  ## solver buys nothing, and irlba(nv=3) aborted below 4 samples.
+  V <- svd(cX, nu = 0, nv = min(3, ncol(cX)))$v
+  ndim <- min(ndim, ncol(V))
+  model <- isotree::isolation.forest(V, ndim=ndim, ntrees=ntrees, nthreads=8)
+  scores <- predict(model, V)
+  scores <- abs(scores - mean(scores, na.rm=TRUE))
+  scores <- scores / sd(scores, na.rm=TRUE)  
+  scores 
+}

R/pgx-preprocess.R

@@ -33,6 +33,9 @@
 #'     \item{impute_method}{Imputation method passed to `imputeMissing`. Default "SVD2".}
 #'     \item{remove_outliers}{Drop outlier samples. Default FALSE.}
 #'     \item{outlier_threshold}{z-score cutoff for `detectOutlierSamples`. Default 3.}
+#'     \item{outlier_methods}{z-score methods for `detectOutlierSamples`: any of
+#'       "z.correlation", "z.distance", "z.features", "z.isoforest". Default the
+#'       first three; "z.isoforest" is opt-in as it fits an isolation forest.}
 #'     \item{meth_type}{Methylation array type for `normalizeMethylation`. Default NULL.}
 #'   }
 #'
@@ -59,6 +62,9 @@ pgx.preprocess <- function(counts,
       impute_method = "SVD2",
       remove_outliers = FALSE,
       outlier_threshold = 3,
+      ## NB: not NULL. detectOutlierSamples() reads NULL as "all methods",
+      ## which would switch on the isoforest behind the caller's back.
+      outlier_methods = c("z.correlation", "z.distance", "z.features"),
       meth_type = NULL
     ),
     options
@@ -194,8 +200,10 @@ pgx.preprocess <- function(counts,
         X <- playbase::imputeMissing(X, method = "SVD2")
       }
     }
-    res <- playbase::detectOutlierSamples(X, plot = FALSE)
-    is.outlier <- (res$z.outlier > opt$outlier_threshold)
+    res <- playbase::detectOutlierSamples(X, methods = opt$outlier_methods, plot = FALSE)
+    ## NA-safe: a non-finite score must not be read as an outlier, and must not
+    ## reach the if() below as NA ("missing value where TRUE/FALSE needed").
+    is.outlier <- !is.na(res$z.outlier) & (res$z.outlier > opt$outlier_threshold)
     if (any(is.outlier) && !all(is.outlier)) {
       X <- X[, which(!is.outlier), drop = FALSE]
       counts <- counts[, colnames(X), drop = FALSE]

tests/testthat/test-pgx-outlier.R

@@ -0,0 +1,102 @@
+## Covers detectOutlierSamples() method selection and the isoforest path,
+## which pgx.preprocess() enables via opt$outlier_methods.
+
+mk_X <- function(n) {
+  set.seed(1)
+  X <- matrix(rnorm(500 * n), 500, n,
+    dimnames = list(paste0("g", 1:500), paste0("s", 1:n)))
+  X[, n] <- X[, n] + 4 ## last sample is the outlier
+  X
+}
+
+test_that("detectOutlierSamples returns one Z column per requested method", {
+  X <- mk_X(8)
+  res <- playbase::detectOutlierSamples(X, plot = FALSE)
+  expect_equal(colnames(res$Z), c("z.correlation", "z.distance", "z.features"))
+  expect_equal(unname(which.max(res$z.outlier)), 8L)
+
+  ## a single method must stay a matrix, else rowMeans() fails
+  res1 <- playbase::detectOutlierSamples(X, methods = "z.distance", plot = FALSE)
+  expect_true(is.matrix(res1$Z))
+  expect_equal(colnames(res1$Z), "z.distance")
+  expect_length(res1$z.outlier, ncol(X))
+})
+
+test_that("detectOutlierSamples rejects an unknown method", {
+  ## used to return NaN z-scores silently, reading as "no outliers"
+  expect_error(
+    playbase::detectOutlierSamples(mk_X(8), methods = "z.correlaton", plot = FALSE),
+    "should be one of"
+  )
+})
+
+test_that("isoforest is off by default and runs when asked", {
+  skip_if_not_installed("isotree")
+  X <- mk_X(8)
+  expect_false("z.isoforest" %in% colnames(playbase::detectOutlierSamples(X, plot = FALSE)$Z))
+
+  ## NB: structure only. The isoforest's score scale and sign are under review
+  ## (it is two-sided and sd-normalized, unlike the mad-based methods), so this
+  ## deliberately does not assert which sample scores highest.
+  res <- playbase::detectOutlierSamples(X, methods = "z.isoforest", plot = FALSE)
+  expect_equal(colnames(res$Z), "z.isoforest")
+  expect_true(all(is.finite(res$z.outlier)))
+
+  ## used to abort in irlba(nv = 3) with fewer than 4 samples
+  small <- playbase::detectOutlierSamples(mk_X(3), methods = "z.isoforest", plot = FALSE)
+  expect_length(small$z.outlier, 3L)
+})
+
+test_that("the mad-based methods score the injected outlier well above threshold", {
+  res <- playbase::detectOutlierSamples(mk_X(8), plot = FALSE)
+  expect_gt(res$z.outlier[["s8"]], 3)
+  expect_true(all(res$z.outlier[paste0("s", 1:7)] < 3))
+})
+
+test_that("tied samples give finite scores instead of NaN", {
+  ## mad() is 0 when over half the values tie, which used to yield 0/0 for
+  ## every method and crash pgx.preprocess with "missing value where
+  ## TRUE/FALSE needed". Same sample uploaded 3x, twice.
+  set.seed(1)
+  X <- matrix(rnorm(500 * 6), 500, 6,
+    dimnames = list(paste0("g", 1:500), paste0("s", 1:6)))
+  X[, 2] <- X[, 1]; X[, 3] <- X[, 1]
+  X[, 5] <- X[, 4]; X[, 6] <- X[, 4]
+
+  res <- playbase::detectOutlierSamples(X, plot = FALSE)
+  expect_true(all(is.finite(res$Z)))
+  expect_true(all(is.finite(res$z.outlier)))
+
+  samples <- data.frame(group = rep(c("a", "b"), each = 3), row.names = colnames(X))
+  out <- playbase::pgx.preprocess(2^X, samples, contrasts = NULL,
+    options = list(remove_outliers = TRUE, outlier_threshold = 3))
+  expect_equal(ncol(out$X), 6L) ## identical samples, none is an outlier
+})
+
+test_that("plotting works on degenerate input and restores par", {
+  pdf(NULL)
+  on.exit({ dev.off(); unlink("Rplots.pdf") }, add = TRUE)
+  before <- graphics::par("mfrow")
+
+  ## non-finite Z used to abort with "need finite 'ylim' values"
+  res <- playbase::detectOutlierSamples(mk_X(3), plot = TRUE)
+  expect_length(res$z.outlier, 3L)
+  expect_equal(graphics::par("mfrow"), before)
+
+  ## all four methods, 6 panels
+  skip_if_not_installed("isotree")
+  expect_length(playbase::detectOutlierSamples(mk_X(8), methods = NULL, plot = TRUE)$z.outlier, 8L)
+  expect_equal(graphics::par("mfrow"), before)
+})
+
+test_that("pgx.preprocess plumbs outlier_methods to detectOutlierSamples", {
+  X <- mk_X(8)
+  counts <- 2^X
+  samples <- data.frame(group = rep(c("a", "b"), each = 4), row.names = colnames(X))
+  run <- function(...) playbase::pgx.preprocess(counts, samples, contrasts = NULL,
+    options = list(remove_outliers = TRUE, outlier_threshold = 3, ...))
+
+  ## the option must reach detectOutlierSamples()' validation, not be ignored
+  expect_error(run(outlier_methods = "nonsense"), "should be one of")
+  expect_no_error(run(outlier_methods = "z.distance"))
+})