---
title: "Simple Example of a Simulation using Nucleotides"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Simple Example of a Simulation using Nucleotides}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, eval = slimr::slim_is_avail())
library(slimr)
```

```{r nuc_sim}

## set some parameters
seed <- 1205
split_prob <- 0.001
max_subpops <- 10

## specify simulation
split_isolate_sim <- slim_script(
  
  slim_block(initialize(), {
  
    setSeed(!!seed);
    
    ## tell SLiM to simulate nucleotides
    initializeSLiMOptions(nucleotideBased=T);
    initializeAncestralNucleotides(randomNucleotides(1000));
    initializeMutationTypeNuc("m1", 0.5, "f", 0.0);
    
    initializeGenomicElementType("g1", m1, 1.0, mmJukesCantor(1e-5));
    initializeGenomicElement(g1, 0, 1000 - 1);
    initializeRecombinationRate(1e-8);
    
  }),
  
  slim_block(1, early(), {
    
    defineGlobal("curr_subpop", 1);
    sim.addSubpop(curr_subpop, 100)
    
  }),
  
  slim_block(2, 10000, late(), {
    
    if(rbinom(1, 1, !!split_prob) == 1) {
      ## split a subpop
      subpop_choose = sample(sim.subpopulations, 1)
      defineGlobal("curr_subpop", curr_subpop + 1)
      
      sim.addSubpopSplit(subpopID = curr_subpop, 
                         size = 100, 
                         sourceSubpop = subpop_choose)
      ## if too many subpops, remove one randomly
      if(size(sim.subpopulations) > !!max_subpops) {
        subpop_del = sample(sim.subpopulations, 1)
        subpop_del.setSubpopulationSize(0)
      }
    }
    
    slimr_output_nucleotides(subpops = TRUE, do_every = 100)
      
  }),
  
  slim_block(10000, late(), {
    sim.simulationFinished()
  })
  
)

results <- slim_run(split_isolate_sim, throw_error = TRUE)

res_data <- slim_results_to_data(results)

res_data

## sequences at generation 10000
image(ape::as.DNAbin(res_data$data[[100]]))

```

And then we can use some other R packages to quickly build a tree based on the simulated nucleotides, to see if it looks like what we would expect from a sequentially splitting population.

```{r nj_tree2, echo=TRUE}
## convert to ape::DNAbin
al <- ape::as.DNAbin(res_data$data[[100]])
dists <- ape::dist.dna(al)
upgma_tree <- ape::as.phylo(hclust(dists, method = "average"))
pal <- paletteer::paletteer_d("RColorBrewer::Paired", 10)
plot(upgma_tree, show.tip.label = FALSE)
ape::tiplabels(pch = 19, col = pal[as.numeric(as.factor(res_data$subpops[[100]]))])
```