Figure_5

#The source code to reproduce Figure 5.
library(tidyverse)
library(furrr)
library(dyno)
library(dyngen)
library(SingleCellExperiment)
library(spectral)
library(Lamian)
library(tradeSeq)
library(tictoc)
theme_set(theme_bw(base_size = 20))

num.targets <- 250
num.hks <- 250
num.features <- 100
de.rate <- 0.5
ko.rate <- 0
num.cores <- 1
cell.range <- c(100, 500, 1000, 5000)

#parameters and functions
min.freq <- 0.2
max.freq <- 1
length.freq <- 101
sim.num <- 100
ko_gene <- "C1_TF1"
LS <- function(g, time1, exp1, time2 = NULL, exp2 = NULL, min.freq = 0, max.freq = 1, length.freq = 101, m = "generalized", sim.num = 100){
  .quiet <- function(x){ 
    sink(tempfile()) 
    on.exit(sink()) 
    invisible(force(x)) 
  }
  time1 <- time1[!is.infinite(time1)]
  exp1 <- exp1[,!is.infinite(time1)]
  if(is.null(time2)){
    seq1 <- exp1[which(rownames(exp1) == g),]
    ls1 <- spec.lomb(x = time1, y = seq1,
                     f = seq(min.freq, max.freq, length = length.freq), mode = m) %>% .quiet
    return(tibble(gene = g, p.dynamic = min(ls1$p), p.de = NA))
  }else{
    time2 <- time2[!is.infinite(time2)]
    exp2 <- exp2[,!is.infinite(time2)]
    seq1 <- exp1[which(rownames(exp1) == g),]
    seq2 <- exp2[which(rownames(exp2) == g),]
    ls1 <- spec.lomb(x = time1, y = seq1,
                     f = seq(min.freq, max.freq, length = length.freq), mode = m) %>% .quiet
    ls2 <- spec.lomb(x = time2, y = seq2,
                     f = seq(min.freq, max.freq, length = length.freq), mode = m) %>% .quiet
    
    null.hypothesis <- NULL
    for(i in seq(sim.num)){
      sim.ls1 <- spec.lomb(x = sample(time1, length(time1)),
                           y = seq1,
                           f = seq(min.freq, max.freq, length = length.freq), mode = m) %>% .quiet
      sim.ls2 <- spec.lomb(x = sample(time2, length(time2)),
                           y = seq2,
                           f = seq(min.freq, max.freq, length = length.freq), mode = m) %>% .quiet
      null.hypothesis <- c(null.hypothesis, dist(rbind(sim.ls1$A, sim.ls2$A), method = "canberra"))
    }
    p <- sum(dist(rbind(ls1$A, ls2$A), method = "canberra") <  null.hypothesis) / sim.num
    return(tibble(gene = g, p.dynamic = min(c(ls1$p, ls2$p)), p.de = p))
  }
}

res <- tibble(method = NA,
              cell.num = NA,
              test = NA,
              time = NA)

model <- "bifurcating"
for(num.cells in cell.range){
  set.seed(8)
  #model generation
  dataset <- dyntoy::generate_dataset(model = model,
                                      num_cells = num.cells,
                                      num_features = num.features,
                                      differentially_expressed_rate = de.rate) %>%
    add_root(root_milestone_id = .$prior_information$start_milestones) %>%
    add_pseudotime()
  save(dataset, file = str_c("./models/dynamicTest_", num.cells, "cells.RData"))

  backbone <- backbone_bifurcating()
  config <- initialise_model(backbone = backbone,
                             num_cells = num.cells,
                             num_tfs = nrow(backbone$module_info),
                             num_targets = num.targets,
                             num_hks = num.hks,
                             num_cores = num.cores,
                             simulation_params = simulation_default(census_interval = 10,
                                                                    ssa_algorithm = ssa_etl(tau = 300 / 3600),
                                                                    experiment_params = simulation_type_wild_type(num_simulations = 100)))
  model_common <- config %>%
    generate_tf_network() %>%
    generate_feature_network() %>% 
    generate_kinetics() %>%
    generate_gold_standard()
  model_wt <- model_common %>% generate_cells()
  model_ko <- model_common
  model_ko$simulation_params$experiment_params <- simulation_type_knockdown(num_simulations = 100,
                                                                            timepoint = 0, 
                                                                            genes = ko_gene,
                                                                            num_genes = 1,
                                                                            multiplier = ko.rate)
  model_ko <- model_ko %>% generate_cells()
  wt <- model_wt %>% generate_experiment() %>% as_dyno() %>% add_root(root_milestone_id = "sA") %>% add_pseudotime(pseudotime = NULL)
  ko <- model_ko %>% generate_experiment() %>% as_dyno() %>% add_root(root_milestone_id = "sA") %>% add_pseudotime(pseudotime = NULL)
  save(wt, file = str_c("./models/DEtest_", num.cells, "_wt.RData"))
  save(ko, file = str_c("./models/DEtest_", num.cells, "_ko.RData"))
}

