bug fixes and improve intitaliser

ococrook · Apr 26, 2024 · 352acbf · 352acbf
1 parent 5df5d4d
commit 352acbf
Show file tree

Hide file tree

Showing 10 changed files with 90 additions and 21 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: RexMS
 Title: Residue-level statistical analysis of HDX-MS data
-Version: 0.99.3
+Version: 0.99.4
 Date: 2024-03-04
 Authors@R: 
     person("Oliver", "Crook", , "oliver.crook@stats.ox.ac.uk", role = c("aut", "cre"),
@@ -47,7 +47,8 @@ Imports:
     vegan,
     plyr,
     ropls,
-    tidyr
+    tidyr,
+    stringr
 Suggests: 
     BiocStyle,
     ggfortify,

diff --git a/NAMESPACE b/NAMESPACE
@@ -2,7 +2,7 @@ import(ggplot2)
 
 importFrom(bio3d, read.pdb, pdbseq)
 importFrom(rlog, log_info)
-importFrom(dplyr, intersect, "%>%", "reframe", "group_by")
+importFrom(dplyr, intersect, "%>%", "reframe", "group_by", "pull", "mutate")
 importFrom(comprehenr, to_list)
 importFrom(RColorBrewer, brewer.pal)
 importFrom(scales, col_bin)
@@ -31,6 +31,8 @@ importFrom("vegan", "procrustes")
 importFrom("plyr", "ldply")
 importFrom("ropls", "opls")
 importFrom("tidyr", "pivot_wider", "pivot_longer")
+importFrom("stringr", "str_count", "str_sub")
+importFrom("utils", "setTxtProgressBar", "txtProgressBar")
 
 exportClasses("RexChains",
               "RexSummary",

diff --git a/NEWS.md b/NEWS.md
@@ -1,5 +1,10 @@
 # RexMS 0.99.3
 
+* fix max uptake bug
+* improved intialiser
+
+# RexMS 0.99.3
+
 * Generate all tutorial and website
 * fix spelling and grammer issues in website
 

diff --git a/R/ReX-utils.R b/R/ReX-utils.R
@@ -44,10 +44,10 @@ prepareIndexes <- function(res) {
 maxUptakes <- function(res) {
   stopifnot("res must be a DataFrame" = is(res, "DFrame"))
 
-    .out <- vapply(seq_along(unique(res$Sequence)),
-        function(z) res$MaxUptake[res$Sequence == unique(res$Sequence)[z]][1],
-        FUN.VALUE = numeric(1)
-    )
+  .out <-  data.frame(res) %>%  
+    mutate(MaxUptake = str_count(str_sub(Sequence, 3, -1)) - 
+             str_count(str_sub(Sequence, 3, -1), "P")) %>% 
+    pull(MaxUptake)
     return(.out)
 }
 
@@ -294,6 +294,7 @@ plotUptake <- function(Uptake,
 ##'  in the model parameters and the error term (e.g. observation level error)
 ##' @param whichChains An integer value indicating which chain to sample from
 ##' @param tCoef A numeric vector of coefficients to use in the prediction
+##' @param range A numeric vector of the range of iterations to sample from
 ##'   
 ##' @return Returns a list of samples 
 ##' @examples
@@ -320,7 +321,8 @@ plotUptake <- function(Uptake,
 ##'                           thin = 1,
 ##'                           range = 1:4,
 ##'                           whichChains = c(1))
-##' samples <- marginalEffect(rex_example)                           
+##' samples <- marginalEffect(rex_example,
+##'                           range= 1:4)                           
 ##' @export
 marginalEffect <- function(params,
                            method = "fitted",
@@ -438,7 +440,7 @@ marginalEffect <- function(params,
 ##'                           thin = 1,
 ##'                           range = 1:4,
 ##'                           whichChains = c(1))
-##' samples <- marginalEffect(rex_example)
+##' samples <- marginalEffect(rex_example, range = 1:4)
 ##' plotUptakeUncertainty(samples)                           
 ##' @export
 plotUptakeUncertainty <- function(samples,

diff --git a/R/Rex-helpers.R b/R/Rex-helpers.R
@@ -11,11 +11,11 @@
 ##'
 ##' @md
 initialiser <- function(uptake_guess,
-    timepoints,
-    numRep,
-    numtimepoints,
-    R,
-    phi) {
+                        timepoints,
+                        numRep,
+                        numtimepoints,
+                        R,
+                        phi) {
     params <- matrix(NA, ncol = 4, nrow = R)
     Err <- vector(mode = "numeric", length = R)
     for (j in seq_len(R)) {
@@ -37,6 +37,23 @@ initialiser <- function(uptake_guess,
                     ),
                     silent = TRUE
                 )
+
+                if (inherits(out, "try-error")) {
+                    out <- try(
+                        minpack.lm::nlsLM(
+                            data = as.data.frame(data),
+                            formula = y ~ phi * ((1 - exp(-b * t^p))),
+                            start = list(b = 0.1, p = 0.1),
+                            lower = c(0, 0),
+                            upper = c(2, 0.99),
+                            algorithm = "LM",
+                            control = nls.control(maxiter = 1000)
+                        ),
+                        silent = TRUE
+                    )
+                }
+
+
             } else {
                 # reduce initialising model if not many data points
                 out <- try(
@@ -66,8 +83,15 @@ initialiser <- function(uptake_guess,
                 Err[j] <- sum((predict(out) - data[, 1])^2)
 
                 if ((numtimepoints > 4)) {
+                  if (length(coef(out)) == 4) {
                     params[j, ] <- coef(out)
-                } else {
+                  } else {
+                    params[j, 1:2] <- coef(out)
+                    params[j, 3] <- 1
+                    params[j, 4] <- 0.001
+                  }
+
+                } else{
                     params[j, 1:2] <- coef(out)
                     params[j, 3] <- 1
                     params[j, 4] <- 0.001

diff --git a/R/Rex-sampler.R b/R/Rex-sampler.R
@@ -116,6 +116,17 @@ resolver <- function(res,
         R = R,
         phi = phi
     )
+
+    # want to set to init_param to b if the initialiser used reduced model
+    condition_for_init_param <- sum(params[seq.int(R_lower, R_upper), 3] == 1,
+                                    na.rm = TRUE)
+    number_of_unique_b <- length(unique(params[seq.int(R_lower, R_upper), 1]))
+
+    if (condition_for_init_param > number_of_unique_b) {
+        init_param <- "b"
+        message("Initial parameter set to b")
+    }
+
 
     currentblong <- params[, 1]
     currentqlong <- params[, 2]
@@ -186,8 +197,16 @@ resolver <- function(res,
 
 
     # sampler here
+    pb <- txtProgressBar(min = 0, max = numIter, style = 3)
+    ._t <- 0
 
     for (j in seq.int(numIter)) {
+
+      if((j %% 100) ==0){
+        setTxtProgressBar(pb, ._t)
+        ._t <- ._t + 1
+      }
+
         currentSigma <- metropolisSigma(
             res = res,
             currentSigma = currentSigma,

diff --git a/R/rex-plotting.R b/R/rex-plotting.R
@@ -496,13 +496,19 @@ plotPeptideError <- function(rex_params,
   rownames(df) <- Rex.peptideError(rex_params)[,1] 
 
   if (relative) {
-    df <- df / maxUptakes(res = HdxData)
+
+    Sequence <- Rex.peptideError(rex_params)[, "Peptide"]
+    norm <- str_count(str_sub(Sequence, 3, -1)) - str_count(str_sub(Sequence, 3, -1), "P")
+    df <- df / norm
   }
 
+  rg <- max(abs(df))
+
   pp <- pheatmap(df, 
                  cluster_rows = FALSE,
                  cluster_cols = FALSE,
-                 show_rownames = TRUE)
+                 show_rownames = TRUE,
+                 breaks = seq(-rg, rg, length.out = 100))
 
 
   return(pp)

diff --git a/man/marginalEffect.Rd b/man/marginalEffect.Rd
diff --git a/man/plotUptakeUncertainty.Rd b/man/plotUptakeUncertainty.Rd
diff --git a/vignettes/ReX.Rmd b/vignettes/ReX.Rmd
@@ -520,7 +520,8 @@ has a number of arguments but the most important are the following:
 ```{r,}
 samples <- marginalEffect(rex_test,
                           method = "predict",
-                          tCoef = c(0, rep(1, numTimepoints - 1)))
+                          tCoef = c(0, rep(1, numTimepoints - 1)),
+                          range = 50:100)
 
 plotUptakeUncertainty(samples)