R/calculate_trajectory_properties.R

In duohongrui/simutils: Common functions used in data simulation and evaluation

#' Calculate Metrics for Comparison of Two Trajectories
#'
#' @param ref_data A gene matrix or the trajectory object generated by `dynwrap::wrap_expression` from real data
#' @param ref_data_grouping The labels of cells in real data
#' @param sim_data A gene matrix or the trajectory object generated by `dynwrap::wrap_expression` from simulated data
#' @param sim_data_grouping The labels of cells in simulated data
#' @param algorithm Optional. Which algorithm used for matching cells in simulated and real data. Improved_Hungarian (default), Hungarian.
#' @param seed Random seed
#' @param verbose Whether the messages are returned to users when processing
#' @importFrom dynwrap is_data_wrapper add_grouping infer_trajectory add_expression
#'
#' @return A list
#' @export
#'
calculate_trajectory_properties <- function(
  ref_data,
  ref_data_grouping = NULL,
  sim_data,
  sim_data_grouping = NULL,
  algorithm = "Hungarian",
  seed = 1,
  verbose = TRUE
){
  ### Check data
  if(is.data.frame(ref_data)){
    ref_data <- as.matrix(ref_data)
  }
  if(is.data.frame(sim_data)){
    sim_data <- as.matrix(sim_data)
  }
  ### Build standard objects
  if(!dynwrap::is_data_wrapper(ref_data)){
    ref_data <- dynwrap::wrap_expression(counts = t(ref_data),
                                         expression = log2(t(ref_data) + 1))
    ref_data_transformation <- TRUE
  }else{
    ref_data_transformation <- FALSE
  }
  if(!dynwrap::is_data_wrapper(sim_data)){
    sim_data <- dynwrap::wrap_expression(counts = t(sim_data),
                                         expression = log2(t(sim_data) + 1))
    sim_data_transformation <- TRUE
  }else{
    sim_data_transformation <- FALSE
  }
  ### grouping
  if(!is.null(ref_data_grouping)){
    ref_data <- dynwrap::add_grouping(dataset = ref_data,
                                      grouping = ref_data_grouping)
  }
  if(!is.null(sim_data_grouping)){
    sim_data <- dynwrap::add_grouping(dataset = sim_data,
                                      grouping = sim_data_grouping)
  }
  ### trajectory inference
  if(!requireNamespace("tislingshot", quietly = TRUE)){
    message("Install tislingshot...")
    devtools::install_github("dynverse/ti_slingshot/package/")
  }
  if(!dynwrap::is_wrapper_with_trajectory(ref_data)){
    message("Performing trajectory inference by Slingshot for reference data...")
    ref_model <- dynwrap::infer_trajectory(dataset = ref_data,
                                           method = tislingshot::ti_slingshot(),
                                           parameters = NULL,
                                           give_priors = NULL,
                                           seed = seed,
                                           verbose = verbose)
  }else{
    ref_model <- ref_data
  }
  if(!dynwrap::is_wrapper_with_trajectory(sim_data)){
    message("Performing trajectory inference by Slingshot for simulated data...")
    sim_model <- dynwrap::infer_trajectory(dataset = sim_data,
                                           method = tislingshot::ti_slingshot(),
                                           parameters = NULL,
                                           give_priors = NULL,
                                           seed = seed,
                                           verbose = verbose)
  }else{
    sim_model <- sim_data
  }
  ### counts and expression
  if(ref_data_transformation){
    ref_model <- dynwrap::add_expression(ref_model,
                                         counts = ref_data$counts,
                                         expression = ref_data$expression)
  }
  if(sim_data_transformation){
    sim_model <- dynwrap::add_expression(sim_model,
                                         counts = sim_data$counts,
                                         expression = sim_data$expression)
  }
  ### Calculate metrics
  result <- calculate_traj_metrics(model_ref = ref_model,
                                   model_sim = sim_model,
                                   algorithm = algorithm)
  return(result)
}