tests/testthat/test-order_cells.R

In cole-trapnell-lab/monocle3: Clustering, Differential Expression, and Trajectory Analysis for Single-Cell RNA-Seq

context("test-order_cells")

skip_not_travis <- function ()
{
  if (identical(Sys.getenv("TRAVIS"), "true")) {
    return(invisible(TRUE))
  }
  skip("Not on Travis")
}

cds <- load_a549()
set.seed(100)
test_that("order_cells error messages work", {
  skip_on_travis()
  expect_error(order_cells(cds), "No dimensionality reduction for UMAP calculated. Please run reduce_dimension with reduction_method = UMAP, cluster_cells, and learn_graph before running order_cells." )
  cds <- estimate_size_factors(cds)
  cds <- preprocess_cds(cds, num_dim = 20)
  cds <- reduce_dimension(cds)
  expect_error(order_cells(cds), "No cell clusters for UMAP calculated. Please run cluster_cells with reduction_method = UMAP and run learn_graph before running order_cells.")
  cds <- cluster_cells(cds)
  expect_error(order_cells(cds), "No principal graph for UMAP calculated. Please run learn_graph with reduction_method = UMAP before running order_cells.")
  cds <- learn_graph(cds)
  expect_error(order_cells(cds, root_cells = c("G07_B02_RT_587"), root_pr_nodes = c("Y_1")), "Please specify either root_pr_nodes or root_cells, not both.")
  expect_error(order_cells(cds, root_cells = c("hannah")), "All provided root_cells must be present in the cell data set.")
  expect_error(order_cells(cds, root_pr_nodes = c("hannah")), "All provided root_pr_nodes must be present in the principal graph.")
  expect_error(order_cells(cds), "(When not in interactive mode, either root_pr_nodes or root_cells must be provided.|No root node was chosen!)")
  expect_error(order_cells(cds, reduction_method = "tSNE"), "Currently only 'UMAP' is accepted as a reduction_method.")
})

cds <- estimate_size_factors(cds)
cds <- preprocess_cds(cds, num_dim = 20)
cds <- reduce_dimension(cds, umap.fast_sgd=FALSE)
cds <- cluster_cells(cds, cluster_method = "louvain")
cds <- learn_graph(cds)

test_that("order_cells works", {
  skip_on_travis()
  cds <- order_cells(cds, root_pr_nodes = "Y_1")
  expect_equal(max(pseudotime(cds)), 11.9, tol = 1e-1)
  expect_equal(min(pseudotime(cds)), 0)
  expect_equal(as.numeric(pseudotime(cds)[1]),  0.0538, tol = 1e-3)
  cds <- order_cells(cds, root_pr_nodes = c("Y_1", "Y_10"))
  expect_equal(max(pseudotime(cds)), 6.34, tol = 1e-2)
  expect_equal(min(pseudotime(cds)), 0)
  expect_equal(as.numeric(pseudotime(cds)[1]), 0.0538, tol = 1e-3)
  cds <- order_cells(cds, root_cells = "G07_B02_RT_587")
  expect_equal(max(pseudotime(cds)), 13.2, tol = 1e-1)
  expect_equal(min(pseudotime(cds)), 0)
  expect_equal(as.numeric(pseudotime(cds)[1]), 1.5, tol = 1e-1)
  cds <- order_cells(cds, root_cells = c("G07_B02_RT_587", "F06_A01_RT_598"))
  expect_equal(max(pseudotime(cds)), 7.26, tol = 1e-2)
  expect_equal(min(pseudotime(cds)), 0)
  expect_equal(as.numeric(pseudotime(cds)[1]), 1.5, tol = 1e-1)
})

cds <- reduce_dimension(cds, max_components = 3, umap.fast_sgd=FALSE)
cds <- cluster_cells(cds, cluster_method = "louvain")
cds <- learn_graph(cds)

