This document gives an introduction to and overview of inferring the clone identity of cells using the cardelino package using a given clonal structure.
cardelino can use variant information extracted from single-cell RNA-seq reads to probabilistically assign single-cell transcriptomes to individual clones.
Briefly, cardelino is based on a Bayesian mixture model with a beta-binomial error model to account for sequencing errors as well as a gene-specific model for allelic imbalance between haplotypes and associated bias in variant detection. Bayesian inference allows the model to account for uncertainty in model parameters and cell assignments.
We assume that clones are tagged by somatic mutations, and that these mutations are known (e.g. from exome sequencing or equivalent). Given a set of known mutations, these sites can be interrogated in scRNA-seq reads to obtain evidence for the presence or absence of each mutation in each cell. As input, the model requires the count of reads supporting the alternative (mutant) allele at each mutation site, the total number of reads overlapping the mutation site (“coverage”).
Typically, coverage of somatic mutations in scRNA-seq data is very sparse (most mutation sites in a given cell have no read coverage), but the cardelino model accounts for this sparsity and aggregates information across all available mutation sites to infer clonal identity.
In many clone ID scenarios, a clonal tree is known. That is, we have been able to infer the clones present in the sampled cell population, for example using bulk or single-cell DNA-seq data, and we know which mutations are expected to be present in which clones.
To infer the clonal identity of cells when a clonal tree is provided, cardelino requires the following input data:
The configuration matrix, Config, can be provided by other tools used to infer the clonal structure of the cell population. For example, the package Canopy can be used to infer a clonal tree from DNA-seq data and the “Z” element of its output is the configuration matrix.
We load the package and the example clone ID dataset distributed with the package in VCF (variant call format) format, which is mostly used for storing genotype data. Here, the
cellSNP.cells.vcf.gz is computed by using cellsnp-lite on a list pre-identified somatic variants from bulk WES.
There are many possible ways to extract the data required by
cardelino from a VCF file, here we show just one approach using convenience functions in
vcf_file <- system.file("extdata", "cellSNP.cells.vcf.gz", package = "cardelino") input_data <- load_cellSNP_vcf(vcf_file)
## Scanning file to determine attributes. ## File attributes: ## meta lines: 37 ## header_line: 38 ## variant count: 112 ## column count: 86 ## Meta line 37 read in. ## All meta lines processed. ## gt matrix initialized. ## Character matrix gt created. ## Character matrix gt rows: 112 ## Character matrix gt cols: 86 ## skip: 0 ## nrows: 112 ## row_num: 0 ## Processed variant: 112 ## All variants processed ##  "112 out of 112 SNPs passed."
Alternatively you can load the
D matrices yourself and combine them into a list, for example
input_data = list('A' = A, 'D' = D).
We can visualize the allele frequency of the mutation allele. As expected, the majority of entries are missing (in grey) due to the high sparsity in scRNA-seq data. For the same reason, even for the non-missing entries, the estimate of allele frequency is of high uncertainty. For this reason, it is crucial to probabilistic clustering with accounting the uncertainty, ideally with guide clonal tree structure from external data.
AF <- as.matrix(input_data$A / input_data$D) p = pheatmap::pheatmap(AF, cluster_rows=FALSE, cluster_cols=FALSE, show_rownames = TRUE, show_colnames = TRUE, labels_row='77 cells', labels_col='112 SNVs', angle_col=0, angle_row=0)
Next, we will load the Canopy tree results for the same individual. The clonal tree inferred by Canopy for this donor consists of three clones, including a “base” clone (“clone1”) that has no subclonal somatic mutations present.
canopy_res <- readRDS(system.file("extdata", "canopy_results.coveraged.rds", package = "cardelino")) Config <- canopy_res$tree$Z
Be careful to ensure that the same variant IDs are used in both data sources.
rownames(Config) <- gsub(":", "_", rownames(Config))
We can visualize the clonal tree structure obtained from
plot_tree(canopy_res$tree, orient = "v")
The included dataset contains the A and D matrices, so combined with the Canopy tree object provided, we have the necessary input to probabilistically assign cells to clones. Note,
min_iter = 800, max_iter = 1200 is used only for quick illustration. Please remove them for the default values or set higher number of iterations to ensure convergence of the Gibbs sampling. Convergence is checked automatically in
clone_id(), using the Geweke z-statistic.
set.seed(7) assignments <- clone_id(input_data$A, input_data$D, Config = Config, min_iter = 800, max_iter = 1200)
## 100% converged. ##  "Converged in 800 iterations." ## DIC: 1398.26 D_mean: 1232.3 D_post: 1201.67 logLik_var: 49.15
##  "theta0" "theta1" "theta0_all" "theta1_all" ##  "element" "logLik" "prob_all" "prob" ##  "prob_variant" "relax_rate" "Config_prob" "Config_all" ##  "relax_rate_all" "DIC" "n_chain"
We can visualise the cell-clone assignment probabilities as a heatmap.
We recommend assigning a cell to the highest-probability clone if the highest posterior probability is greater than 0.5 and leaving cells “unassigned” if they do not reach this threshold. The
assign_cells_to_clones function conveniently assigns cells to clones based on a threshold and returns a data.frame with the cell-clone assignments.
df <- assign_cells_to_clones(assignments$prob) head(df)
## cell clone prob_max ## 1 ERR2806034 clone2 1.0000000 ## 2 ERR2806035 clone1 0.9999962 ## 3 ERR2806036 clone1 0.9998902 ## 4 ERR2806037 clone1 0.9996968 ## 5 ERR2806038 clone1 0.9999999 ## 6 ERR2806039 clone2 1.0000000
## ## clone1 clone2 clone3 unassigned ## 44 24 7 2
Also, Cardelino will update the guide clonal tree Config matrix (as a prior) and return a posterior estimate. In the figure below, negative value means the probability of a certain variant presents in a certain clone is reduced in posterior compared to prior (i.e., the input Config). Vice verse for the positive values.
heat_matrix(t(assignments$Config_prob - Config)) + scale_fill_gradient2() + ggtitle('Changes of clonal Config') + labs(x='Clones', y='112 SNVs') + theme(axis.text.y = element_blank(), legend.position = "right")
Finally, we can visualize the results cell assignment and updated mutations clonal configuration at the raw allele frequency matrix:
AF <- as.matrix(input_data$A / input_data$D) cardelino::vc_heatmap(AF, assignments$prob, Config, show_legend=TRUE)