Example dataset

library(magrittr)
library(ggplot2)

library(ggtree)

library(sclanternR)

We illustrate sclanternR on a subset of the HAP1 cell-line data from our preprint, consisting of the top 2000 loci from the full analysis in our preprint.

We start with a matrix of allele counts in each locus in each cell:

dim(hap1_subset)
#> [1] 9168 2369

hap1_subset[1:5,1:5]
#> 5 x 5 sparse Matrix of class "dgRMatrix"
#>                     D_AACCTAGACGCTGCTG D_AACCTAGACTCACCCT D_AACCTAGACTCTCCAC
#> chr1_100084142_T                     .                  .                  .
#> chr1_100084142_TA                    .                  .                  .
#> chr1_100084142_TAA                   .                  .                  .
#> chr1_100084142_TAAA                  .                  1                  .
#> chr1_100187700_A                     .                  .                  .
#>                     D_AACGACTTGACTTGAC D_AACGACTTGGACTCAT
#> chr1_100084142_T                     .                  .
#> chr1_100084142_TA                    .                  .
#> chr1_100084142_TAA                   .                  .
#> chr1_100084142_TAAA                  .                  .
#> chr1_100187700_A                     .                  .

The rows of this matrix are are locus-allele pairs, and the columns are cells. The rownames should follow the format CHROM_POS_ALLELE, e.g. chr1_100084142_TAAA.

Construct cell metadata

For our example dataset, the cells (column) names contain the ground truth clone and the cell barcode. We extract this metadata here so we can make plots with it later.

cell_meta <- stringr::str_match(colnames(hap1_subset), "(.*)_(.*)")
colnames(cell_meta) <- c("cell_id", "clone", "cell_barcode")
cell_meta %<>% as.data.frame()
rownames(cell_meta) <- cell_meta$cell_id

head(cell_meta)
#>                               cell_id clone     cell_barcode
#> D_AACCTAGACGCTGCTG D_AACCTAGACGCTGCTG     D AACCTAGACGCTGCTG
#> D_AACCTAGACTCACCCT D_AACCTAGACTCACCCT     D AACCTAGACTCACCCT
#> D_AACCTAGACTCTCCAC D_AACCTAGACTCTCCAC     D AACCTAGACTCTCCAC
#> D_AACGACTTGACTTGAC D_AACGACTTGACTTGAC     D AACGACTTGACTTGAC
#> D_AACGACTTGGACTCAT D_AACGACTTGGACTCAT     D AACGACTTGGACTCAT
#> D_AACGACTTGTCACAGT D_AACGACTTGTCACAGT     D AACGACTTGTCACAGT

clone_colors <- pals::brewer.paired(12)[c(7,5,6,1,2)]
names(clone_colors) <- sort(unique(cell_meta$clone))

Data object construction

We start by converting the count matrix into a CellAlleleCounts object:

cac <- CellAlleleCounts(hap1_subset)
cac
#> CellAlleleCounts: 2000 loci, 2369 cells, 9168 features

Alternatively, if using PacBio long-read data processed by our sclantern-nf pipeline, the CellAlleleCounts can be created by calling CellAlleleCounts_from_dataframe() on the long table of counts produced by that pipeline.

Data slicing

Loci can be sliced along rows and cells can be sliced along columns:

cac[1:5, 1:5]
#> CellAlleleCounts: 5 loci, 5 cells, 27 features

It can also be alternatively treated like a list; the ith element is a matrix of the counts at the ith locus:

head(names(cac))
#> [1] "chr1_100084142" "chr1_100187700" "chr1_100291929" "chr1_103551885" "chr1_10376834" 
#> [6] "chr1_10420389"

first <- cac[[1]] # same as cac[["chr1_100084142"]]

dim(first)
#> [1]    4 2369

first[, 1:5]
#> 4 x 5 sparse Matrix of class "dgRMatrix"
#>                     D_AACCTAGACGCTGCTG D_AACCTAGACTCACCCT D_AACCTAGACTCTCCAC
#> chr1_100084142_T                     .                  .                  .
#> chr1_100084142_TA                    .                  .                  .
#> chr1_100084142_TAA                   .                  .                  .
#> chr1_100084142_TAAA                  .                  1                  .
#>                     D_AACGACTTGACTTGAC D_AACGACTTGGACTCAT
#> chr1_100084142_T                     .                  .
#> chr1_100084142_TA                    .                  .
#> chr1_100084142_TAA                   .                  .
#> chr1_100084142_TAAA                  .                  .