test_that("order_cells works 3d", {
  skip_on_travis()
  cds <- order_cells(cds, root_pr_nodes = "Y_1")
  expect_equal(max(pseudotime(cds)), 10.0, tol = 1e-1)
  expect_equal(min(pseudotime(cds)), 0)
  expect_equal(as.numeric(pseudotime(cds)[1]),  0.664, tol = 1e-3)
  cds <- order_cells(cds, root_pr_nodes = c("Y_1", "Y_10"))
  expect_equal(max(pseudotime(cds)),  8.64, tol = 1e-2)
  expect_equal(min(pseudotime(cds)), 0)
  expect_equal(as.numeric(pseudotime(cds)[1]),  0.664, tol = 1e-3)
  cds <- order_cells(cds, root_cells = "G07_B02_RT_587")
  expect_equal(max(pseudotime(cds)), 10.4, tol = 1e-1)
  expect_equal(min(pseudotime(cds)), 0)
  expect_equal(as.numeric(pseudotime(cds)[1]), 0.664, tol = 1e-3)
  cds <- order_cells(cds, root_cells = c("G07_B02_RT_587", "F06_A01_RT_598"))
  expect_equal(max(pseudotime(cds)), 10.2, tol = 1e-1)
  expect_equal(min(pseudotime(cds)), 0)
  expect_equal(as.numeric(pseudotime(cds)[1]), 0.664, tol = 1e-3)
})

cds <- cluster_cells(cds, random_seed = 100)
cds <- learn_graph(cds)

test_that("order_cells works leiden", {
  skip_on_travis()
  cds <- order_cells(cds, root_pr_nodes = "Y_1")
  expect_equal(max(pseudotime(cds)), 9.94, tol = 1e-2)
  expect_equal(min(pseudotime(cds)), 0)
  expect_equal(as.numeric(pseudotime(cds)[1]), 4.49, tol = 1e-2)
  cds <- order_cells(cds, root_pr_nodes = c("Y_1", "Y_2"))
  expect_equal(max(pseudotime(cds)), 4.72, tol = 1e-2)
  expect_equal(min(pseudotime(cds)), 0)
  expect_equal(as.numeric(pseudotime(cds)[1]), 4.35, tol = 1e-2)
  cds <- order_cells(cds, root_cells = "G07_B02_RT_587")
  expect_equal(max(pseudotime(cds)), 8.03, tol = 1e-2)
  expect_equal(min(pseudotime(cds)), 0)
  expect_equal(as.numeric(pseudotime(cds)[1]), 0.121 , tol = 1e-3)
  cds <- order_cells(cds, root_cells = c("G07_B02_RT_587", "F06_A01_RT_598"))
  expect_equal(max(pseudotime(cds)), 5.85, tol = 1e-2)
  expect_equal(min(pseudotime(cds)), 0)
  expect_equal(as.numeric(pseudotime(cds)[1]), 0.121 , tol = 1e-3)
})

cds <- reduce_dimension(cds, max_components = 3, umap.fast_sgd=FALSE)
cds <- cluster_cells(cds)
cds <- learn_graph(cds)

test_that("order_cells works leiden 3d", {
  skip_on_travis()
  cds <- order_cells(cds, root_pr_nodes = "Y_1")
  expect_equal(max(pseudotime(cds)), 9.94, tol = 1e-2)
  expect_equal(min(pseudotime(cds)), 0)
  expect_equal(as.numeric(pseudotime(cds)[1]),  4.49, tol = 1e-2)
  cds <- order_cells(cds, root_pr_nodes = c("Y_1", "Y_2"))
  expect_equal(max(pseudotime(cds)), 4.72, tol = 1e-2)
  expect_equal(min(pseudotime(cds)), 0)
  expect_equal(as.numeric(pseudotime(cds)[1]),  4.35, tol = 1e-2)
  cds <- order_cells(cds, root_cells = "G07_B02_RT_587")
  expect_equal(max(pseudotime(cds)), 8.03, tol = 1e-2)
  expect_equal(min(pseudotime(cds)), 0)
  expect_equal(as.numeric(pseudotime(cds)[1]), 0.121, tol = 1e-3)
  cds <- order_cells(cds, root_cells = c("G07_B02_RT_587", "F06_A01_RT_598"))
  expect_equal(max(pseudotime(cds)), 5.85, tol = 1e-2)
  expect_equal(min(pseudotime(cds)), 0)
  expect_equal(as.numeric(pseudotime(cds)[1]), 0.121, tol = 1e-3)
})




