## Error in get(paste0(generic, ".", class), envir = get_method_env()) :
## object 'type_sum.accel' not found
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 [556 × 1] (S3: tbl_df/tbl/data.frame)
## $ name: chr [1:556] "183-1" "7-1" "164-1" "7-2" ...
## # A tibble: 554 × 4
## from to birth_year capture_year
## <chr> <chr> <int> <int>
## 1 183-1 7-1 40 40
## 2 164-1 7-2 40 40
## 3 131-2 8-1 40 40
## 4 205-2 8-2 40 40
## 5 266-1 15-1 40 40
## 6 259-1 15-2 40 40
## 7 274-2 16-1 40 40
## 8 269-1 16-2 40 40
## 9 131-2 3-1 39 39
## 10 115-1 3-2 39 39
## # ℹ 544 more rowsVisualize 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 1 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 0
## 11-1 0 0 0
## 11-2 0 0 0
## 12-1 0 0 0
## 12-2 0 0 0
## 1-1 0 0 0
## 1-2 0 0 0
## 2-1 0 0 0
## 2-2 0 0 0
## 5-1 0 0 0
## 5-2 0 0 0
## 6-1 0 0 0
## 6-2 1 0 0
## 9-1 0 0 0
## 9-2 0 0 0
## 13-1 0 0 1
## 13-2 0 0 0
## 14-1 0 0 0
## 14-2 0 0 0
## 10-1 0 0 0
## 10-2 0 0 0Save 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 1 2 5 6 9 13 14 10
## . . . . . . . . GT 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 0|0
## . . . . . . . . GT 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 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|0 1|0 0|0 0|0You 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 1 2 5 6 9 13 14 10
## chr7 1 chr7-1 . . . . . 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
## chr7 2 chr7-2 . . . . . GT 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 0|0
## chr7 3 chr7-3 . . . . . GT 0|1 0|0 0|0 0|0 0|0 1|0 0|0 0|0 0|0 0|0 0|0 1|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 0|0 0|0 0|0 0|0 0|1 0|0 0|0 0|0 0|0
## chr8 1 chr8-1 . . . . . 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
## chr8 2 chr8-2 . . . . . GT 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|0
## chr8 3 chr8-3 . . . . . 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|1 0|0 0|0 0|0