tests/testthat/test_150normalize_batch.R

In elsayed-lab/hpgltools: A pile of (hopefully) useful R functions

start <- as.POSIXlt(Sys.time())
context("150normalize_batch.R")

pombe_expt <- make_pombe_expt(annotation = FALSE)

testing <- normalize_expt(pombe_expt, filter = TRUE, batch = "limma")
test_counts <- exprs(testing)
actual <- as.numeric(test_counts[1:10, 1])
expected <- c(27.21933, 33.08824, 147.85699, 19.53751, 89.46974,
              206.14546, 39.54991, 117.83365, 139.79517, 20.77142)
test_that("limma batch modification provides expected values?", {
  expect_equal(expected, actual, tolerance = 0.0001)
})

testing <-  normalize_expt(pombe_expt, filter = TRUE, batch = "svaseq", surrogates = 1)
test_counts <- exprs(testing)
actual <- as.numeric(test_counts[1:10, 1])
## 20181210 These values once again changed.
##expected <- c(18.09460, 18.17614, 91.76242, 12.54344, 65.79543,
##              136.00556, 22.60825, 68.96622, 95.12627, 10.06193)
## 20190304 These values once again changed, though minimally this time.
## Probably by a small enough margin that I could just adjust the tolerance.
## But I want some record of what is going on with these...
expected <- c(18.08580, 18.18667, 91.93633, 12.54363, 65.98463,
              136.14987, 22.61450, 68.84573, 95.34932, 10.08210)
test_that("svaseq batch modification provides expected values?", {
  expect_equal(expected, actual, tolerance = 0.0001)
})

testing <-  normalize_expt(pombe_expt, filter = TRUE, batch = "fsva", surrogates = 1)
test_counts <- exprs(testing)
actual <- as.numeric(test_counts[1:10, 1])
##expected <- c(17.72009, 17.93641, 89.86064, 12.38509, 64.24830,
##              133.14076, 22.31810, 68.05949, 93.19210, 10.02845)
## These also changed very slightly.
expected <- c(17.69286, 17.92328, 89.89976, 12.37286, 64.33528,
              133.08574, 22.29550, 67.87838, 93.26926, 10.03672)
test_that("fsva batch modification provides expected values?", {
  expect_equal(expected, actual, tolerance = 0.0001)
})

testing <-  sm(normalize_expt(pombe_expt, filter = TRUE, batch = "ssva"))
test_counts <- exprs(testing)
actual <- as.numeric(test_counts[1:10, 1])
##expected <- c(25.082365, 26.421625, 87.254755, 34.719536, 60.721234,
##              124.360685, 48.678028, 72.210289, 67.940531, 9.452095)
## Same 20190304
expected <- c(25.076090, 26.506346, 87.250748, 34.717952, 60.786558,
              124.520088, 48.690709, 72.372897, 67.975395, 9.450728)
test_that("ssva batch modification provides expected values?", {
  expect_equal(expected, actual, tolerance = 0.0001)
})

testing <-  normalize_expt(pombe_expt, filter = TRUE, batch = "ruvg", surrogates = 1, thresh = 1)
test_counts <- exprs(testing)
actual <- as.numeric(test_counts[1:10, 1])
##expected <- c(9.786953, 17.456471, 17.656921, 21.515695, 86.856397,
##              13.031170, 61.307727, 124.076674, 23.202769, 63.794729)
## 20190304
expected <- c(9.797174, 17.451404, 17.659941, 21.550785, 86.959849,
              13.008831, 61.374017, 124.203114, 23.175993, 63.807827)
test_that("ruvg batch modification provides expected values?", {
  expect_equal(expected, actual, tolerance = 0.0001)
})

testing <-  normalize_expt(pombe_expt, filter = TRUE, batch = "ruv")
test_counts <- exprs(testing)
actual <- as.numeric(test_counts[1:10, 1])
##expected <- c(25.291192, 24.423036, 84.240174, 30.607522, 54.163433,
##              114.476214, 47.969103, 87.116043, 70.473976, 7.592006)
## 20190304
expected <- c(25.240082, 24.546607, 84.302026, 30.692295, 54.435697,
              115.006634, 47.900441, 87.172670, 70.529875, 7.653674)