#### TRAVIS ####


cds <- load_a549()
set.seed(100)
test_that("order_cells error messages work", {
  skip_not_travis()
  expect_error(order_cells(cds), "No dimensionality reduction for UMAP calculated. Please run reduce_dimension with reduction_method = UMAP, cluster_cells, and learn_graph before running order_cells." )
  cds <- estimate_size_factors(cds)
  cds <- preprocess_cds(cds, num_dim = 20)
  cds <- reduce_dimension(cds)
  expect_error(order_cells(cds), "No cell clusters for UMAP calculated. Please run cluster_cells with reduction_method = UMAP and run learn_graph before running order_cells.")
  cds <- cluster_cells(cds)
  expect_error(order_cells(cds), "No principal graph for UMAP calculated. Please run learn_graph with reduction_method = UMAP before running order_cells.")
  cds <- learn_graph(cds)
  expect_error(order_cells(cds, root_cells = c("G07_B02_RT_587"), root_pr_nodes = c("Y_1")), "Please specify either root_pr_nodes or root_cells, not both.")
  expect_error(order_cells(cds, root_cells = c("hannah")), "All provided root_cells must be present in the cell data set.")
  expect_error(order_cells(cds, root_pr_nodes = c("hannah")), "All provided root_pr_nodes must be present in the principal graph.")
  expect_error(order_cells(cds), paste("When not in interactive mode, either",
                                       "root_pr_nodes or root_cells must be",
                                       "provided."))
  expect_error(order_cells(cds, reduction_method = "tSNE"), "Currently only 'UMAP' is accepted as a reduction_method.")
})

cds <- estimate_size_factors(cds)
cds <- preprocess_cds(cds, num_dim = 20)
cds <- reduce_dimension(cds, umap.fast_sgd=FALSE)
cds <- cluster_cells(cds, cluster_method = "louvain")
cds <- learn_graph(cds)

test_that("order_cells works", {
  skip_not_travis()
  cds <- order_cells(cds, root_pr_nodes = "Y_1")
  expect_equal(max(pseudotime(cds)), 11.9, tol = 1e-1)
  expect_equal(min(pseudotime(cds)), 0)
  expect_equal(as.numeric(pseudotime(cds)[1]), 0.0538, tol = 1e-4)
  cds <- order_cells(cds, root_pr_nodes = c("Y_1", "Y_10"))
  expect_equal(max(pseudotime(cds)), 6.34, tol = 1e-2)
  expect_equal(min(pseudotime(cds)), 0)
  expect_equal(as.numeric(pseudotime(cds)[1]), 0.0538, tol = 1e-4)
  cds <- order_cells(cds, root_cells = "G07_B02_RT_587")
  expect_equal(max(pseudotime(cds)), 13.2, tol = 1e-1)
  expect_equal(min(pseudotime(cds)), 0)
  expect_equal(as.numeric(pseudotime(cds)[1]), 1.5, tol = 1e-1)
  cds <- order_cells(cds, root_cells = c("G07_B02_RT_587", "F06_A01_RT_598"))
  expect_equal(max(pseudotime(cds)), 7.26, tol = 1e-2)
  expect_equal(min(pseudotime(cds)), 0)
  expect_equal(as.numeric(pseudotime(cds)[1]), 1.5, tol = 1e-1)
})

cds <- reduce_dimension(cds, max_components = 3, umap.fast_sgd=FALSE)
cds <- cluster_cells(cds, cluster_method = "louvain")
cds <- learn_graph(cds)

