R-package coala: how to create a stepping stone migration model

Question

[dsw.dejonge]

I am using the coala R-package to do coalescent simulations, and I was wondering if someone knows how to easily implement a stepping-stone migration model?

A reproducible example: 4 linearly distributed populations exchange migrants according to stepping-stone pattern (only the adjacent populations).

model <- coal_model(sample_size = c(5, 5, 5, 5), # sample 5 ind from every pop
                                   loci_number = 1,
                                   loci_length = 10,
                                    ploidy = 1) +
feat_mutation(rate = mut_rate, # e.g. 0.1
                       model = "HKY",
                       base_frequencies = c(0.25,0.25,0.25,0.25),
                       tstv_ratio = 4) +
feat_migration(mig_rate, 1, 2) + # mig_rate can be e.g. 0.5
feat_migration(mig_rate, 2, 1) +
feat_migration(mig_rate, 2, 3) +
feat_migration(mig_rate, 3, 2) +
feat_migration(mig_rate, 3, 4) +
feat_migration(mig_rate, 4, 3) +
sumstat_dna(name = "dna", transformation = identity)

This example works, but the downside is that I have to write two new lines for every pair of populations that exchange individuals. It is fine for a small number of populations, but I want to do a large simulation with about 70 populations. Does someone has a good idea how to automate this? The documentation has not helped me so far.

I tried two things that didn't work:

feat_migration(mig_rate, c(1,2,2,3,3,4), c(2,1,3,2,4,3))

and something like this:

migration_model <- function(){
  for(i in 1:n_pops){
    feat_migration(mig_rate, i, i+1) +
    feat_migration(mig_rate, i+1, i))
}

In the latter case, I don't really know how I can correctly create and parse all functions correctly into my model.

Good ideas are very welcome! :)

[Ram]

You can look into do.call as a means to build a function call. Lexical scoping works very well in R so you shouldn't have a problem generating a function call as plain text and then executing it as an actual function call.

[dsw.dejonge]

Thank you for your suggestion! I've never used this function before, and I tried to implement it as follows:

model <- coal_model(<params>) +
               feat_mutation(<params>) +
               do.call(feat_migration, list(mig_rate, 
                           pop_from = c(1,2,2,3,3,4), 
                           pop_to = c(2,1,3,2,4,3))) +
               sumstat_dna(<params>)

However, I get an error from the package "Error: population not identical to "all" or length(population) not equal to expected_length" (no such error with the original code). Am I implementing the do.call function right? Should I do something about quoting the arguments or the environment maybe?

[Ram]

Play around with it. You might need to use apply to ensure only one of each of the pop_from and pop_to lists are passed per feat_mutation call.

[Ram]

Hello dsw.dejonge!

It appears that your post has been cross-posted to another site: https://stackoverflow.com/questions/54042332/how-to-automatically-set-up-and-add-functions-to-a-model-in-r

This is typically not recommended as it runs the risk of annoying people in both communities.

Answer 1

[dsw.dejonge]

Instead of using do.call() I used Map() and Reduce() as proposed by user Parfait on Stackoverflow.

n_pops <- 4    
start_pts <- as.vector(sapply(seq(n_pops-1), function(x) c(x, x+1)))  
end_pts <- as.vector(sapply(seq(n_pops-1), function(x) c(x+1, x)))

# LIST OF feat_migration()
feats <- Map(function(x, y) feat_migration(mig_rate, x, y), start_pts, end_pts)

# LIST OF FUNCTIONS
funcs <- c(list(coal_model(sample_size = c(5, 5, 5, 5),
                           loci_number = 1,
                           loci_length = 10,
                           ploidy = 1),
                feat_mutation(rate = mut_rate, # e.g. 0.1
                              model = "HKY",
                              base_frequencies = c(0.25,0.25,0.25,0.25),
                              tstv_ratio = 4),
                sumstat_dna(name = "dna", transformation = identity)),

            feats)
           )

# MODEL CALL     
model <- Reduce(`+`, funcs)