knitr::opts_chunk$set(echo = TRUE) library(qfm)
In this example, we demonstrate quantitative fate mapping on lineage barcode data measured with cell types. First, we apply Phylotime to reconstruct time-scaled phylogeny from lineage barcodes. Second, ICE-FASE is applied to extract parameters of progenitor population dynamics, including commitment times, progenitor population sizes and commitment biases.
We start with a character matrix, the example data is in the format below. An cell id is required, and is provided in the first column.
# read mutant allele matrix cell_allele_table = readr::read_csv("../data/example/mutant_allele_matrix_with_missing.csv") cell_allele_table[1:5]
First, we convert the table to a character matrix
chr_mat = as.matrix(cell_allele_table[-1]) rownames(chr_mat) = cell_allele_table$cell
Generally, mutation matrix contains missing data. We will use impute_characeters
to impute the missing data.
chr_mat_imputed = impute_characters(chr_mat, nrounds = 100)
The input to Phylotime is a character matrix, with rows corresponding to cells and columns corresponding to barcoding sites, tibble / data.frame are convert to a character matrix. The matrix needs to have rownames, which will be used as the names of the cells.
Within each barcoding site (column), the same character denotes a unique mutant allele. The unmutated allele needs to be coded as the character "0". (zero) The user also needs to input the total amount of active barcoding time, which is the duration for which mutations are being actively generated.
Phylotime estimates mutation rate from data provided and assumes a uniform allele emergence probability by default. More on how to set user-defined mutation parameters when prior information is available.
By default, Phylotime uses future to parallel pairwise distance estimations, use plan to set up parallel method and number of workers. Alternatively, set 'parallel=F'.
In this example, barcoding activates at 0.6 days, and samples are collected at 15 days.
library(furrr) plan(multisession, workers = 8) tr = phylotime(chr_mat_imputed, t_total = 15. - 0.6)
If mutation rate or allele emergence probabilities are known a priori, they can be provided to phylotime
mut_rate_tb = read_csv("../data/example/mut_p_mut_rate.csv") mut_rate_tb
allele_prob_tb = read_csv("../data/example/mut_p_recur_vec.csv") allele_prob_tb
tr = phylotime(chr_mat_imputed, mut_rate = mut_rate_tb$rate, recur_vec = make_recur_vec_list(allele_prob_tb), t_total = 15. - 0.6)
Now to visualize the reconstructed tree with the lineage barcodes:
plot_barcodes(chr_mat, tr, show_column_names = F)
To plot the reconstrcuted tree only:
plot_tr(tr)
Next, to apply ICE-FASE, we need, in additional, the terminal cell state classifications. This is provided as a named vector, where the name correspond to cell names. We can visualzie the tree with terminal cell states:
cell_type_tb = read_csv("../data/example/cell_type.csv") sc_celltypes = cell_type_tb$type names(sc_celltypes) = cell_type_tb$cell plot_barcodes(chr_mat[, 1:50], tr, tip_celltype = sc_celltypes, show_column_names = F)
To run ICE-FASE: ('root_time' is the time until barcode activation.)
res = ice_fase(tr, sc_celltypes, total_time = 15 - 0.6, root_time = 0.6)
With ICE-FASE fitted, a number of progenitor states are identified as nodes in the fate map topology. We can plot the reconstructed fate map topology with commitment times now. Note that the inferred progenitor states (iP) are denoted as "Node-x" below.
plot_topology(gr = res$gr, total_time = 15, gr_node_time = res$gr_trans_time, type_col = res$col_pal, show_node_label = T)
To visualize the internal node state assignment in reconstructed phylogeny:
plot_tr(res$tr, node_types = res$tr_node_assign, type_col = res$col_pal)
To visualize distribution of Inferred Commitment Event (ICE) times:
plot_ice_times(res)
To visualize progenitor population sizes:
plot_node_sizes(res)
Finally, we can get a summary table for the progenitor state estimates:
output_estimates(res)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.