GRN without repressions

[koenvandenberge]

Dear all,

I was hoping to simulate from a GRN without repressing interactions. I have used the following backbone and code to simulate:

backbone <- bblego(
    bblego_start("A", type = "simple"),
    bblego_linear("A", "B", type = "simple", num_modules = 5),
    bblego_linear("B", "C", type = "simple", num_modules = 5),
    bblego_end("C", type = "simple")
  )

model <- initialise_model(
  backbone = backbone,
  num_tfs = 20,
  num_targets = 500,
  num_hks = 10,
  tf_network_params = tf_network_default(sample_num_regulators = function() 2),
  feature_network_params = feature_network_default(max_in_degree=20),
  download_cache_dir = "~/.cache/dyngen",
  num_cores = 8,
  simulation_params = simulation_default(
    census_interval = 1,
    compute_cellwise_grn = TRUE
  )
)
dataset <- model %>%
  generate_tf_network() %>%
  generate_feature_network() %>%
  generate_kinetics() %>%
  generate_gold_standard() %>%
  generate_cells() %>%
  generate_experiment()

where each building block is of type simple which according to the documentation should only consider upregulations. However, when I take a look at the regulatory network I get a significant number of negative interactions:

> table(dataset$regulatory_network$effect)
 -1   1 
228 678 

Am I missing something?

Thanks! Koen

[rcannood]

Hello Koen!

The backbone indeed does not contain any repressions in your example, but the interactions between the genes from the backbone to randomly sampled targets do. You can see this in effect by running the functions separately:

library(tidyverse)
library(dyngen)

backbone <- bblego(
  bblego_start("A", type = "simple"),
  bblego_linear("A", "B", type = "simple", num_modules = 5),
  bblego_linear("B", "C", type = "simple", num_modules = 5),
  bblego_end("C", type = "simple")
)

model <- initialise_model(
  backbone = backbone,
  download_cache_dir = "~/.cache/dyngen/"
) %>%
  generate_tf_network()

model$feature_network$effect %>% table(effect = .)
# effect
#  1 
# 19 

model <- model %>% 
  generate_feature_network() %>% 
  generate_kinetics()

model$feature_network$effect %>% table(effect = .)
# effect
#  -1   1 
#  45 143 

You can override this behaviour by changing the kinetics function. Take a look at the function dyngen::kinetics_default():

function() {
  # satisfy r cmd check
  transcription_rate <- translation_rate <- mrna_halflife <- protein_halflife <- 
    independence <- splicing_rate <- effect <- strength <- hill <- NULL

  sampler_tfs <-
    function(feature_info, feature_network, cache_dir = NULL, verbose = FALSE) {
      feature_info %>% mutate(
        transcription_rate = transcription_rate %|% runif(n(), 10, 20),
        translation_rate = translation_rate %|% runif(n(), 100, 150),
        mrna_halflife = mrna_halflife %|% runif(n(), 2.5, 5),
        protein_halflife = protein_halflife %|% runif(n(), 5, 10),
        independence = independence %|% 1,
        splicing_rate = splicing_rate %|% (log(2) / 2)
      )
    }

  sampler_interactions <-
    function(feature_info, feature_network, cache_dir = NULL, verbose = FALSE) {
      feature_network %>% 
        mutate(
          effect = effect %|% sample(c(-1L, 1L), n(), replace = TRUE, prob = c(.25, .75)),
          strength = strength %|% 10 ^ runif(n(), log10(1), log10(100)),
          hill = hill %|% rnorm_bounded(n(), 2, 2, min = 1, max = 10)
        )
    }

  lst(sampler_tfs, sampler_nontfs = sampler_tfs, sampler_interactions)
}
<bytecode: 0x5612cd9fb620>
<environment: namespace:dyngen>

There are two ways of overriding the kinetics_default logic. One is to copy paste the whole function and changing the parts you would like to see changed. Another is to make a wrapper function which will slightly modify the kinetics_default() function:

my_kinetics <- kinetics_default()
my_kinetics$sampler_interactions <- 
  carrier::crate( # use carrier::crate to avoid accidental infinite recursion
    function(...) {
      df <- parent_fun(...)
      df$effect <- 1L
      df
    }, parent_fun = my_kinetics$sampler_interactions
  )

model <- initialise_model(
  backbone = backbone,
  download_cache_dir = "~/.cache/dyngen/",
  kinetics_params = my_kinetics
) %>%
  generate_tf_network() %>% 
  generate_feature_network() %>% 
  generate_kinetics()

model$feature_network$effect %>% table(effect = .)
# effect
#   1
# 183

Does this solution work for you?

Robrecht

[koenvandenberge]

Many thanks, Robrecht, for yet another clear and very useful response. This indeed provides what I was looking for.

Best, Koen