knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

library(tumopp)
set.seed(24601L)
result = tumopp::tumopp("-N20000 -D3 -Chex -Lconst -k10 --seed=42")
population = result$population[[1L]]
graph = result$graph[[1L]]

Workflow

Perform simulation and extract population data:

library(tumopp)
library(dplyr)
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()

Estimate $F_{ST}$ from VAF.

tumopp::dist_vaf(vaf, ncell) |> tumopp::fst()

# True FST from cell genealogy
tumopp::dist_genealogy(subgraph, regions$id) |> tumopp::fst()

Summarize and visualize VAF:

vaf_tidy = vaf |>
  tumopp::filter_detectable(0.01) |>
  tumopp::sort_vaf() |>
  tumopp::longer_vaf() |>
  print()

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


heavywatal/rtumopp documentation built on March 30, 2024, 11:08 a.m.