knitr::opts_chunk$set(echo = TRUE)
library(qfm)

Read example data

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

Imputing missing data

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)

Running Phylotime

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)

Running ICE-FASE

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)

ICE-FASE results

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)


Kalhor-Lab/QFM documentation built on Nov. 25, 2024, 10:18 p.m.