Trouble running Lamian using UMAP embeddings

[DavidJDeGeorge]

Firstly, thank you for putting together such a useful tool!

Apologies in advance if this is a naive issue as I am relatively new to R and scRNAseq analysis. I am trying to use Lamian to compare pseudotime on a single cell dataset derived from two mouse models as part of my PhD. I have used Harmony to integrate my separate samples and am working with a Seurat object.

Following the online examples for Lamian, I have managed to put together the following code that allows me to create the equivalent of "res" in the online example and to evaluate branch uncertainty when I use the Harmony embeddings - this is as far as I have gone so far so this seems to work.

umap_coords <- Embeddings(seurat_obj, reduction = "umap")
umap_matrix <- as.matrix(umap_coords)

log_norm_data <- as.matrix(seurat_obj@assays[["RNA"]]@data)

metadata <- seurat_obj@meta.data
genotype_info <- metadata[["genotype"]]

cellanno <- data.frame(
  cell = rownames(metadata),
  sample = genotype_info
)

lamian_input <- list(
  pca = umap_matrix,
  expression = log_norm_data,
  cellanno = cellanno
)

set.seed(12345)
res = infer_tree_structure(pca = lamian_input[['pca']],
                           expression = lamian_input[['expression']],
                           cellanno = lamian_input[['cellanno']],
                           origin.marker = c('Myb'),
                           number.cluster = 10,
                           xlab='Principal component 1',
                           ylab = 'Principal component 2')

Lamian::plotmclust(res, cell_point_size = 0.5,
                   x.lab = 'Principal component 1',
                   y.lab = 'Principal component 2')

result <- evaluate_uncertainty(res, n.permute=100)