Run analysis pipeline with default parameters

The sclanternR procedure involves several steps, consisting of an initial clustering, selection of informative loci, smoothing of allele frequencies, and dimension reduction. The function fitCellDistancePipeline() runs each of these steps with automatic defaults.

fit <- fitCellDistancePipeline(
    cac, verbose=1
)
#> [07:45:06] Step 1: Fitting FMM with 20 clusters
#> [07:45:06] Fitting combined multinomial initialization
#> [07:45:06] Computing hclust initialization
#> [07:45:07] Finished hclust initialization
#> [07:45:07] Initializing FMM
#> [07:45:07] Finished initializing FMM
#> [07:45:07] Starting EM (2000 loci, 100 iterations max)
#> [07:45:08] Building model object
#> [07:45:08] Initializing FMM
#> [07:45:08] Finished initializing FMM
#> [07:45:08] Starting EM (2000 loci, 100 iterations max)
#> [07:45:09] Building model object
#> [07:45:09] Step 2: CLR-PCA on topic frequencies
#> [07:45:09] Using 5 PCs
#> [07:45:09] Step 3: Computing mutual information across 2000 loci
#> [07:45:09] Automatic MI cutoff (fdr): 18.1, 300 loci selected
#> [07:45:09] Selected 300 informative loci
#> [07:45:09] Step 4: Predicting frequencies (pseudocount = 1)
#> [07:45:10] Step 5: PCA: reducing 1345 features to 5 PCs
#> [07:45:11] Step 6: Computing pairwise distances

Examine PCA

We plot the PCA embedding produced by the pipeline.

# shuffle order of cells for plotting
set.seed(12345)
shuffle_order <- sample(1:ncol(cac))

ggplot(
    cbind(
        cell_meta,
        fit$cell_pca
    )[shuffle_order,],
    aes(x=PC1, y=PC2, color=clone)
) +
    geom_point(shape=1) +
    scale_color_manual(values=clone_colors) +
    theme_bw(base_size=16)

Tree building

For building distance-based phylogenies (e.g. with neighbor-joining), we recommend to use the squared Euclidean distances in the PCA space. These are also stored in the dist field of the fitted object.

tr <- ape::nj(fit$dist)

ggtree(tr, layout='ape') %<+%
    cell_meta +
    ggrastr::rasterise(geom_point(
        data = function(d) {
            d %>% dplyr::filter(isTip) %>% .[shuffle_order,]
        },
        #shape=1,
        aes(color = clone)
    ), dpi=300) +
    scale_color_manual(values=clone_colors) +
    theme(legend.position='bottom')

Rerunning with better parameters

fitCellDistancePipeline() has several tuning parameters, and it tries to set reasonable defaults automatically. However, the automatic heuristics are imperfect, and it’s recommended to visualize these cutoffs and set them manually.

Number of loci

The first cutoff to consider is the number of informative loci to retain. We use a scaled version of the mutual information (G-statistic) to rank the loci.

For the automatic cutoff, we test for the loci being independent of the clusters, and select loci at 5% FDR, with bounds of at least 10 loci or at most 300 loci selected.

Below, we plot the scaled mutual information scores and the cutoff selected. The maximum cutoff of 300 was chosen; however this was likely too high, as we can visually see a sharp drop at between 10-30 loci. Also, as a rule of thumb, we’ve found that selecting between 30-100 loci tends to give the best results.

plot(fit$mi_cutoff)

Number of principal components

The second cutoff to consider is the number of PCs to keep.

By default, we use the Gavish-Donoho procedure to select the number of PCs. But it’s recommended to check the scree plot to see if there is a more obvious inflection point visually. In this case, there wasn’t a more obvious cutoff, and we stick with the default selected 5 PCs.

plot(fit$scree)

Rerunning the pipeline

The pipeline can be rerun with the manually selected cutoffs. We can pass the previous fit in order to skip intermediate steps that are the same between the 2 runs (namely, the initial clustering procedure).

