Simulate SNP pattern • tekkamaki

library(tekkamaki)
library(ggplot2)
grDevices::palette("Okabe-Ito")
options(
  ggplot2.continuous.colour = "viridis",
  ggplot2.continuous.fill = "viridis",
  ggplot2.discrete.colour = grDevices::palette()[-1],
  ggplot2.discrete.fill = grDevices::palette()[-1]
)

Generate a gene genealogy

result = tekka("--seed 42 -y40 -K100 -r2 -l2 --sa 2,2 --sj 2,2")
samples = result$sample_family[[1L]]
set.seed(666L)
genealogy = make_gene_genealogy(samples) |> print()

## $V  tibble [362 × 1] (S3: tbl_df/tbl/data.frame)
##  $ name: chr [1:362] "145-2" "7-1" "142-2" "7-2" ...
## # A tibble: 336 × 4
##    from  to    birth_year capture_year
##    <chr> <chr>      <int>        <int>
##  1 145-2 7-1           40           40
##  2 142-2 7-2           40           40
##  3 164-2 8-1           40           40
##  4 150-1 8-2           40           40
##  5 140-1 15-1          40           40
##  6 98-2  15-2          40           40
##  7 140-2 16-1          40           40
##  8 181-1 16-2          40           40
##  9 99-1  3-1           39           39
## 10 98-2  3-2           39           39
## # ℹ 326 more rows

Visualize the generated gene genealogy

plot(genealogy) +
  theme_void() +
  theme(legend.position = "top")

Generate SNPs

Randomly place a fixed number of mutations on a given genealogy without recombination.

snp = place_mutations(genealogy, 3L) |> print()

##      [,1] [,2] [,3]
## 7-1     0    0    0
## 7-2     0    0    0
## 8-1     0    0    0
## 8-2     0    0    0
## 15-1    1    0    0
## 15-2    0    0    0
## 16-1    0    0    0
## 16-2    0    0    0
## 3-1     0    0    0
## 3-2     0    0    0
## 4-1     0    0    1
## 4-2     0    0    0
## 11-1    0    0    0
## 11-2    0    1    1
## 12-1    0    0    1
## 12-2    0    0    0
## 5-1     0    0    0
## 5-2     0    0    0
## 6-1     0    0    0
## 6-2     0    0    0
## 1-1     0    0    0
## 1-2     0    0    1
## 2-1     0    0    1
## 2-2     0    0    0
## 13-1    0    0    0
## 13-2    0    0    0
## 14-1    0    0    0
## 14-2    0    0    0
## 9-1     1    0    0
## 9-2     0    0    1
## 10-1    1    0    0
## 10-2    1    0    0

Save SNP patterns as a VCF file:

write_vcf(snp, some_file)

## ##fileformat=VCFv4.5
## #CHROM   POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  7   8   15  16  3   4   11  12  5   6   1   2   13  14  9   10
## .    .   .   .   .   .   .   .   GT  0|0 0|0 1|0 0|0 0|0 0|0 0|0 0|0 0|0 0|0 0|0 0|0 0|0 0|0 1|0 1|1
## .    .   .   .   .   .   .   .   GT  0|0 0|0 0|0 0|0 0|0 0|0 0|1 0|0 0|0 0|0 0|0 0|0 0|0 0|0 0|0 0|0
## .    .   .   .   .   .   .   .   GT  0|0 0|0 0|0 0|0 0|0 1|0 0|1 1|0 0|0 0|0 0|1 1|0 0|0 0|0 0|1 0|0

You can also generate a list of SNP matrices for chromosomes with recombination in parallel. Use RNGkind("L'Ecuyer-CMRG") and set.seed() for reproducibility:

options(mc.cores = 2L) # parallel::detectCores(logical = FALSE)
RNGkind("L'Ecuyer-CMRG")
set.seed(666L)
snp = make_snp(samples, ss = c(chr7 = 4L, chr8 = 3L))
vcf = as_vcf(snp) |> add_pos_id()
write_vcf(vcf, some_file)

## ##fileformat=VCFv4.5
## #CHROM   POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  7   8   15  16  3   4   11  12  5   6   1   2   13  14  9   10
## chr7 1   chr7-1  .   .   .   .   .   GT  1|0 0|0 0|0 0|0 0|0 0|0 0|0 0|0 0|0 0|0 0|0 0|0 0|0 0|0 0|0 0|0
## chr7 2   chr7-2  .   .   .   .   .   GT  0|0 0|0 0|1 0|0 0|0 0|0 0|0 0|0 0|0 0|0 0|0 1|0 0|0 0|0 0|0 0|0
## chr7 3   chr7-3  .   .   .   .   .   GT  0|0 0|0 0|0 0|0 0|0 1|0 0|0 0|0 0|0 0|0 0|0 0|0 0|0 0|0 0|0 0|0
## chr7 4   chr7-4  .   .   .   .   .   GT  1|1 0|1 1|1 1|0 0|1 1|1 1|0 1|0 1|0 1|0 0|1 1|1 1|1 1|0 1|1 1|1
## chr8 1   chr8-1  .   .   .   .   .   GT  0|0 0|0 0|0 1|0 0|0 0|0 0|0 0|0 0|0 1|1 0|0 0|0 0|0 0|0 1|0 0|0
## chr8 2   chr8-2  .   .   .   .   .   GT  0|0 0|0 0|0 0|0 0|1 0|0 0|0 1|0 0|0 0|0 0|0 0|0 0|0 0|0 0|0 0|0
## chr8 3   chr8-3  .   .   .   .   .   GT  0|0 0|1 0|0 0|0 0|0 0|0 0|0 0|0 0|0 0|0 0|0 0|0 0|0 0|0 0|0 0|0