for(num.cells in cell.range){
  #dynamic test
  load(str_c("./models/dynamicTest_", num.cells, "cells.RData"))
  counts <- t(as.matrix(dataset$counts))
  expression <- t(as.matrix(dataset$expression))
  sce <- SingleCellExperiment(assays = List(counts = counts, log2 = log2(counts)))
  pseudotime <- dataset$pseudotime
  
  #LS
  .pseudotime <- pseudotime / max(pseudotime)
  plan(sequential)
  tic()
  ls.res <- future_map_dfr(rownames(expression), LS, .pseudotime, expression,
                           .options = furrr_options(seed = 8))
  time <- toc()
  res <- res %>% add_row(method = "LS",
                         cell.num = num.cells,
                         test = "Dynamic test",
                         time = time$callback_msg %>% str_split(pattern = " ", simplify = TRUE) %>% .[,1] %>% as.numeric)
  
  #tradeSeq
  cellWeights <- rep(1, num.cells)
  tic()
  sce <- fitGAM(counts = counts, pseudotime = pseudotime, cellWeights = cellWeights,
                nknots = 6, verbose = FALSE)
  assoRes <- associationTest(sce)
  time <- toc()
  res <- res %>% add_row(method = "tradeSeq (association)",
                         cell.num = num.cells,
                         test = "Dynamic test",
                         time = time$callback_msg %>% str_split(pattern = " ", simplify = TRUE) %>% .[,1] %>% as.numeric)

  tic()
  sce <- fitGAM(counts = counts, pseudotime = pseudotime, cellWeights = cellWeights,
                nknots = 6, verbose = FALSE)
  sveRes <- startVsEndTest(sce)
  time <- toc()
  res <- res %>% add_row(method = "tradeSeq (startVsEnd)",
                         cell.num = num.cells,
                         test = "Dynamic test",
                         time = time$callback_msg %>% str_split(pattern = " ", simplify = TRUE) %>% .[,1] %>% as.numeric)
    
  #Lamian
  if(num.cells %% 2 == 0){
    .cellanno <- data.frame(Cell = colnames(expression), Sample = rep(c("WT1", "WT2"), each = num.cells / 2))
  }else{
    .cellanno <- data.frame(Cell = colnames(expression), Sample = rep(c("WT1", "WT2"), c(round(num.cells / 2), round(num.cells / 2) - 1)))
  }
  
  .design <- matrix(c(1, 1, 0, 1), ncol = 2, dimnames = list(c("WT1", "WT2"), c("intercept", "group")))
  lamian.res <- NULL
  tic()
  try(lamian.res <- lamian_test(expr = expression,
                                cellanno = .cellanno,
                                pseudotime = pseudotime,
                                design = .design,
                                test.type = 'time',
                                permuiter = 100,
                                ncores = num.cores))
  time <- toc()
  
  if(!is.null(lamian.res)){
    res <- res %>% add_row(method = "Lamian",
                           cell.num = num.cells,
                           test = "Dynamic test",
                           time = time$callback_msg %>% str_split(pattern = " ", simplify = TRUE) %>% .[,1] %>% as.numeric)
  }else{
    res <- res %>% add_row(method = "Lamian",
                           cell.num = num.cells,
                           test = "Dynamic test",
                           time = NA)
  }
  
  
  #DE test
  load(str_c("./models/DEtest_", num.cells, "_wt.RData"))
  load(str_c("./models/DEtest_", num.cells, "_ko.RData"))
  wt.counts <- wt$counts %>% as.matrix() %>% t()
  ko.counts <- ko$counts %>% as.matrix() %>% t()
  wt.expression <- wt$expression %>% as.matrix() %>% t()
  wt.counts <- wt$counts %>% as.matrix() %>% t()
  wt.pseudotime <- wt$pseudotime
  ko.expression <- ko$expression %>% as.matrix() %>% t()
  ko.counts <- ko$counts %>% as.matrix() %>% t()
  ko.pseudotime <- ko$pseudotime
  counts <- cbind(wt.counts, ko.counts)
  
  #LS
  plan(sequential)
  tic()
  res.LS <- future_map_dfr(rownames(wt.expression), LS, wt.pseudotime, wt.expression,
                           ko.pseudotime, ko.expression, .options = furrr_options(seed = 8))
  time <- toc()
  res <- res %>% add_row(method = "LS",
                         cell.num = num.cells,
                         test = "DE test",
                         time = time$callback_msg %>% str_split(pattern = " ", simplify = TRUE) %>% .[,1] %>% as.numeric)
  
  #Lamian
  .expression <- cbind(wt.expression, ko.expression)
  colnames(.expression) <- c(str_c("WT_", colnames(wt.expression)), str_c("KO_", colnames(ko.expression)))
  .pseudotime <- c(wt.pseudotime, ko.pseudotime)
  names(.pseudotime) <- c(str_c("WT_", colnames(wt.expression)), str_c("KO_", colnames(ko.expression)))
  .cellanno <- data.frame(Cell = colnames(.expression),
                          Sample = rep(c("WT", "KO"), each = num.cells / 2))
  .design <- matrix(c(1, 1, 0, 1), ncol = 2, dimnames = list(c("WT", "KO"), c("intercept", "group")))
  tic()
  res.lamian <- lamian_test(expr = .expression,
                            cellanno = .cellanno,
                            pseudotime = .pseudotime,
                            design = .design,
                            test.type = 'variable',
                            overall.only = TRUE,
                            test.method = "chisq",
                            permuiter = 100,
                            ncores = num.cores)
  time <- toc()
  res <- res %>% add_row(method = "Lamian",
                         cell.num = num.cells,
                         test = "DE test",
                         time = time$callback_msg %>% str_split(pattern = " ", simplify = TRUE) %>% .[,1] %>% as.numeric)
  
  
  #tradeSeq
  colnames(counts) <- colnames(.expression)
  .condition <- factor(rep(c("WT", "KO"), each = num.cells), levels = c("WT", "KO"))
  .pseudotime <-  matrix(c(wt.pseudotime, rep(0, num.cells), rep(0, num.cells), ko.pseudotime), ncol = 2)
  .cellweights <- matrix(c(rep(1, num.cells), rep(0, num.cells), rep(0, num.cells), rep(1, num.cells)), ncol = 2)
  
  tic()
  gam <- fitGAM(counts = counts, pseudotime = .pseudotime,
                cellWeights = .cellweights,
                nknots = 6, parallel = FALSE)
  res.tradeseq.diffEnd <- diffEndTest(gam)
  time <- toc()
  res <- res %>% add_row(method = "tradeSeq (diffEndTest)",
                         cell.num = num.cells,
                         test = "DE test",
                         time = time$callback_msg %>% str_split(pattern = " ", simplify = TRUE) %>% .[,1] %>% as.numeric)
  
  tic()
  gam <- fitGAM(counts = counts, pseudotime = .pseudotime,
                cellWeights = .cellweights,
                nknots = 6, parallel = FALSE)
  res.tradeseq.pattern <- patternTest(gam)
  time <- toc()
  res <- res %>% add_row(method = "tradeSeq (patternTest)",
                         cell.num = num.cells,
                         test = "DE test",
                         time = time$callback_msg %>% str_split(pattern = " ", simplify = TRUE) %>% .[,1] %>% as.numeric)
  
  tic()
  res.tradeseq.earlyDE <- earlyDETest(gam, knots = c(1, 3))
  gam <- fitGAM(counts = counts, pseudotime = .pseudotime,
                cellWeights = .cellweights,
                nknots = 6, parallel = FALSE)
  time <- toc()
  res <- res %>% add_row(method = "tradeSeq (earlyDETest)",
                         cell.num = num.cells,
                         test = "DE test",
                         time = time$callback_msg %>% str_split(pattern = " ", simplify = TRUE) %>% .[,1] %>% as.numeric)
  
}
res$test <- factor(res$test, levels = c("Dynamic test", "DE test"))
res <- res[-1,]
save(res, file = "./res.RData")
res$method <- factor(res$method, levels = c("LS", "tradeSeq (startVsEnd)", "tradeSeq (association)", "tradeSeq (earlyDETest)", "tradeSeq (diffEndTest)", "tradeSeq (patternTest)", "Lamian"))

g <- ggplot(res, aes(x = cell.num %>% as.factor, y = time, fill = method)) +
  geom_bar(stat = "identity", position = "dodge") +
  xlab("cell number") + ylab("Time (sec)") +
  facet_wrap(. ~ test, scales = "free") +
  theme(legend.title = element_blank(), legend.position = "bottom", legend.text = element_text(size = 20)) +
  geom_text(aes(label = str_c(round(time, 1), " sec"), y = time + 3), position = position_dodge(width = 0.9))
g
ggsave(g, file = "./time_label.pdf", dpi = 300)