dynverse / dyngen

@@ -3,6 +3,8 @@
Loading
3 3
#' [generate_kinetics()] samples the kinetics of genes in the feature network for which 
4 4
#'   the kinetics have not yet been defined.
5 5
#' [kinetics_default()] is used to configure parameters pertaining this process.
6 +
#' [kinetics_random_distributions()] will do the same, but the distributions are also 
7 +
#' randomised.
6 8
#' 
7 9
#' @param model A dyngen intermediary model for which the feature network has been generated with [generate_feature_network()].
8 10
#' 
@@ -56,9 +58,6 @@
Loading
56 58
generate_kinetics <- function(model) {
57 59
  model <- .add_timing(model, "4_kinetics", "checks")
58 60
  
59 -
  # satisfy r cmd check
60 -
  burn <- mol_premrna <- mol_mrna <- mol_protein <- val <- NULL
61 -
  
62 61
  assert_that(
63 62
    !is.null(model$feature_info),
64 63
    !is.null(model$feature_network)
@@ -100,10 +99,10 @@
Loading
100 99
  # determine variables to be used during burn in
101 100
  burn_variables <- 
102 101
    model$feature_info %>% 
103 -
    filter(burn) %>% 
104 -
    select(mol_premrna, mol_mrna, mol_protein) %>% 
105 -
    gather(col, val) %>% 
106 -
    pull(val)
102 +
    filter(.data$burn) %>% 
103 +
    select(.data$mol_premrna, .data$mol_mrna, .data$mol_protein) %>% 
104 +
    gather("col", "val") %>% 
105 +
    pull(.data$val)
107 106
    
108 107
  # return system
109 108
  model <- .add_timing(model, "4_kinetics", "create output")
@@ -122,41 +121,65 @@
Loading
122 121
#' @rdname generate_kinetics
123 122
#' @importFrom stats runif
124 123
kinetics_default <- function() {
125 -
  # satisfy r cmd check
126 -
  transcription_rate <- translation_rate <- mrna_halflife <- protein_halflife <- 
127 -
    independence <- splicing_rate <- effect <- strength <- hill <- NULL
128 -
  
129 -
  sampler_tfs <-
130 -
    function(feature_info, feature_network, cache_dir = NULL, verbose = FALSE) {
131 -
      feature_info %>% mutate(
132 -
        transcription_rate = transcription_rate %|% runif(n(), 10, 20),
133 -
        translation_rate = translation_rate %|% runif(n(), 100, 150),
134 -
        mrna_halflife = mrna_halflife %|% runif(n(), 2.5, 5),
135 -
        protein_halflife = protein_halflife %|% runif(n(), 5, 10),
136 -
        independence = independence %|% 1,
137 -
        splicing_rate = splicing_rate %|% (log(2) / 2)
138 -
      )
139 -
    }
124 +
  lst(
125 +
    sampler_tfs = .kinetics_default_genes_sampler, 
126 +
    sampler_nontfs = .kinetics_default_genes_sampler, 
127 +
    sampler_interactions = .kinetics_default_interactions_sampler
128 +
  )
129 +
}
130 +
131 +
.kinetics_default_genes_sampler <- function(feature_info, feature_network, cache_dir = NULL, verbose = FALSE) {
132 +
  feature_info %>% mutate(
133 +
    transcription_rate = .data$transcription_rate %|% runif(n(), 10, 20),
134 +
    translation_rate = .data$translation_rate %|% runif(n(), 100, 150),
135 +
    mrna_halflife = .data$mrna_halflife %|% runif(n(), 2.5, 5),
136 +
    protein_halflife = .data$protein_halflife %|% runif(n(), 5, 10),
137 +
    independence = .data$independence %|% 1,
138 +
    splicing_rate = .data$splicing_rate %|% (log(2) / 2)
139 +
  )
140 +
}
141 +
142 +
.kinetics_default_interactions_sampler <- function(feature_info, feature_network, cache_dir = NULL, verbose = FALSE) {
143 +
  feature_network %>% mutate(
144 +
    effect = .data$effect %|% sample(c(-1L, 1L), n(), replace = TRUE, prob = c(.25, .75)),
145 +
    strength = .data$strength %|% 10 ^ runif(n(), log10(1), log10(100)),
146 +
    hill = .data$hill %|% rnorm_bounded(n(), 2, 2, min = 1, max = 10)
147 +
  )
148 +
}
149 +
150 +
#' @export
151 +
#' @rdname generate_kinetics
152 +
kinetics_random_distributions <- function() {
153 +
  .kinetics_default_genes_sampler <- function(feature_info, feature_network, cache_dir = NULL, verbose = FALSE) {
154 +
    feature_info %>% mutate(
155 +
      transcription_rate = .data$transcription_rate %|% runif_subrange(n(), 5, 30),
156 +
      translation_rate = .data$translation_rate %|% runif_subrange(n(), 50, 200),
157 +
      mrna_halflife = .data$mrna_halflife %|% runif_subrange(n(), 2, 6),
158 +
      protein_halflife = .data$protein_halflife %|% runif_subrange(n(), 4, 12),
159 +
      independence = .data$independence %|% 1,
160 +
      splicing_rate = .data$splicing_rate %|% (log(2) / 2)
161 +
    )
162 +
  }
140 163
  
141 -
  sampler_interactions <-
142 -
    function(feature_info, feature_network, cache_dir = NULL, verbose = FALSE) {
143 -
      feature_network %>% 
144 -
        mutate(
145 -
          effect = effect %|% sample(c(-1L, 1L), n(), replace = TRUE, prob = c(.25, .75)),
146 -
          strength = strength %|% 10 ^ runif(n(), log10(1), log10(100)),
147 -
          hill = hill %|% rnorm_bounded(n(), 2, 2, min = 1, max = 10)
148 -
        )
149 -
    }
164 +
  .kinetics_default_interactions_sampler <- function(feature_info, feature_network, cache_dir = NULL, verbose = FALSE) {
165 +
    effect_prob <- runif(1, .2, .8)
166 +
    
167 +
    feature_network %>% mutate(
168 +
      effect = .data$effect %|% sample(c(-1L, 1L), n(), replace = TRUE, prob = c(effect_prob, 1-effect_prob)),
169 +
      strength = .data$strength %|% 10 ^ runif_subrange(n(), log10(.5), log10(200)),
170 +
      hill = .data$hill %|% rnorm_bounded(n(), mean = 2, sd = 2, min = 1, max = 10)
171 +
    )
172 +
  }
150 173
  
151 -
  lst(sampler_tfs, sampler_nontfs = sampler_tfs, sampler_interactions)
174 +
  lst(
175 +
    sampler_tfs = .kinetics_default_genes_sampler,
176 +
    sampler_nontfs = .kinetics_default_genes_sampler,
177 +
    sampler_interactions = .kinetics_default_interactions_sampler
178 +
  )
152 179
}
153 180
154 181
155 182
.kinetics_generate_gene_kinetics <- function(model) {
156 -
  # satisfy r cmd check
157 -
  is_tf <- mrna_halflife <- protein_halflife <- to <- feature_id <- effect <- 
158 -
    basal <- basal_2 <- num_molecules <- mult <- id <- NULL
159 -
  
160 183
  if (model$verbose) cat("Generating kinetics for ", nrow(model$feature_info), " features\n", sep = "")
161 184
  params <- model$kinetics_params
162 185
  
@@ -168,13 +191,13 @@
Loading
168 191
  
169 192
  # generate relatively stable kinetics for TFs
170 193
  feature_info_tf <- params$sampler_tfs(
171 -
    feature_info %>% filter(is_tf),
194 +
    feature_info %>% filter(.data$is_tf),
172 195
    feature_network
173 196
  )
174 197
  
175 198
  # sample kinetics from dataset for non-tfs
176 199
  feature_info_nontf <- params$sampler_tfs(
177 -
    feature_info %>% filter(!is_tf),
200 +
    feature_info %>% filter(!.data$is_tf),
178 201
    feature_network, 
179 202
    cache_dir = model$download_cache_dir,
180 203
    verbose = model$verbose
@@ -184,8 +207,8 @@
Loading
184 207
  feature_info <- 
185 208
    bind_rows(feature_info_tf, feature_info_nontf) %>% 
186 209
    mutate(
187 -
      mrna_decay_rate = log(2) / mrna_halflife,
188 -
      protein_decay_rate = log(2) / protein_halflife
210 +
      mrna_decay_rate = log(2) / .data$mrna_halflife,
211 +
      protein_decay_rate = log(2) / .data$protein_halflife
189 212
    )
190 213
  
191 214
  # sample network kinetics from distributions
@@ -201,19 +224,19 @@
Loading
201 224
    left_join(
202 225
      feature_info,
203 226
      feature_network %>% 
204 -
        rename(feature_id = to) %>% 
205 -
        group_by(feature_id) %>% 
227 +
        rename(feature_id = .data$to) %>% 
228 +
        group_by(.data$feature_id) %>% 
206 229
        summarise(
207 -
          basal_2 = .kinetics_calculate_basal(effect)
230 +
          basal_2 = .kinetics_calculate_basal(.data$effect)
208 231
        ),
209 232
      by = "feature_id"
210 233
    ) %>% 
211 234
    mutate(
212 235
      # 1 for genes that are not being regulated by any other genes,
213 236
      # yet did not already have a value for 'basal' defined
214 -
      basal = basal %|% basal_2 %|% 1 
237 +
      basal = .data$basal %|% .data$basal_2 %|% 1 
215 238
    ) %>% 
216 -
    select(-basal_2)
239 +
    select(-.data$basal_2)
217 240
  
218 241
  model$feature_info <- feature_info
219 242
  model$feature_network <- feature_network
@@ -223,10 +246,6 @@
Loading
223 246
224 247
#' @importFrom GillespieSSA2 reaction
225 248
.kinetics_generate_formulae <- function(model) {
226 -
  # satisfy r cmd check
227 -
  from <- to <- `.` <- NULL
228 -
  
229 -
  
230 249
  if (model$verbose) cat("Generating formulae\n")
231 250
  
232 251
  # add helper information to feature info
@@ -234,14 +253,13 @@
Loading
234 253
    model$feature_info %>% 
235 254
    left_join(
236 255
      model$feature_network %>% 
237 -
        group_by(feature_id = to) %>% 
238 -
        do({tibble(regulators = list(.))}) %>% 
239 -
        ungroup(),
256 +
        group_by(feature_id = .data$to) %>% 
257 +
        summarise(regulators = list(data.frame(from = .data$from, effect = .data$effect, strength = .data$strength)), .groups = "drop"),
240 258
      by = "feature_id"
241 259
    ) %>% 
242 260
    left_join(
243 261
      model$feature_network %>% 
244 -
        group_by(feature_id = from) %>% 
262 +
        group_by(feature_id = .data$from) %>% 
245 263
        summarise(num_targets = n()),
246 264
      by = "feature_id"
247 265
    )
@@ -356,23 +374,22 @@
Loading
356 374
}
357 375
358 376
.kinetics_extract_parameters_as_df <- function(feature_info, feature_network) {
359 -
  # satisfy r cmd check
360 -
  feature_id <- transcription_rate <- splicing_rate <- translation_rate <- mrna_decay_rate <- protein_decay_rate <-
361 -
    basal <- independence <- param <- value <- id <- from <- to <- dissociation <- hill <- strength <- `.` <- from <- NULL
362 -
  
363 377
  # extract production / degradation rates, ind and bas
364 378
  feature_params <- 
365 379
    feature_info %>% 
366 -
    select(feature_id, transcription_rate, splicing_rate, translation_rate, mrna_decay_rate, protein_decay_rate, bas = basal, ind = independence) %>% 
367 -
    gather(param, value, -feature_id) %>% 
368 -
    mutate(id = paste0(param, "_", feature_id), type = "feature_info")
380 +
    select(
381 +
      .data$feature_id, .data$transcription_rate, .data$splicing_rate, .data$translation_rate,
382 +
      .data$mrna_decay_rate, .data$protein_decay_rate, bas = .data$basal, ind = .data$independence
383 +
    ) %>% 
384 +
    gather("param", "value", -.data$feature_id) %>% 
385 +
    mutate(id = paste0(.data$param, "_", .data$feature_id), type = "feature_info")
369 386
  
370 387
  # extract dis, hill, str
371 388
  edge_params <- 
372 389
    feature_network %>% 
373 -
    select(from, to, dis = dissociation, hill, str = strength) %>% 
374 -
    gather(param, value, -from, -to) %>% 
375 -
    mutate(id = paste0(param, "_", from, "_", to), type = "feature_network")
390 +
    select(.data$from, .data$to, dis = .data$dissociation, .data$hill, str = .data$strength) %>% 
391 +
    gather("param", "value", -.data$from, -.data$to) %>% 
392 +
    mutate(id = paste0(.data$param, "_", .data$from, "_", .data$to), type = "feature_network")
376 393
  
377 394
  bind_rows(feature_params, edge_params)
378 395
}
@@ -393,10 +410,6 @@
Loading
393 410
}
394 411
395 412
.kinetics_calculate_dissociation <- function(feature_info, feature_network) {
396 -
  # satisfy r cmd check
397 -
  transcription_rate <- mrna_decay_rate <- splicing_rate <- max_premrna <- translation_rate <- 
398 -
    protein_decay_rate <- max_mrna <- feature_id <- max_protein <- NULL
399 -
  
400 413
  remove <- c("max_premrna", "max_mrna", "max_protein", "dissociation", "k", "max_protein")
401 414
  
402 415
  feature_info <- feature_info[, !colnames(feature_info) %in% remove]
@@ -405,16 +418,16 @@
Loading
405 418
  feature_info <- 
406 419
    feature_info %>%
407 420
    mutate(
408 -
      max_premrna = transcription_rate / (mrna_decay_rate + splicing_rate),
409 -
      max_mrna = splicing_rate / mrna_decay_rate * max_premrna,
410 -
      max_protein = translation_rate / protein_decay_rate * max_mrna
421 +
      max_premrna = .data$transcription_rate / (.data$mrna_decay_rate + .data$splicing_rate),
422 +
      max_mrna = .data$splicing_rate / .data$mrna_decay_rate * .data$max_premrna,
423 +
      max_protein = .data$translation_rate / .data$protein_decay_rate * .data$max_mrna
411 424
    )
412 425
  
413 426
  feature_network <- 
414 427
    feature_network %>% 
415 -
    left_join(feature_info %>% select(from = feature_id, max_protein), by = "from") %>% 
428 +
    left_join(feature_info %>% select(from = .data$feature_id, .data$max_protein), by = "from") %>% 
416 429
    mutate(
417 -
      dissociation = max_protein / 2
430 +
      dissociation = .data$max_protein / 2
418 431
    )
419 432
  
420 433
  lst(feature_info, feature_network)

@@ -188,13 +188,7 @@
Loading
188 188
      lib_size <- lib_size_new[[cell_i]]
189 189
      
190 190
      # sample 'lib_size' molecules for each of the genes, weighted by 'gene_vals'
191 -
      tryCatch({
192 191
      rmultinom(1, lib_size, gene_vals)
193 -
      }, error = function(e) {
194 -
        cat("cell_i: ", cell_i, ", lib_size: ", lib_size, "\n", sep = "")
195 -
        print(gene_vals)
196 -
        print(e)
197 -
      })
198 192
    } else {
199 193
      integer(0)
200 194
    }

@@ -16,7 +16,13 @@
Loading
16 16
    if (!dir.exists(dirname(file))) {
17 17
      dir.create(dirname(file), recursive = TRUE)
18 18
    }
19 -
    download.file(url, destfile = file, quiet = !verbose)
19 +
    
20 +
    
21 +
    status <- suppressWarnings(
22 +
      utils::download.file(url, destfile = file, quiet = !verbose)
23 +
    )
24 +
    
25 +
    if (status != 0) stop("Cannot download file from ", url, call. = FALSE)
20 26
  }
21 27
  
22 28
  readRDS(file)

@@ -0,0 +1,43 @@
Loading
1 +
#' A bounded version of rnorm
2 +
#' 
3 +
#' @inheritParams stats::rnorm
4 +
#' @param min lower limits of the distribution.
5 +
#' @param max upper limits of the distribution.
6 +
#'
7 +
#' @importFrom stats pnorm qnorm
8 +
#' 
9 +
#' @return Generates values with rnorm, bounded by \[min, max\]
10 +
#' 
11 +
#' @export 
12 +
#' 
13 +
#' @examples
14 +
#' rnorm_bounded(10)
15 +
rnorm_bounded <- function(n, mean = 0, sd = 1, min = -Inf, max = Inf) {
16 +
  unif_min <- pnorm(min, mean = mean, sd = sd)
17 +
  unif_max <- pnorm(max, mean = mean, sd = sd)
18 +
  quan <- runif(n, unif_min, unif_max)
19 +
  qnorm(quan, mean = mean, sd = sd)
20 +
}
21 +
22 +
#' A subrange version of runif
23 +
#' 
24 +
#' Will generate numbers from a random subrange within the given range. 
25 +
#' For example, if \[min, max\]` is set to \[0, 10\], this function
26 +
#' could decide to generate `n` numbers between 2 and 6.
27 +
#' 
28 +
#' @param n Number of observations
29 +
#' @param min Lower limits of the distribution.
30 +
#' @param max Upper limits of the distribution.
31 +
#' 
32 +
#' @importFrom stats runif
33 +
#' 
34 +
#' @return Generates values with runif, bounded by a range drawn from `sort(runif(2, min, max))`.
35 +
#' 
36 +
#' @export
37 +
#' 
38 +
#' @examples 
39 +
#' runif_subrange(20, 0, 10)
40 +
runif_subrange <- function(n, min, max) {
41 +
  range <- sort(runif(2, min, max))
42 +
  runif(n, range[[1]], range[[2]])
43 +
}
Files Coverage
R 67.07%
Project Totals (22 files) 67.07%

No yaml found.

Create your codecov.yml to customize your Codecov experience

Sunburst
The inner-most circle is the entire project, moving away from the center are folders then, finally, a single file. The size and color of each slice is representing the number of statements and the coverage, respectively.
Icicle
The top section represents the entire project. Proceeding with folders and finally individual files. The size and color of each slice is representing the number of statements and the coverage, respectively.
Grid
Each block represents a single file in the project. The size and color of each block is represented by the number of statements and the coverage, respectively.
Loading