Introduction

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.

Clone ID with a clonal tree provided

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:

• A: a variant x cell matrix of integer counts, providing the number of reads supporting the alternative allele for each variant in each cell;
• D: a variant x cell matrix of integer counts, providing the total number of reads overlapping each variant in each cell;
• Config: a variant x clone “configuration” matrix of binary values providing the clone-variant configuration by indicating which mutations are expected to be present in which clones.

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.

Here, we demonstrate the use of cardelino to assign 77 cells to clones identified with Canopy using 112 somatic mutations.

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.

library(ggplot2)
library(cardelino)

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 cardelino:

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
## [1] "112 out of 112 SNPs passed."

Alternatively you can load the A and 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 Canopy:

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.
## [1] "Converged in 800 iterations."
## DIC: 1384.96 D_mean: 1230.68 D_post: 1203.9 logLik_var: 45.27

names(assignments)

##  [1] "theta0"         "theta1"         "theta0_all"     "theta1_all"    
##  [5] "element"        "logLik"         "prob_all"       "prob"          
##  [9] "prob_variant"   "relax_rate"     "Config_prob"    "Config_all"    
## [13] "relax_rate_all" "DIC"            "n_chain"

We can visualise the cell-clone assignment probabilities as a heatmap.

prob_heatmap(assignments$prob)

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.9999971
## 3 ERR2806036 clone1 0.9998761
## 4 ERR2806037 clone1 0.9996611
## 5 ERR2806038 clone1 0.9999999
## 6 ERR2806039 clone2 1.0000000

table(df$clone)

## 
##     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)

Clone ID without input clonal tree

As shown above, the Config can be updated by the observations from scRNA-seq data. A step further is to not including the Config entirely. This can be possible, as there could be no clonal tree available. In this case, you can use cardelino in its de-novo by set Config=NULL and set a number for n_clones. Note, by default we will keep the first clone as base clone, i.e., no mutations. You can turn it off by set keep_base_clone=FALSE.

{denovo-cell-assign} assignments <- clone_id(input_data$A, input_data$D, Config=NULL, n_clone = 3)

AD_file <- system.file("extdata", "passed_ad.mtx", package = "cardelino")
DP_file <- system.file("extdata", "passed_dp.mtx", package = "cardelino")
id_file <- system.file("extdata", "passed_variant_names.txt", 
                       package = "cardelino")

AD <- Matrix::readMM(AD_file)
DP <- Matrix::readMM(DP_file)
var_ids <- read.table(id_file, )
rownames(AD) <- rownames(DP) <- var_ids[, 1]
colnames(AD) <- colnames(DP) <- paste0('Cell', seq(ncol(DP)))

AF_mt <- as.matrix(AD / DP)
pheatmap::pheatmap(AF_mt)

set.seed(7)
assign_mtClones <- clone_id(AD, DP, Config=NULL, n_clone = 3, 
                            keep_base_clone=FALSE)

## Config is NULL: de-novo mode is in use.

## 100% converged.
## [1] "Converged in 5000 iterations."
## DIC: NaN D_mean: Inf D_post: Inf logLik_var: NaN

Config_mt <- assign_mtClones$Config_prob
Config_mt[Config_mt >= 0.5] = 1
Config_mt[Config_mt <  0.5] = 0
cardelino::vc_heatmap(AF_mt, assign_mtClones$prob, Config_mt, show_legend=TRUE)