However, I want to visualise branches in the UMAP space that I have used for the rest of my analysis. When I change my input to instead include the UMAP embeddings (I have done this by simply replacing "harmony" in the above code with "umap" for now, I am able to generate the "res" equivalent and generate a plot with branches overlaid on my UMAP space. Except now when I try to evaluate uncertainty I receive the following error message:

Error in branchcomb[tmpcombid, 2] : subscript out of bounds

As far as I can tell, the structures of the two "res" objects (when either using "umap" embeddings or "harmony" embeddings are quite similar). I am possibly being too naive and it might just not be possible to run this on UMAP embeddings, however, it would be good to know a way of showing the branching in the context of something more relatable to the rest of my analysis, regardless.

Thank you for any assistance or clarifications that can be provided.

[Alexis-Varin]

Hi, I am using UMAP embeddings and passing evaluate_uncertainty() fine so you can definitely use them. I am not quite sure what you mean by "harmony embeddings", harmony isn't a low dimensional space like UMAP, you only have 3 choices, PCA, t-SNE, or UMAP. Embeddings(seurat_obj, reduction = "harmony") will work because Seurat 5 now stores the integration results as corrected PCA coordinates instead of the integrated Assay before but these reductions should not be used for any downstream analysis, only for UMAP calculation.

Can you do c(ncol(log_norm_data), length(umap_matrix)/2, length(umap_coords)/2, length(cellanno$cell))

You should obtain 4 times the same number.

[jingsday]

Hello Alexis,

I am encountering similar issues as I am using scVI embeddings as indicated in the paper under result section: The input for Lamian includes (1) a low-dimensional representation of cells, such as principal components (PCs) or other low-dimensional embeddings of the scRNA-seq data from multiple samples that have been harmonized into a common space using methods such as Seurat[32], Harmony[33] or scVI.

Could you help confirm which is the correct input for evaluate_uncertainty()? From this, a design output will be generated for further trajectory differential tests. Would it be possible to bypass this to get design calculated? Thank you for you attention and help.

Regards, Jing

Error:

result <- evaluate_uncertainty(res, n.permute=100) Error in [<-(*tmp*, rownames(ctcomp.logit), colnames(ctcomp.logit), : subscript out of bounds

Part of my res resulting from infer_tree_structure(): str(res) List of 11 $ pcareduceres: num [1:171836, 1:5] -0.197 -0.294 -0.725 0.4 -0.356 ... ..- attr(, "dimnames")=List of 2 .. ..$ : chr [1:171836] "AAAGTAGGTCCGTTAA.7.NEC-2" "AAGGAGCCATGTTCCC.7.NEC-2" "ACACCAAGTCCATCCT.7.NEC-2" "ACACCGGTCGAGCCCA.7.NEC-2" ... .. ..$ : chr [1:5] "scvi_1" "scvi_2" "scvi_3" "scvi_4" ... $ MSTtree :Class 'igraph' hidden list of 10 ..$ : num 6 ..$ : logi FALSE ..$ : num [1:5] 1 2 4 5 5 ..$ : num [1:5] 0 1 1 2 3 ..$ : NULL ..$ : NULL ..$ : NULL ..$ : NULL ..$ :List of 4 .. ..$ : num [1:3] 1 0 1 .. ..$ : Named list() .. ..$ :List of 1 .. .. ..$ name: chr [1:6] "1" "2" "3" "4" ... .. ..$ :List of 1 .. .. ..$ weight: num [1:5] 1.035 0.764 0.826 0.661 0.704 ..$ :<environment: 0x5af89f611f48> $ clusterid : Named int [1:171836] 3 3 1 3 2 1 1 2 1 5 ... ..- attr(, "names")= chr [1:171836] "AAAGTAGGTCCGTTAA.7.NEC-2" "AAGGAGCCATGTTCCC.7.NEC-2" "ACACCAAGTCCATCCT.7.NEC-2" "ACACCGGTCGAGCCCA.7.NEC-2" ... $ clucenter : num [1:6, 1:5] -0.6679 -0.0795 0.0754 -0.2952 0.4157 ... $ pseudotime : Named int [1:288244] 1 2 3 4 5 6 7 8 9 10 ... ..- attr(, "names")= chr [1:288244] "SRR10784712-13" "SRR10788545-13" "SRR10791476-13" "SRR10782733-13" ... $ branch :List of 3 ..$ : int [1:2] 2 1 ..$ : int [1:4] 4 6 3 2 ..$ : int [1:2] 5 2 $ js.cut : Named num [1:3] 0.169 0.609 0.157 ..- attr(, "names")= chr [1:3] "99%" "99%" "99%" $ oc.cut : Named num [1:3] 0.289 0.757 0.271 ..- attr(, "names")= chr [1:3] "99%" "99%" "99%" $ pca : num [1:171836, 1:5] -0.197 -0.294 -0.725 0.4 -0.356 ... ..- attr(, "dimnames")=List of 2 .. ..$ : chr [1:171836] "AAAGTAGGTCCGTTAA.7.NEC-2" "AAGGAGCCATGTTCCC.7.NEC-2" "ACACCAAGTCCATCCT.7.NEC-2" "ACACCGGTCGAGCCCA.7.NEC-2" ... .. ..$ : chr [1:5] "scvi_1" "scvi_2" "scvi_3" "scvi_4" ... $ order :List of 2 ..$ backbone 4,6,3,2,1: chr [1:144886] "SRR10784712-13" "SRR10788545-13" "SRR10791476-13" "SRR10782733-13" ... ..$ branch: 4,6,3,2,5 : chr [1:143358] "SRR10784712-13" "SRR10788545-13" "SRR10791476-13" "SRR10782733-13" ... $ allsample : Factor w/ 13 levels "cell type 1",..: 12 12 9 12 12 9 9 12 9 11 ... ..- attr(*, "names")= chr [1:171836] "AAAGTAGGTCCGTTAA.7.NEC-2" "AAGGAGCCATGTTCCC.7.NEC-2" "ACACCAAGTCCATCCT.7.NEC-2" "ACACCGGTCGAGCCCA.7.NEC-2" ...

[Alexis-Varin]

Hello @jingsday, Just like David using Harmony, I don't think scVI has its own low dimensional space embeddings, all analysis pipelines I have seen just now with a quick check (I am really not familiar with scVI or Harmony) use a UMAP calculation step, on the scVI results, to visualize cell clusters. So, I don't know what these 5 scVI dimensions represent, but I am not sure these are coordinates and represent something that can be input for pseudotime inference (even though you do get results...), or if these could be PCA but named after scVI (which would be misleading).

Does it work fine if you plot the results ? Like with :

Lamian::plotmclust(res, cell_point_size = 0.5,
                   x.lab = 'Principal component 1',
                   y.lab = 'Principal component 2')