test_that("order_cells works 3d", {
  skip_not_travis()
  cds <- order_cells(cds, root_pr_nodes = "Y_1")
  expect_equal(max(pseudotime(cds)), 10.4, tol = 1e-1)
  expect_equal(min(pseudotime(cds)), 0)
  expect_equal(as.numeric(pseudotime(cds)[1]), 0.664, tol = 1e-3)
  cds <- order_cells(cds, root_pr_nodes = c("Y_1", "Y_10"))
  expect_equal(max(pseudotime(cds)),  8.64, tol = 1e-2)
  expect_equal(min(pseudotime(cds)), 0)
  expect_equal(as.numeric(pseudotime(cds)[1]), 0.664, tol = 1e-3)
  cds <- order_cells(cds, root_cells = "G07_B02_RT_587")
  expect_equal(max(pseudotime(cds)), 10.4, tol = 1e-1)
  expect_equal(min(pseudotime(cds)), 0)
  expect_equal(as.numeric(pseudotime(cds)[1]), 0.664, tol = 1e-3)
  cds <- order_cells(cds, root_cells = c("G07_B02_RT_587", "F06_A01_RT_598"))
  expect_equal(max(pseudotime(cds)), 10.4, tol = 1e-1)
  expect_equal(min(pseudotime(cds)), 0)
  expect_equal(as.numeric(pseudotime(cds)[1]), 0.664, tol = 1e-3)
})

cds <- cluster_cells(cds, random_seed = 100)
cds <- learn_graph(cds)

test_that("order_cells works leiden", {
  skip_not_travis()
  cds <- order_cells(cds, root_pr_nodes = "Y_1")
  expect_equal(max(pseudotime(cds)), 9.94, tol = 1e-2)
  expect_equal(min(pseudotime(cds)), 0)
  expect_equal(as.numeric(pseudotime(cds)[1]), 4.49, tol = 1e-2)
  cds <- order_cells(cds, root_pr_nodes = c("Y_1", "Y_2"))
  expect_equal(max(pseudotime(cds)), 4.72, tol = 1e-2)
  expect_equal(min(pseudotime(cds)), 0)
  expect_equal(as.numeric(pseudotime(cds)[1]), 4.35, tol = 1e-2)
  cds <- order_cells(cds, root_cells = "G07_B02_RT_587")
  expect_equal(max(pseudotime(cds)), 8.03, tol = 1e-2)
  expect_equal(min(pseudotime(cds)), 0)
  expect_equal(as.numeric(pseudotime(cds)[1]), 0.121, tol = 1e-3)
  cds <- order_cells(cds, root_cells = c("G07_B02_RT_587", "F06_A01_RT_598"))
  expect_equal(max(pseudotime(cds)), 6.1, tol = 1e-1)
  expect_equal(min(pseudotime(cds)), 0)
  expect_equal(as.numeric(pseudotime(cds)[1]), 0.121, tol = 1e-3)
})

cds <- reduce_dimension(cds, max_components = 3, umap.fast_sgd=FALSE)
cds <- cluster_cells(cds)
cds <- learn_graph(cds)

test_that("order_cells works leiden 3d", {
  skip_not_travis()
  cds <- order_cells(cds, root_pr_nodes = "Y_1")
  expect_equal(max(pseudotime(cds)), 9.94, tol = 1e-2)
  expect_equal(min(pseudotime(cds)), 0)
  expect_equal(as.numeric(pseudotime(cds)[1]), 4.49, tol = 1e-2)
  cds <- order_cells(cds, root_pr_nodes = c("Y_1", "Y_2"))
  expect_equal(max(pseudotime(cds)),  4.72, tol = 1e-2)
  expect_equal(min(pseudotime(cds)), 0)
  expect_equal(as.numeric(pseudotime(cds)[1]), 4.35, tol = 1e-2)
  cds <- order_cells(cds, root_cells = "G07_B02_RT_587")
  expect_equal(max(pseudotime(cds)), 8.03, tol = 1e-2)
  expect_equal(min(pseudotime(cds)), 0)
  expect_equal(as.numeric(pseudotime(cds)[1]), 0.121, tol = 1e-3)
  cds <- order_cells(cds, root_cells = c("G07_B02_RT_587", "F06_A01_RT_598"))
  expect_equal(max(pseudotime(cds)), 6.1, tol = 1e-1)
  expect_equal(min(pseudotime(cds)), 0)
  expect_equal(as.numeric(pseudotime(cds)[1]), 0.121, tol = 1e-3)
})