fix bug that caused inf estimates for Ri or Ci ranges

R/06_optimization.R

@@ -174,7 +174,7 @@
 #' x is a vector containing first the RiLog values and then the Ci values
 #' (length 2*m if no Ci are fixed).
 #'
-#' For the equality constraint (eqfun with corresponding bound eqB), the 
+#' For the equality constraint (eqfun with corresponding bound eqB), the
 #' difference of the sum of the Ci values and 1 is calculated. The bound is set
 #' to 0, i.e. forcing the sum of the Ci to be 1. In case of fixed Ci values, the
 #' sum of the fixed and the free Ci values is considered.
@@ -460,8 +460,8 @@
 #' G <- graphs$sample1_sample2[[2]]
 #' # small example with a small grid size
 #' iterateOverCi(G, gridSize = 100)
-#' 
-#' @importFrom checkmate assertClass assertFlag assertIntegerish assertList 
+#'
+#' @importFrom checkmate assertClass assertFlag assertIntegerish assertList
 #'  assertNumeric checkTRUE
 #' @importFrom igraph is_bipartite V
 #' @importFrom  pbapply pbmapply pboptions
@@ -524,12 +524,12 @@ iterateOverCi <- function(G,
 #'                                      The data.frame resulting from the
 #'                                      [bppg::iterateOverCi()] function,
 #'                                      filtered for a specific protein.
-#' @param error_tol \strong{numeric(1)} \cr 
+#' @param error_tol \strong{numeric(1)} \cr
 #'                  tolerance for the error term.
-#' @param ratioLog_tol \strong{numeric(1)} \cr 
+#' @param ratioLog_tol \strong{numeric(1)} \cr
 #'                      tolerance for the log protein ratios.
 #'
-#' @returns Vector with minimal error, estimate for Ri (single value 
+#' @returns Vector with minimal error, estimate for Ri (single value
 #' or min/max), estimate for Ci (single value or min/max).
 #'
 .analyseResultSingleProt <- function(protNr, resProt,
@@ -553,7 +553,7 @@ iterateOverCi <- function(G,
             res_tmp <- c(mean(e_tmp), RLog_tmp_mean_rounded, NA, NA, NA,
                 min(C_tmp), max(C_tmp), NA)
             # all(R[ind_min_tol] == 0) |???
-            res_tmp[8] <- ifelse(length(indMinError) == nrow(resProt), 1, 4) 
+            res_tmp[8] <- ifelse(length(indMinError) == nrow(resProt), 1, 4)
         } else {
             # case 2 and 5: range solution for Rlog
             res_tmp <- c(mean(e_tmp), NA, min(RLog_tmp), max(RLog_tmp), NA,
@@ -604,7 +604,7 @@ iterateOverCi <- function(G,
 #'                                      Tolerance for a constant log2 protein
 #'                                      ratios. The default is 1e-6.
 #'
-#' @return A SummarizedExperiment object with one row for each protein in the 
+#' @return A SummarizedExperiment object with one row for each protein in the
 #'  assay.
 #' @export
 #'
@@ -638,7 +638,18 @@ automatedAnalysisIteratedCi <- function(G,
         job.id <- NA
     }
 
-    f <- function(x, res, error_tol, ratioLog_tol, 
+    # filter res for potential NaNs in the error column (could occur if a Ci
+    # is estimated as 0)
+    ind_error_NA <- which(is.na(res$error))
+    if (length(ind_error_NA) > 0) {
+        if (verbose) {
+            message(paste0(length(ind_error_NA),
+                           " grid points with NA or NaN error term were removed."))
+        }
+        res <- res[-ind_error_NA, ]
+    }
+
+    f <- function(x, res, error_tol, ratioLog_tol,
         use_results_from_other_proteins) {
 
         cols <- c("protein", "error", paste0("RLog", x), paste0("C", x))
@@ -650,7 +661,7 @@ automatedAnalysisIteratedCi <- function(G,
             resProt2 <- res[res$protein != x, cols]
             colnames(resProt2) <- c("protein", "error", "RLog", "C")
             # remove too extreme Ci
-            resProt2 <- resProt2[resProt2$C > 0.01 | resProt2$C < 0.99,] 
+            resProt2 <- resProt2[resProt2$C > 0.01 & resProt2$C < 0.99,]
             resProt <- rbind(resProt, resProt2)
         }
         return(.analyseResultSingleProt(x, resProt, error_tol = error_tol,
@@ -664,15 +675,15 @@ automatedAnalysisIteratedCi <- function(G,
         use_results_from_other_proteins = use_results_from_other_proteins)
 
     # comparison=rep(comparison, n),
-    RES_info <- data.frame(accession = accessions, 
+    RES_info <- data.frame(accession = accessions,
         graphID = rep(graphID, n), proteinNr = seq_len(n))
 
     RES <- cbind(RES_info, as.data.frame(t(RES)))
     rownames(RES) <- RES$accession
     RES <- RES[, -1] # remove accession, as it is now in rownames
 
     RES_SE <- SummarizedExperiment::SummarizedExperiment(
-        assays = list(results = RES), 
+        assays = list(results = RES),
         rowData = data.frame(accession = accessions),
         colData = data.frame(colnames = colnames(RES)))
     return(RES_SE)

tests/testthat/_snaps/06_optimization.md

@@ -66,3 +66,15 @@
       prot_3           0.1    0.1    4
       prot_4;prot_5    0.9    0.9    4
 
+---
+
+    Code
+      SummarizedExperiment::assay(res5)
+    Output
+                    graphID proteinNr   error_min     RiLog RiLog_min RiLog_max  Ci
+      prot_3             NA         1 0.006932845 0.9779656        NA        NA 0.9
+      prot_4;prot_5      NA         2 0.006635250 0.9835544        NA        NA 0.9
+                    Ci_min Ci_max case
+      prot_3            NA     NA    3
+      prot_4;prot_5     NA     NA    3
+

tests/testthat/test_06_optimization.R

@@ -129,6 +129,12 @@ test_that("test automated analysis", {
   # 2         NA      NA         2 0.00663525 0.983554        NA        NA NA    0.9    0.9    4
   # TODO: it switches to case 4 because for both proteins an almost identical result is achieved
 
+  ## test case when an error term is NaN (which may happen during the optimization,
+  ## if a Ci is estimated as 0)
+  res_snap$error[c(1,5)] <- NaN
+  res5 <- automatedAnalysisIteratedCi(G_M, res_snap)
+  expect_snapshot(SummarizedExperiment::assay(res5))
+
 })