fix: use modified times for peak start and end

added test for correct_peaks
ethanbass · Dec 18, 2023 · ca1dd4a · ca1dd4a
1 parent 2536d6a
commit ca1dd4a
Show file tree

Hide file tree

Showing 3 changed files with 37 additions and 28 deletions.
diff --git a/R/correct_rt.R b/R/correct_rt.R
@@ -276,7 +276,11 @@ correct_peaks <- function(peak_list, mod_list, chrom_list){
              if (nrow(profile) > 0){
                cbind(profile,
                      rt.cor = c(predict.ptw(mod, profile[, "rt"], what = "time",
-                                                  RTref = ref_times)))
+                                                  RTref = ref_times)),
+                     start.cor = c(predict.ptw(mod, profile[, "start"], what = "time",
+                                            RTref = ref_times)),
+                     end.cor = c(predict.ptw(mod, profile[, "end"], what = "time",
+                                            RTref = ref_times)))
              } else {
                cbind(profile, rt.cor = rep(0, 0))
              }

diff --git a/R/get_peaktable.R b/R/get_peaktable.R
@@ -100,6 +100,8 @@ get_peaktable <- function(peak_list, chrom_list, response = c("area", "height"),
   clust <- match.arg(clust, c("rt","sp.rt"))
   out <- match.arg(out, c('data.frame', 'matrix'))
   rt <- ifelse(use.cor, "rt.cor", "rt")
+  start <- ifelse(use.cor, "start.cor", "start")
+  end <- ifelse(use.cor, "end.cor", "end")
   if (!inherits(peak_list, "peak_list"))
     stop("Peak list must be of the associated class.")
   if (clust == "sp.rt"){
@@ -150,7 +152,7 @@ get_peaktable <- function(peak_list, chrom_list, response = c("area", "height"),
       sing <- which(pkcenters.cl == 0)
       pkcenters.cl[sing] <- max(pkcenters.cl) + seq_along(sing)
     }
-    vars <- c(rt, "start", "end", "sd", "width", "tau", "FWHM", "r.squared", "purity")
+    vars <- c(rt, start, end, "sd", "width", "tau", "FWHM", "r.squared", "purity")
     vars <- vars[vars %in% colnames(xx)]
     vars.idx <- match(vars, colnames(xx))
     cl.centers <- aggregate(xx[, vars.idx], by = list(pkcenters.cl), FUN = "mean",

diff --git a/tests/testthat/test-pipeline.R b/tests/testthat/test-pipeline.R
@@ -150,32 +150,35 @@ test_that("get_peaktable works", {
   expect_equal(class(pk_tab$args), "list")
 })
 
-# dat.pr <- Sa_pr
-# test_that("correct_peaks works", {
-#   skip_on_cran()
-#   pks <- get_peaks(dat.pr, lambdas = 210)
-#   warping.models <- correct_rt(chrom_list = dat.pr,lambdas = 210, what="models")
-#   pks.cor <- correct_peaks(pks, warping.models)
-#   pktab_cor <- get_peaktable(pks.cor, use.cor = TRUE)
-#   
-#   ptw_warp <- correct_rt(dat.pr, models=warping.models)
-#   pks_warp <- get_peaks(ptw_warp, lambdas = 210, fit = "egh", smooth_type = "none",
-#                        show_progress = FALSE)
-#   pktab_warp <- get_peaktable(pks_warp)
-#   
-#   pks_reg <- get_peaks(dat.pr, lambdas = 210, show_progress = FALSE)
-#   pktab_reg <- get_peaktable(pks_reg, plot_it = TRUE)
-#   
-#   mean(colMeans(abs(pktab_cor$tab - pktab_warp$tab)))
-#   mean(colMeans(abs(pktab_cor$tab - pktab_reg$tab)))
-#   
-#   expect_equal(rownames(pk_tab$tab), names(dat.pr))
-#   expect_equal(colnames(pk_tab$tab), colnames(pk_tab$pk_meta))
-#   expect_equal(class(pk_tab), "peak_table")
-#   expect_equal(class(pk_tab$tab), "data.frame")
-#   expect_equal(class(pk_tab$pk_meta), "data.frame")
-#   expect_equal(class(pk_tab$args), "list")
-# })
+test_that("correct_peaks works", {
+  skip_on_cran()
+  data(Sa_pr)
+  pks <- get_peaks(Sa_pr, lambdas = 210, show_progress = FALSE)
+
+  warping.models <- correct_rt(chrom_list = Sa_pr, lambdas = 210, 
+                               what = "models", show_progress = FALSE)
+  pks_cor <- correct_peaks(pks, warping.models)
+  pktab_cor <- get_peaktable(pks_cor, use.cor = TRUE)
+
+  ptw_warp <- correct_rt(Sa_pr, models = warping.models)
+  pks_warp <- get_peaks(ptw_warp, lambdas = 210, show_progress = FALSE)
+  pktab_warp <- get_peaktable(pks_warp)
+
+  expect_equal(pks_cor[[1]][[1]]$rt.cor[-33], pks_warp[[1]][[1]]$rt, 
+               tolerance = .001)
+  expect_equal(pks_cor[[2]][[1]]$rt.cor, pks_warp[[2]][[1]]$rt[-6], 
+               tolerance = .001)
+
+  pks_reg <- get_peaks(Sa_pr, lambdas = 210, show_progress = FALSE)
+  pktab_reg <- get_peaktable(pks_reg)
+
+  expect_equal(row.names(pktab_cor), names(Sa_pr))
+  expect_equal(row.names(pktab_reg), names(Sa_pr))
+  expect_equal(row.names(pktab_warp), names(Sa_pr))
+
+  expect_lt(ncol(pktab_warp), ncol(pktab_reg))
+  expect_lt(ncol(pktab_cor), ncol(pktab_reg))
+})
 
 test_that("strip plot works", {
   skip_on_cran()