Simulate SNP pattern • tekkamaki

library(conflicted)
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

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 [590 × 1] (S3: tbl_df/tbl/data.frame)
##  $ name: chr [1:590] "196-1" "7-1" "186-1" "7-2" ...
## # A tibble: 588 × 4
##    from  to    birth_year capture_year
##    <chr> <chr>      <int>        <int>
##  1 196-1 7-1           40           40
##  2 186-1 7-2           40           40
##  3 220-2 8-1           40           40
##  4 208-2 8-2           40           40
##  5 277-1 15-1          40           40
##  6 271-1 15-2          40           40
##  7 288-2 16-1          40           40
##  8 279-1 16-2          40           40
##  9 131-2 3-1           39           39
## 10 115-1 3-2           39           39
## # ℹ 578 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()

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