Restrict features to a subgroup of loci (#181)

* dump version

* add infrastructure for storing locus groups

* add locus_group argument to features

* restrict group models to their features

* add example

* fix lints

* update NEWS

DESCRIPTION

@@ -1,5 +1,5 @@
 Package: coala
-Version: 0.4.1.9005
+Version: 0.4.1.9006
 License: MIT + file LICENSE
 Title: A Framework for Coalescent Simulation
 Authors@R: c(

NEWS.md

@@ -4,6 +4,9 @@ coala 0.5.0
 * Major interal refactoring on how simulators interface with coala (#174).
 * Support for calculating an expanded version of MCMF (#173, #179). This
   feature was contributed by Jorge E. Amaya Romero (@jorgeamaya).
+* Introduces the optional `locus_group` argument for features. Using it, 
+  features can be defined only for a subset of the loci in the model
+  (#161, #181). Thanks to @andrewparkermorgan for suggesting this feature.
 
 
 

R/feature.R

@@ -6,12 +6,17 @@ feature_class <- R6Class("feature", inherit = model_part,
     population = NULL,
     time = NULL,
     rate = NA,
+    locus_group = NULL,
     set_population = function(population, expected_length = 1) {
       assert_that(identical(population, "all") || is.numeric(population))
       assert_that(identical(population, "all") ||
                     length(population) == expected_length)
       private$population <- population
     },
+    set_locus_group = function(locus_group) {
+      assert_that(identical(locus_group, "all") || is.numeric(locus_group))
+      private$locus_group <- locus_group
+    },
     add_parameter = function(parameter, required = TRUE, add_par = TRUE) {
       expr <- NA
       if (is.numeric(parameter) && length(parameter) == 1) {
@@ -31,14 +36,16 @@ feature_class <- R6Class("feature", inherit = model_part,
     }
   ),
   public = list(
-    initialize = function(rate, population, time) {
+    initialize = function(rate, population, time, locus_group) {
       if (!missing(rate)) private$rate <- private$add_parameter(rate)
       if (!missing(time)) private$time <- private$add_parameter(time)
       if (!missing(population)) private$set_population(population)
+      if (!missing(locus_group)) private$set_locus_group(locus_group)
     },
     get_parameters = function() private$parameter,
     get_population = function() private$population,
     get_time = function() private$time,
+    get_locus_group = function() private$locus_group,
     reset_parameters = function() private$parameter <- list(),
     print = function() {
       cat("Feature of type", private$feature_table$type[1], "\n")

R/feature_growth.R

@@ -30,7 +30,7 @@ growth_class <- R6Class("growth", inherit = feature_class,
 #' @export
 #' @seealso For instantaneous population size
 #'          changes: \code{\link{feat_size_change}}
-#' @family features
+#' @template feature
 #' @examples
 #' # Simulate a haploid population that has been expanding for
 #' # the last 2*Ne generations
@@ -40,8 +40,14 @@ growth_class <- R6Class("growth", inherit = feature_class,
 #'   feat_mutation(10) +
 #'   sumstat_sfs()
 #' simulate(model)
-feat_growth <- function(rate, population = "all", time="0") {
-  growth_class$new(rate, population, time)
+feat_growth <- function(rate,
+                        population = "all",
+                        time = "0",
+                        locus_group = "all") {
+
+  growth_class$new(rate, population, time,
+                   locus_group = locus_group)
+
 }
 
 

R/feature_ignore_singletons.R

@@ -13,9 +13,8 @@ ign_singletons_class <- R6Class("ign_singletons", inherit = feature_class)
 #' allele is observed exactly once in all sequences, regardless of the
 #' population structure.
 #'
-#' @return The feature, which can be added to a model using `+`.
 #' @export
-#' @family features
+#' @template feature
 #' @examples
 #' model <- coal_model(2, 1) +
 #'   feat_mutation(10) +
@@ -24,8 +23,8 @@ ign_singletons_class <- R6Class("ign_singletons", inherit = feature_class)
 #' # In this model, all mutations are singletons. Therefore,
 #' # the number of mutations is 0 when removing singletons:
 #' simulate(model)$n_mut
-feat_ignore_singletons <- function() {
-  ign_singletons_class$new()
+feat_ignore_singletons <- function(locus_group = "all") {
+  ign_singletons_class$new(locus_group = locus_group)
 }
 
 

R/feature_migration.R

@@ -1,8 +1,9 @@
 migration_class <- R6Class("migration", inherit = feature_class,
   private = list(rate = NA),
   public = list(
-    initialize = function(rate, pop_from, pop_to, time, symmetric = FALSE) {
-      super$initialize(rate = rate, time = time)
+    initialize = function(rate, pop_from, pop_to, time,
+                          symmetric = FALSE, locus_group) {
+      super$initialize(rate = rate, time = time, locus_group = locus_group)
 
       if (symmetric) {
         private$population <- "all"
@@ -48,7 +49,7 @@ migration_class <- R6Class("migration", inherit = feature_class,
 #'        rate is set. The rate applies to the time past warts
 #'        of the time point, until it is changed again.
 #' @export
-#' @family features
+#' @template feature
 #'
 #' @examples
 #' # Asymmetric migration between two populations:
@@ -67,15 +68,18 @@ migration_class <- R6Class("migration", inherit = feature_class,
 #'   sumstat_sfs()
 #' simulate(model)
 feat_migration <- function(rate, pop_from = NULL, pop_to = NULL,
-                           symmetric = FALSE, time = "0") {
+                           symmetric = FALSE, time = "0",
+                           locus_group = "all") {
   if (symmetric) {
     if (!(is.null(pop_from) && is.null(pop_to))) {
       warning("Ignoring 'pop_form' and 'pop_to' because 'symmetric' is TRUE")
     }
-    return(migration_class$new(rate, time = time, symmetric = TRUE))
+    return(migration_class$new(rate, time = time, symmetric = TRUE,
+                               locus_group = locus_group))
   }
 
-  migration_class$new(rate, pop_from, pop_to, time)
+  migration_class$new(rate, pop_from, pop_to, time,
+                      locus_group = locus_group)
 }
 
 

R/feature_mutation.R

@@ -9,8 +9,10 @@ mutation_class <- R6Class("mutation", inherit = feature_class,
   ),
   public = list(
     initialize = function(rate, model, base_frequencies,
-                          tstv_ratio, gtr_rates, fixed) {
+                          tstv_ratio, gtr_rates, fixed,
+                          locus_group) {
       private$rate <- private$add_parameter(rate, add_par = FALSE)
+      private$set_locus_group(locus_group)
 
       assert_that(length(model) == 1)
       assert_that(any(model == c("IFS", "HKY", "GTR")))
@@ -91,7 +93,7 @@ mutation_class <- R6Class("mutation", inherit = feature_class,
 #' @seealso For using rates that variate between the loci in a model:
 #'   \code{\link{par_variation}}, \code{\link{par_zero_inflation}}
 #' @seealso For adding recombination: \code{\link{feat_recombination}}.
-#' @family features
+#' @template feature
 #'
 #' @section Mutation Models:
 #' The infinite sites mutation (\strong{IFS}) model is a frequently used simplification
@@ -119,6 +121,16 @@ mutation_class <- R6Class("mutation", inherit = feature_class,
 #' model <- coal_model(5, 1) + feat_mutation(5) + sumstat_seg_sites()
 #' simulate(model)
 #'
+#' # A model with a mutation of 5.0 for the first 10 loci, and 7.5 for the
+#' # second 10 loci
+#' model <- coal_model(4) +
+#'   locus_averaged(10, 100) +
+#'   locus_averaged(10, 100) +
+#'   feat_mutation(5.0, locus_group = 1) +
+#'   feat_mutation(7.5, locus_group = 2) +
+#'   sumstat_seg_sites()
+#' simulate(model)
+#'
 #' # A model with 7 mutations per locus:
 #' model <- coal_model(5, 1) +
 #'   feat_mutation(7, fixed = TRUE) +
@@ -147,10 +159,12 @@ feat_mutation <- function(rate,
                           base_frequencies = NA,
                           tstv_ratio = NA,
                           gtr_rates = NA,
-                          fixed_number = FALSE) {
+                          fixed_number = FALSE,
+                          locus_group = "all") {
 
   mutation_class$new(rate, model, base_frequencies,
-                     tstv_ratio, gtr_rates, fixed_number)
+                     tstv_ratio, gtr_rates, fixed_number,
+                     locus_group = locus_group)
 }
 
 is_feat_mutation <- function(feat) any("mutation" == class(feat))

R/feature_outgroup.R

@@ -1,7 +1,8 @@
 outgroup_class <- R6Class("outgroup", inherit = feature_class,
   public = list(
-    initialize = function(population) {
+    initialize = function(population, locus_group) {
       private$set_population(population)
+      private$set_locus_group(locus_group)
     },
     print = function() {
       cat("Outgroup: Population", private$population, "\n")
@@ -19,16 +20,16 @@ outgroup_class <- R6Class("outgroup", inherit = feature_class,
 #'
 #' @param population The population that is marked as outgroup.
 #' @export
-#' @family features
+#' @template feature
 #' @examples
 #' # A simple finite sites model
 #' model <- coal_model(c(4, 6, 1), 2, 10) +
 #'    feat_outgroup(3) +
 #'    feat_pop_merge(0.5, 2, 1) +
 #'    feat_pop_merge(2, 3, 1) +
 #'    feat_mutation(5, model = "GTR", gtr_rates = 1:6)
-feat_outgroup <- function(population) {
-  outgroup_class$new(population)
+feat_outgroup <- function(population, locus_group = "all") {
+  outgroup_class$new(population, locus_group = locus_group)
 }
 
 is_feat_outgroup <- function(feat) inherits(feat, "outgroup")

R/feature_pop_merge.R

@@ -1,8 +1,9 @@
 pop_merge_class <- R6Class("pop_merge", inherit = feature_class,
   public = list(
-    initialize = function(time, pop_from, pop_to) {
+    initialize = function(time, pop_from, pop_to, locus_group) {
       private$time <- private$add_parameter(time)
       private$set_population(c(from = pop_from, to = pop_to), 2)
+      private$set_locus_group(locus_group)
     },
     print = function() {
       cat("Merge of pop", private$population[1],
@@ -30,7 +31,7 @@ pop_merge_class <- R6Class("pop_merge", inherit = feature_class,
 #'        the population in which the speciation event occurs.
 #' @param time The time at which the merge occurs.
 #' @export
-#' @family features
+#' @template feature
 #' @examples
 #' # Two population which merge after 0.5 time units:
 #' model <- coal_model(c(25,25), 1) +
@@ -43,8 +44,9 @@ pop_merge_class <- R6Class("pop_merge", inherit = feature_class,
 #' model_iwm <- model +
 #'   feat_migration(.75, symmetric = TRUE)
 #' simulate(model_iwm)
-feat_pop_merge <- function(time, pop_source, pop_target) {
-  pop_merge_class$new(time, pop_source, pop_target)
+feat_pop_merge <- function(time, pop_source, pop_target, locus_group = "all") {
+  pop_merge_class$new(time, pop_source, pop_target,
+                      locus_group = locus_group)
 }
 
 #' @describeIn conv_to_ms_arg Feature conversion

R/feature_recombination.R

@@ -21,7 +21,7 @@ recombination_class <- R6Class("recombination", inherit = feature_class,
 #'        recombination event within the locus occurs in one generation.
 #' @return The feature, which can be added to a model using `+`.
 #' @export
-#' @family features
+#' @template feature
 #' @seealso For adding recombination: \code{\link{feat_mutation}}.
 #'
 #' @examples
@@ -31,8 +31,8 @@ recombination_class <- R6Class("recombination", inherit = feature_class,
 #'   feat_mutation(5) +
 #'   sumstat_sfs()
 #' simulate(model)
-feat_recombination <- function(rate) {
-  recombination_class$new(rate)
+feat_recombination <- function(rate, locus_group = "all") {
+  recombination_class$new(rate, locus_group = locus_group)
 }
 
 #' @describeIn conv_to_ms_arg Feature conversion

R/feature_sample.R

@@ -2,6 +2,8 @@ sample_class <- R6Class("sample", inherit = feature_class,
   private = list(size = NA, ploidy = NA),
   public = list(
     initialize = function(sizes, ploidy, time) {
+      private$set_locus_group("all")
+
       assert_that(is.numeric(sizes))
       assert_that(length(sizes) >= 1)
       private$size <- sizes
@@ -36,8 +38,8 @@ sample_class <- R6Class("sample", inherit = feature_class,
 #' @param ploidy The number of chromosomes that will be simulated per
 #'   individual.
 #' @param time The time at which the sample is taken.
-#' @return The feature, which can be added to a model using `+`.
-#' @keywords internal
+#' @keywords interna
+# @template feature
 feat_sample <- function(individuals, ploidy = 1, time = "0") {
   if (time != "0")
     stop("Samples at time different from 0 at not supported at the moment")

R/feature_selection.R

@@ -13,7 +13,8 @@ selection_class <- R6Class("selection", inherit = feature_class,
   public = list(
     initialize = function(strength_AA, strength_Aa, strength_aa,
                           population, time, additive = FALSE,
-                          start, start_frequency, Ne, position, force_keep) {
+                          start, start_frequency, Ne, position, force_keep,
+                          locus_group) {
 
       assert_that(is.logical(additive) && length(additive) == 1)
       private$additive <- additive
@@ -41,6 +42,7 @@ selection_class <- R6Class("selection", inherit = feature_class,
 
       private$time <- private$add_parameter(time)
       private$set_population(population)
+      private$set_locus_group(locus_group)
     },
     get_strength_AA = function() private$strength_AA,
     get_strength_Aa = function() private$strength_Aa,
@@ -116,7 +118,7 @@ selection_class <- R6Class("selection", inherit = feature_class,
 #' @seealso For summary statistics that are sensitive for selection:
 #'   \code{\link{sumstat_tajimas_d}}, \code{\link{sumstat_ihh}},
 #'   \code{\link{sumstat_omega}}, \code{\link{sumstat_mcmf}}
-#' @family features
+#' @template  feature
 #' @examples
 #' # Positive additive selection in population 2:
 #' model <- coal_model(c(10, 13), 1, 10000) +
@@ -138,7 +140,8 @@ feat_selection <- function(strength_AA = 0,
                            start_frequency = 0.0005,
                            Ne = 10000,
                            position = 0.5,
-                           force_keep = TRUE) {
+                           force_keep = TRUE,
+                           locus_group = "all") {
 
   if (!is.null(strength_A)) {
     return(selection_class$new(strength_Aa = strength_A,
@@ -149,11 +152,13 @@ feat_selection <- function(strength_AA = 0,
                                start_frequency = start_frequency,
                                Ne = Ne,
                                position = position,
-                               force_keep = force_keep))
+                               force_keep = force_keep,
+                               locus_group = locus_group))
   }
 
   selection_class$new(strength_AA, strength_Aa, strength_aa, population, time,
-                      FALSE, start, start_frequency, Ne, position, force_keep)
+                      FALSE, start, start_frequency, Ne, position, force_keep,
+                      locus_group = locus_group)
 }
 
 

R/feature_size_change.R

@@ -26,7 +26,7 @@ size_change_class <- R6Class("size_change", inherit = feature_class,
 #' @return The feature, which can be added to a model using `+`.
 #' @export
 #' @seealso For continuous size changes over time: \code{\link{feat_growth}}.
-#' @family features
+#' @template feature
 #' @examples
 #' # A model with one smaller population:
 #' model <- coal_model(c(20, 5), 3) +
@@ -43,8 +43,10 @@ size_change_class <- R6Class("size_change", inherit = feature_class,
 #'   feat_mutation(20) +
 #'   sumstat_sfs()
 #' simulate(model)
-feat_size_change <- function(new_size, population = 1, time = "0") {
-  size_change_class$new(new_size, population, time)
+feat_size_change <- function(new_size, population = 1, time = "0",
+                             locus_group = "all") {
+  size_change_class$new(new_size, population, time,
+                        locus_group = locus_group)
 }
 
 #' @describeIn conv_to_ms_arg Feature conversion

R/feature_sumstats.R

@@ -1,5 +1,6 @@
 segsites_feat_class <- R6Class("seg_sites_feat", inherit = feature_class,
   public = list(
+    initialize = function() super$initialize(locus_group = "all"),
     print = function() cat("Generating Seg. Sites\n")
   )
 )
@@ -29,6 +30,7 @@ conv_to_seqgen_arg.seg_sites_feat <- conv_to_ms_arg.seg_sites_feat #nolint
 
 trees_feat_class <- R6Class("trees_feat", inherit = feature_class,
   public = list(
+    initialize = function() super$initialize(locus_group = "all"),
     print = function() cat("Generating Trees\n")
   )
 )
@@ -55,6 +57,7 @@ conv_to_seqgen_arg.trees_feat <- function(feature, model) {
 
 files_feat_class <- R6Class("files_feat", inherit = feature_class,
   public = list(
+    initialize = function() super$initialize(locus_group = "all"),
     print = function() cat("Generating Files\n")
   )
 )