fit2 <- fitCellDistancePipeline(
    cac,
    n_loci = 30,
    n_pcs = 5,
    init = fit,
    verbose=1
)
#> [07:45:19] Step 1: Using FMM from previous pipeline result
#> [07:45:19] Step 2: CLR-PCA on topic frequencies
#> [07:45:19] Using 5 PCs
#> [07:45:19] Step 3: Computing mutual information across 2000 loci
#> [07:45:19] Selected 30 informative loci
#> [07:45:19] Step 4: Predicting frequencies (pseudocount = 1)
#> [07:45:19] Step 5: PCA: reducing 64 features to 5 PCs
#> [07:45:19] Step 6: Computing pairwise distances

Plot updated PCA

We plot the updated PCA from the rerun fit.

ggplot(
    cbind(
        cell_meta,
        fit2$cell_pca
    )[shuffle_order,],
    aes(x=PC1, y=PC2, color=clone)
) +
    geom_point(shape=1) +
    scale_color_manual(values=clone_colors) +
    theme_bw(base_size=16)

Plot updated tree

We recompute the neighbor joining tree and plot it.

tr <- ape::nj(fit2$dist)

ggtree(tr, layout='ape') %<+%
    cell_meta +
    ggrastr::rasterise(geom_point(
        data = function(d) {
            d %>% dplyr::filter(isTip) %>% .[shuffle_order,]
        },
        #shape=1,
        aes(color = clone)
    ), dpi=300) +
    scale_color_manual(values=clone_colors) +
    theme(legend.position='bottom')

Plotting heatmap of informative loci

It is also useful to plot the estimated cell-level allele frequencies at the selected loci. The fit object contains these frequencies in the freq_list field, which can be plotted with heatmapLocusMatrixList().

heatmapLocusMatrixList(
    fit2$freq_list,
    cell_meta[, 'clone', drop=FALSE],
    cell_meta_colors=list(
        clone=clone_colors
    )
)
#> `use_raster` is automatically set to TRUE for a matrix with more than 2000 rows.
#> You can control `use_raster` argument by explicitly setting TRUE/FALSE to it.
#> 
#> Set `ht_opt$message = FALSE` to turn off this message.
#> 'magick' package is suggested to install to give better rasterization.
#> 
#> Set `ht_opt$message = FALSE` to turn off this message.

By default the cells are sorted by a hierarchical clustering, but the order can also be manually specified by setting cluster_rows=FALSE. Below, we sort the cells by ground truth clone.

o <- order(cell_meta$clone)

heatmapLocusMatrixList(
    lapply(fit2$freq_list, function(x) x[,o]),
    cell_meta[, 'clone', drop=FALSE],
    cell_meta_colors=list(
        clone=clone_colors
    ),
    cluster_rows = FALSE
)
#> `use_raster` is automatically set to TRUE for a matrix with more than 2000 rows.
#> You can control `use_raster` argument by explicitly setting TRUE/FALSE to it.
#> 
#> Set `ht_opt$message = FALSE` to turn off this message.
#> 'magick' package is suggested to install to give better rasterization.
#> 
#> Set `ht_opt$message = FALSE` to turn off this message.

heatmapLocusMatrixList() can also be used to plot other list of matrices. Below, we demonstrate plotting the log10-scaled counts:

log_cnt_list <- lapply(cac[names(fit2$freq_list)],
                       function(x) log10(as.matrix(x) + 1))

heatmapLocusMatrixList(
    log_cnt_list,
    cell_meta[, 'clone', drop=FALSE],
    cell_meta_colors=list(
        clone=clone_colors
    ),
    cluster_rows = FALSE,
    colors = rev(pals::cubehelix(10))
)
#> `use_raster` is automatically set to TRUE for a matrix with more than 2000 rows.
#> You can control `use_raster` argument by explicitly setting TRUE/FALSE to it.
#> 
#> Set `ht_opt$message = FALSE` to turn off this message.
#> 'magick' package is suggested to install to give better rasterization.
#> 
#> Set `ht_opt$message = FALSE` to turn off this message.

sclanternR vignette

Jack Kamm

Revised: 2026-08-21