test_that("ruv batch modification provides expected values?", {
  expect_equal(expected, actual, tolerance = 0.0001)
})

##message("FIXME! I think something is wrong with my isva invocation.")
##testing <-  normalize_expt(pombe_expt, filter = TRUE, batch = "isva")
##test_counts <- exprs(testing)
##actual <- as.numeric(test_counts[1:10, 1])
##expected <- c(-6.171629, 15.531229, -13.421936, 7.192514, -20.826294,
##              -78.238150, 77.874809, -88.632286, 3.597047, -11.672935)
##test_that("isva batch modification provides expected values?", {
##  expect_equal(expected, actual, tolerance = 0.0001)
##})

testing <-  normalize_expt(pombe_expt, filter = TRUE, batch = "combat")
test_counts <- exprs(testing)
actual <- as.numeric(test_counts[1:10, 1])
expected <- c(20.72062, 32.99341, 150.40940, 12.93767, 87.57504,
              177.59722, 42.29628, 99.52522, 146.12594, 20.49443)
test_that("combat batch modification provides expected values?", {
  ## expect_equal(expected, actual, tolerance = 0.0001)
  expect_equal(expected, actual, tolerance = 0.1)
})

## Now I think on it, compare_surrogate_estimates runs all of the others for me...
pombe_result <- compare_surrogate_estimates(pombe_expt, do_catplots = TRUE)
## Hmm I am not sure what to test in this.
expected <- 3
actual <- ncol(pombe_result[["adjustments"]][["pca_adjust"]])
test_that("Do we get expected results from compare_surrogate_estimates()?", {
  expect_equal(actual, expected)
})

pombe_filt <- normalize_expt(pombe_expt, filter = TRUE)
adjust_test_sva <- sm(all_adjusters(pombe_filt, estimate_type = "ssva"))
adjust_test_svaseq <- sm(all_adjusters(pombe_filt, estimate_type = "svaseq"))
comparison <- cor(adjust_test_sva[["new_counts"]][, 1],
                  adjust_test_svaseq[["new_counts"]][, 1])
test_that("Do we get similar sva/svaseq results?", {
  expect_gt(comparison, 0.99)
})

adjust_test_smartsva <- all_adjusters(pombe_filt, estimate_type = "smartsva")
comparison <- cor(adjust_test_smartsva[["new_counts"]][, 1],
                  adjust_test_sva[["new_counts"]][, 1])
test_that("Do we get similar sva/isva results?", {
  expect_gt(comparison, 0.99)
})

adjust_test_ruv <- all_adjusters(pombe_filt, estimate_type = "ruv")
comparison <- cor(adjust_test_ruv[["new_counts"]][, 1],
                  adjust_test_sva[["new_counts"]][, 1])
test_that("Do we get similar sva/isva results?", {
  expect_gt(comparison, 0.99)
})

## If I invoke RUVg, it returns surrogates + counts.
## I attempted to implement their solver, I therefore want to test that
## a matrix acquired from their solver + their model gets similar/identical
## counts to that acquired entirely from RUVg.  If this is true, then I think
## that gives me some confidence that I can use the solver for other methods.
ruv_estimate <- all_adjusters(pombe_filt, estimate_type = "ruvg")
ruv_counts <- ruv_estimate[["source_counts"]]
ruv_model <- ruv_estimate[["model_adjust"]]
my_counts <- counts_from_surrogates(pombe_filt, adjust = ruv_model, method = "ruv")
## 12
test_that("Do we get to a similar end point from ruv with/without solving?", {
  expect_equal(ruv_counts, my_counts, tolerance = 0.5)
})

prev_counts <- counts_from_surrogates(pombe_filt, adjust = ruv_model, method = "cbcb_subtract")


end <- as.POSIXlt(Sys.time())
elapsed <- round(x = as.numeric(end - start))
message("\nFinished 120model_surrogates.R in ", elapsed,  " seconds.")