R/2_tf_network.R

In dyngen: A Multi-Modal Simulator for Spearheading Single-Cell Omics Analyses

#' Generate a transcription factor network from the backbone
#' 
#' [generate_tf_network()] generates the transcription factors (TFs) that 
#' drive the dynamic process a cell undergoes. 
#' [tf_network_default()] is used to configure parameters pertaining this process.
#' 
#' @param model A dyngen initial model created with [initialise_model()].
#' @param min_tfs_per_module The number of TFs to generate per module in the backbone.
#' @param sample_num_regulators A function to generate the number of TFs per module each TF will be regulated by.
#' @param weighted_sampling When determining what TFs another TF is regulated by, whether to perform weighted
#'  sampling (by rank) or not.
#' 
#' @export
#' @rdname generate_tf_network
#' 
#' @return A dyngen model.
#' 
#' @seealso [dyngen] on how to run a complete dyngen simulation
#' 
#' @examples
#' model <- 
#'   initialise_model(
#'     backbone = backbone_bifurcating()
#'   )
#' model <- model %>%
#'   generate_tf_network() 
#'   
#' \donttest{
#' plot_feature_network(model)
#' }
generate_tf_network <- function(
  model
) {
  if (model$verbose) cat("Generating TF network\n")
  
  model <- model %>% 
    .add_timing("2_tf_network", "generate_tf_info") %>% 
    .generate_tf_info()
  
  model <- model %>% 
    .add_timing("2_tf_network", "generate_tf_network") %>% 
    .generate_tf_network() 
  
  .add_timing(model, "2_tf_network", "end")
}

#' @export
#' @rdname generate_tf_network
tf_network_default <- function(
  min_tfs_per_module = 1L,
  sample_num_regulators = function() 2,
  weighted_sampling = FALSE
) {
  lst(
    min_tfs_per_module,
    sample_num_regulators,
    weighted_sampling
  )
}

.generate_tf_info <- function(model) {
  # satisfy r cmd check
  module_id <- num_tfs <- feature_id <- NULL
  
  module_info <- model$backbone$module_info
  module_network <- model$backbone$module_network
  numbers <- model$numbers
  
  module_info <- 
    module_info %>% mutate(
      num_tfs = .generate_partitions(
        num_elements = max(numbers$num_modules, numbers$num_tfs),
        num_groups = numbers$num_modules, 
        min_elements_per_group = model$tf_network_params$min_tfs_per_module
      ),
      feature_id = map2(module_id, num_tfs, function(module_id, num_tfs) paste0(module_id, "_TF", seq_len(num_tfs))),
      is_tf = TRUE,
      is_hk = FALSE
    )
  
  model$feature_info <- 
    module_info %>% 
    unnest(feature_id) %>% 
    select(feature_id, everything(), -num_tfs)
  
  model
}

.generate_tf_network <- function(model) {
  # satisfy r cmd check
  to <- from <- module_id <- feature_id <- NULL
  
  module_network <- model$backbone$module_network
  tf_info <- model$feature_info
  
  # initialise model structures
  tf_network <- map(tf_info$feature_id, ~list()) %>% set_names(tf_info$feature_id)
  num_targets <- rep(0, length(tf_network)) %>% set_names(tf_info$feature_id)
  
  # go over each tf to find their regulators
  for (i in sample.int(nrow(tf_info))) {
    fi <- tf_info$feature_id[[i]]
    mi <- tf_info$module_id[[i]]
    
    # what are the regulating modules, make sure that each regulating
    # module has at least one regulator
    regulating_modules <- module_network %>% filter(to == mi) %>% pull(from)
    
    # generate regulating tfs
    edges <- 
      map_df(
        regulating_modules,
        function(mreg) {
          # what tfs are in this module
          candidate_regulating_tfs <- tf_info %>% filter(module_id == mreg) %>% pull(feature_id)
          
          # how many regulators will we sample?
          num_regulating_tfs <-
            model$tf_network_params$sample_num_regulators() %>% 
            min(length(candidate_regulating_tfs)) %>%
            max(1L)
          
          weights <- 
            if (model$tf_network_params$weighted_sampling) {
              # do weighted sampling based on the number of
              # targes the candidate regulator already has
              num_targets[candidate_regulating_tfs] + 1
            } else {
              NULL
            }
          regulating_tfs <- sample(candidate_regulating_tfs, num_regulating_tfs, replace = FALSE, prob = weights)
          
          tibble(
            from = regulating_tfs,
            to = fi,
            from_module = mreg,
            to_module = mi
          )
        }
      )
    
    if (nrow(edges) > 0) {
      num_targets[edges$from] <- num_targets[edges$from] + 1
      tf_network[[fi]] <- edges
    }
  }
  
  model$feature_network <- 
    bind_rows(unname(tf_network)) %>% 
    left_join(module_network %>% rename(from_module = from, to_module = to), by = c("from_module", "to_module"))
  
  model
}