Simulate SNP pattern • tekkamaki

library(conflicted)
library(tekkamaki)
library(ggplot2)
palette_discrete = grDevices::palette.colors(palette = "Okabe-Ito")[-1]
theme_palette = ggplot2::theme(
  palette.colour.continuous = "viridis",
  palette.fill.continuous = "viridis",
  palette.colour.discrete = palette_discrete,
  palette.fill.discrete = palette_discrete
)

Generate a gene genealogy

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

## $V  tibble [722 × 1] (S3: tbl_df/tbl/data.frame)
##  $ name: chr [1:722] "211-1" "5-1" "192-1" "5-2" ...
## # A tibble: 720 × 4
##    from  to    birth_year capture_year
##    <chr> <chr>      <int>        <int>
##  1 211-1 5-1           40           40
##  2 192-1 5-2           40           40
##  3 242-2 6-1           40           40
##  4 219-2 6-2           40           40
##  5 322-1 13-1          40           40
##  6 314-1 13-2          40           40
##  7 342-2 14-1          40           40
##  8 337-1 14-2          40           40
##  9 87-2  1-1           39           39
## 10 17-1  1-2           39           39
## # ℹ 710 more rows

Visualize the generated gene genealogy

plot(genealogy) +
  theme_void() +
  theme_palette +
  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()

##           
## 5-1  0 0 0
## 5-2  0 0 0
## 6-1  0 0 0
## 6-2  0 0 0
## 13-1 1 0 0
## 13-2 0 1 0
## 14-1 0 0 0
## 14-2 0 0 0
## 1-1  0 0 0
## 1-2  0 0 0
## 2-1  0 0 0
## 2-2  0 0 0
## 9-1  0 0 0
## 9-2  0 0 0
## 10-1 0 0 0
## 10-2 0 0 0
## 7-1  0 0 0
## 7-2  0 0 0
## 8-1  0 0 0
## 8-2  0 0 0
## 3-1  0 0 1
## 3-2  0 0 0
## 4-1  0 0 0
## 4-2  0 0 0
## 12-1 0 0 0
## 12-2 0 0 0
## 15-1 0 0 0
## 15-2 0 0 0
## 16-1 0 0 0
## 16-2 0 0 0
## 11-1 0 0 0
## 11-2 0 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  5   6   13  14  1   2   9   10  7   8   3   4   12  15  16  11
## .    .   .   .   .   .   .   .   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 0|0 0|0
## .    .   .   .   .   .   .   .   GT  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|0 0|0 0|0
## .    .   .   .   .   .   .   .   GT  0|0 0|0 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 0|0

You can also generate a list of SNP data.frame 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  5   6   13  14  1   2   9   10  7   8   3   4   12  15  16  11
## chr7 1   chr7-1  .   .   .   .   .   GT  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 0|0
## chr7 2   chr7-2  .   .   .   .   .   GT  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 0|0
## chr7 3   chr7-3  .   .   .   .   .   GT  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|1
## chr7 4   chr7-4  .   .   .   .   .   GT  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|0
## chr8 1   chr8-1  .   .   .   .   .   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
## chr8 2   chr8-2  .   .   .   .   .   GT  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 0|0 0|0 0|0
## chr8 3   chr8-3  .   .   .   .   .   GT  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 0|0 0|0 0|0