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.

Filter loci

We remove A>G and T>C mutations due to RNA editing. Since the data didn’t have annotated REF or ALT alleles, we use the major allele as the reference.

cac <- collapse_rna_editing(cac, reference="major")
cac
#> CellAlleleCounts: 1887 loci, 2369 cells, 8936 features

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"
#> [5] "chr1_10376834"  "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,
    n_loci = 30, # default
    verbose = 1
)
#> [07:22:28] Step 1: Fitting FMM with 20 clusters
#> [07:22:28] Computing hclust initialization
#> [07:22:31] Finished hclust initialization
#> [07:22:31] Initializing FMM
#> [07:22:31] Finished initializing FMM
#> [07:22:31] Starting EM (1887 loci, 100 iterations max)
#> [07:22:33] Building model object
#> [07:22:33] Refitting FMM on 30 loci (round 1 of 5)
#> [07:22:33] Initializing FMM
#> [07:22:33] Finished initializing FMM
#> [07:22:33] Starting EM (30 loci, 100 iterations max)
#> [07:22:34] Building model object
#> [07:22:34] Refitting FMM on 30 loci (round 2 of 5)
#> [07:22:34] Initializing FMM
#> [07:22:34] Finished initializing FMM
#> [07:22:34] Starting EM (30 loci, 100 iterations max)
#> [07:22:34] Building model object
#> [07:22:34] Refitting FMM on 30 loci (round 3 of 5)
#> [07:22:34] Initializing FMM
#> [07:22:34] Finished initializing FMM
#> [07:22:34] Starting EM (30 loci, 100 iterations max)
#> [07:22:34] Building model object
#> [07:22:34] Refitting FMM on 30 loci (round 4 of 5)
#> [07:22:34] Initializing FMM
#> [07:22:34] Finished initializing FMM
#> [07:22:34] Starting EM (30 loci, 100 iterations max)
#> [07:22:35] Building model object
#> [07:22:35] Refitting FMM on 30 loci (round 5 of 5)
#> [07:22:35] Initializing FMM
#> [07:22:35] Finished initializing FMM
#> [07:22:35] Starting EM (30 loci, 100 iterations max)
#> [07:22:35] Building model object
#> [07:22:35] Selected 30 informative loci
#> [07:22:35] Step 3: CLR-PCA on topic frequencies
#> [07:22:35] Using 7 PCs
#> [07:22:35] Step 4: Predicting frequencies (pseudocount = 1)
#> [07:22:35] Step 5: PCA: reducing 68 features to 7 PCs
#> [07:22:35] 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 adjusted 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 consider rerunning the pipeline with manual settings.

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.

The default uses 30 loci. But we can plot the mutualinfo scores to see if there is a better cutoff.

plot(sort(fit$mutualinfo, decreasing=TRUE)[1:100])

# plot without first 2 loci which dominate
plot(sort(fit$mutualinfo, decreasing=TRUE)[3:100])

Number of principal components

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

By default, we use the Marchenko-Pasteur 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.

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 slow intermediate steps that are the same between the 2 runs (namely, the initial clustering procedure). We demonstrate here, rerunning with fewer loci and PCs.

fit2 <- fitCellDistancePipeline(
    cac,
    n_loci = 20,
    n_pcs = 5,
    init = fit,
    verbose=1
)
#> [07:22:43] Step 1: Using FMM from previous pipeline result
#> [07:22:43] Refitting FMM on 20 loci (round 1 of 5)
#> [07:22:43] Initializing FMM
#> [07:22:43] Finished initializing FMM
#> [07:22:43] Starting EM (20 loci, 100 iterations max)
#> [07:22:44] Building model object
#> [07:22:44] Refitting FMM on 20 loci (round 2 of 5)
#> [07:22:44] Initializing FMM
#> [07:22:44] Finished initializing FMM
#> [07:22:44] Starting EM (20 loci, 100 iterations max)
#> [07:22:44] Building model object
#> [07:22:44] Refitting FMM on 20 loci (round 3 of 5)
#> [07:22:44] Initializing FMM
#> [07:22:44] Finished initializing FMM
#> [07:22:44] Starting EM (20 loci, 100 iterations max)
#> [07:22:45] Building model object
#> [07:22:45] Refitting FMM on 20 loci (round 4 of 5)
#> [07:22:45] Initializing FMM
#> [07:22:45] Finished initializing FMM
#> [07:22:45] Starting EM (20 loci, 100 iterations max)
#> [07:22:45] Building model object
#> [07:22:45] Refitting FMM on 20 loci (round 5 of 5)
#> [07:22:45] Initializing FMM
#> [07:22:45] Finished initializing FMM
#> [07:22:45] Starting EM (20 loci, 100 iterations max)
#> [07:22:46] Building model object
#> [07:22:46] Selected 20 informative loci
#> [07:22:46] Step 3: CLR-PCA on topic frequencies
#> [07:22:46] Using 5 PCs
#> [07:22:46] Step 4: Predicting frequencies (pseudocount = 1)
#> [07:22:46] Step 5: PCA: reducing 42 features to 5 PCs
#> [07:22:46] Step 6: Computing pairwise distances

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(
    fit$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.

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(fit$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.

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(fit$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.

sclanternR vignette

Jack Kamm

Revised: 2026-10-24