Workflow
Perform simulation and extract population data:
library(ggplot2)
library(tumopp)
set.seed(24601L)
result = tumopp::tumopp("-N20000 -D3 -Chex -Lconst -k10")
population = result$population[[1L]]
graph = result$graph[[1L]]
Sample cells and put neutral mutations on the lineages:
extant = population |> tumopp::filter_extant()
ncell = 200L
regions = tumopp::sample_uniform_regions(extant, nsam = 4L, ncell = ncell)
subgraph = tumopp::subtree(graph, unlist(regions$id))
vaf = tumopp::make_vaf(subgraph, regions$id, mu = 8.0) |> print()
## # A tibble: 29,094 × 4
## `1` `2` `3` `4`
## <dbl> <dbl> <dbl> <dbl>
## 1 0 0 0.005 0
## 2 0 0 0.005 0
## 3 0 0 0.005 0
## 4 0 0 0.005 0
## 5 0 0 0.005 0
## 6 0 0 0.005 0
## 7 0 0 0.005 0
## 8 0 0 0.005 0
## 9 0 0 0.005 0
## 10 0 0 0.005 0
## # ℹ 29,084 more rows
Estimate from VAF.
tumopp::dist_vaf(vaf, ncell) |> tumopp::fst()
## [1] 0.0440782
# True FST from cell genealogy
tumopp::dist_genealogy(subgraph, regions$id) |> tumopp::fst()
## [1] 0.04560559
Summarize and visualize VAF:
vaf_tidy = vaf |>
tumopp::filter_detectable(0.01) |>
tumopp::sort_vaf() |>
tumopp::longer_vaf() |>
print()
ggplot(vaf_tidy) +
aes(sample, site) +
geom_tile(aes(fill = frequency)) +
scale_fill_distiller(palette = "Spectral", limit = c(0, 1), guide = FALSE) +
coord_cartesian(expand = FALSE)