Introduction

Latent Embedding Multivariate Regression (LEMUR)

lemur is a first-line tool for the analysis of multi-condition single-cell data. If you have collected a single-cell RNA-seq dataset from more than one condition (e.g., different experimental treatments, different patient samples), lemur predicts for each cell and gene what the expression data would be in all of the other conditions. lemur then uses these predictions to find groups of cells that show similar patterns of differential expression across the conditions. We call these groups “neighbourhoods”. You can subsequently overlap these neighbourhoods with pre-existing cell type annotations, or in fact use the discovered neighbourhoods to refine the pre-existing cell type and state annotation for your data.

The differential expression results are statistically validated using a pseudo-bulk differential expression test on hold-out data using glmGamPoi or edgeR.

lemur implements a novel framework to disentangle the effects of known covariates, latent cell states, and their interactions. At the core, is a combination of matrix factorization and regression analysis implemented as geodesic regression on Grassmann manifolds. We call this latent embedding multivariate regression. For more details see our paper.

Schematic of the matrix decomposition at the core of LEMUR

Schematic of the matrix decomposition at the core of LEMUR

Installation

You can install lemur directly from Bioconductor. Just paste the following snippet into your R console:

if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
BiocManager::install("lemur")

This will install the latest release version. Alternatively, you can install the development version from GitHub using devtools:

devtools::install_github("const-ae/lemur")

A note to users

We continue working on improvements to this package. We are delighted if you decide to try out the package. Please use the Bioconductor forum or open an issue in the GitHub repository if you think you found a bug, have an idea for a cool feature, or have any questions about how LEMUR works.

Overview

A basic lemur workflow is as easy as the following.

# ... sce is a SingleCellExperiment object with your data 
fit <- lemur(sce, design = ~ patient_id + condition, n_embedding = 15)
fit <- align_harmony(fit)   # This step is optional
fit <- compute_contrasts(fit, contrast = cond(condition = "ctrl") - cond(condition = "panobinostat"))
nei <- find_de_neighborhoods(fit, group_by = vars(patient_id, condition))

We will now go through these steps one by one.

A worked through example

We demonstrate lemur using data by Zhao et al. (2021). The data consist of tumour biopsies from five glioblastomas, each of which was treated with the drug panobinostat and with a control. Accordingly, we look at ten samples in a paired experimental design.

We start by loading required packages1.

library("tidyverse")
library("SingleCellExperiment")
library("lemur")
set.seed(42)
data("glioblastoma_example_data", package = "lemur")
glioblastoma_example_data
class: SingleCellExperiment 
dim: 300 5000 
metadata(0):
assays(2): counts logcounts
rownames(300): ENSG00000210082 ENSG00000118785 ... ENSG00000167468
  ENSG00000139289
rowData names(6): gene_id symbol ... strand. source
colnames(5000): CGCCAGAGCGCA AGCTTTACTGCG ... TGAACAGTGCGT TGACCGGAATGC
colData names(10): patient_id treatment_id ... sample_id id
reducedDimNames(0):
mainExpName: NULL
altExpNames(0):

As you can see above, the lemur package ships with a reduced-size version of the glioblastoma data, with 300 genes and 5000 cells. A UMAP visualisation indicates that the cells’ data separate them by the covariates patient_id and condition.

orig_umap <- uwot::umap(as.matrix(t(logcounts(glioblastoma_example_data))))

as_tibble(colData(glioblastoma_example_data)) |>
  mutate(umap = orig_umap) |>
  ggplot(aes(x = umap[,1], y = umap[,2])) +
    geom_point(aes(colour = patient_id, shape = condition), size = 0.5) +
    labs(title = "UMAP of logcounts") + coord_fixed()
Figure 1: UMAP of glioblastoma_example_data prior to between-condition alignment.

We fit the LEMUR model by calling the function lemur. We provide the experimental design using a formula. The elements of the formula refer to columns of the colData of the SingleCellExperiment object.

We need to set the number of latent dimensions, n_embedding, which has a similar interpretation as the number of dimensions in a principal component analysis (PCA). Its choice is a trade-off between data set size (larger data sets allow for higher values of n_embedding), biological complexity (samples with more different cell types require higher values of n_embedding), and compute resources (higher values of n_embedding require more memory and compute cycles). We recommend trying out a few different choices and verifying that the results do not substantially depend on it.

The test_fraction argument sets the fraction of cells that are held out from the fitting of the LEMUR model and will only be used later for the differential expression testing2.

fit <- lemur(glioblastoma_example_data, design = ~ patient_id + condition, 
             n_embedding = 15, test_fraction = 0.5)
fit
class: lemur_fit 
dim: 300 5000 
metadata(9): n_embedding design ... use_assay row_mask
assays(2): counts logcounts
rownames(300): ENSG00000210082 ENSG00000118785 ... ENSG00000167468
  ENSG00000139289
rowData names(6): gene_id symbol ... strand. source
colnames(5000): CGCCAGAGCGCA AGCTTTACTGCG ... TGAACAGTGCGT TGACCGGAATGC
colData names(10): patient_id treatment_id ... sample_id id
reducedDimNames(2): linearFit embedding
mainExpName: NULL
altExpNames(0):

The lemur function returns an object of class lemur_fit, which extends the SingleCellExperiment class. It supports subsetting and all the usual accessor methods (e.g., nrow, assay, colData, rowData). In addition, lemur overloads the $ operator to allow easy access to additional fields produced by the LEMUR model. For example, the low-dimensional embedding can be accessed using fit$embedding:

fit$embedding |> str()
 num [1:15, 1:5000] 8.077 0.208 2.317 -4.313 0.868 ...
 - attr(*, "dimnames")=List of 2
  ..$ : NULL
  ..$ : chr [1:5000] "CGCCAGAGCGCA" "AGCTTTACTGCG" "AGGATGACCGCA" "AGAACTATTTTT" ...

Optionally, we can do what is often called “batch correction”: tweak the data such that corresponding cells from different patients and conditions are closer together in the embedding space (fit$embedding). This can be done either using explicit guidance from a cell type annotation (through the align_by_grouping function) or using an automated alignment procedure, e.g., via the align_harmony function, which calls the eponymous method from the harmony package.

fit <- align_harmony(fit)

Figure 2 shows a UMAP visualisation of fit$embedding. This is similar to looking at an integrated PCA space in other single-cell analysis workflows.

umap <- uwot::umap(t(fit$embedding))

as_tibble(fit$colData) |>
  mutate(umap = umap) |>
  ggplot(aes(x = umap[,1], y = umap[,2])) +
    geom_point(aes(colour = patient_id), size = 0.5) +
    facet_wrap(vars(condition)) + coord_fixed()
Figure 2: UMAPs of fit$embedding. The points are shown separately for the two conditions, but reside in the same latent space. In contrast to Figure 1, the Harmony alignment procedure has been applied with the aim of making corresponding cells from the different patient samples and conditions lie next to each other.

Next, let’s predict the effect of the panobinostat treatment for each cell and each gene – even for the cells that were observed in the control condition. The compute_contrasts function takes a lemur_fit object and returns the object with a new slot (in SummarizedExperiment parlance: assay) called DE. This slot contains the predicted logarithmic fold changes between the two conditions specified in contrast. Note that lemur implements a special notation for contrasts. Instead of providing a contrast vector or design matrix column names, you provide for each condition the levels, and lemur automatically forms the contrast vector. This is intended to make the notation more readable.

fit <- compute_contrasts(fit, 
         contrast = cond(condition = "panobinostat") - cond(condition = "ctrl"))

We can pick a gene, say GAP43, which in our data is represented by its Ensembl gene ID ENSG00000172020, and show its differential expression pattern on the UMAP plot:

df <- tibble(umap = umap) |>
  mutate(de = assay(fit, "DE")["ENSG00000172020", ])
 
ggplot(df, aes(x = umap[,1], y = umap[,2])) +
    geom_point(aes(colour = de)) +
    scale_colour_gradient2(low = "#FFD800", high= "#0056B9") + coord_fixed()

ggplot(df, aes(x = de)) + geom_histogram(bins = 100)
(a) Superimposed on the same UMAP plot as in Figure 2
(b) The histogram shows that the values are negative for most cells.
Figure 3: Differential expression (log fold changes) of GAP43.

A gene that is down-regulated in most cells upon treatment is interesting, but perhaps even more interesting are genes that show differential expression only in certain subsets of the cells. We iterate through all genes and identify what we call neighbourhoods (sets of cells that are close together in latent space) that show a similar differential expression pattern (assay(fit, "DE")). The function find_de_neighborhoods then statistically screens the results of our search with pseudobulk differential expression testing. For that, it uses the test data (fit$test_data) that were put aside when we called lemur(). In addition, find_de_neighborhoods assesses if the difference between the conditions is significantly larger for the cells inside the neighbourhood than for the cells outside the neighbourhood (see columns starting with did, short for difference-in-difference).

The group_by argument determines how the pseudobulk samples are formed. It specifies the columns in fit$colData that are used to define a sample and is inspired by the group_by function in dplyr. Typically, you provide the covariates that were used for the experimental design. 

neighborhoods <- find_de_neighborhoods(fit, group_by = vars(patient_id, condition))

We print the results for the 5 neighbourhoods with the smallest p values. Here, we use left_join and relocate to add the gene symbol from the fit object.

as_tibble(neighborhoods) |>
  left_join(as_tibble(rowData(fit)[,1:2]), by = c("name" = "gene_id")) |>
  relocate(symbol, .before = "name") |>
  arrange(pval) |>
  head(5)
# A tibble: 5 × 14
  symbol name      neighborhood n_cells sel_statistic    pval adj_pval f_statistic   df1
  <chr>  <chr>     <I<list>>      <int>         <dbl>   <dbl>    <dbl>       <dbl> <int>
1 MT1X   ENSG0000… <chr>           2285         -43.2 1.22e-5  0.00365       124.      1
2 KNG1   ENSG0000… <chr>           4661          49.2 7.51e-5  0.0113         70.4     1
3 POLR2L ENSG0000… <chr>           3816          95.1 1.20e-4  0.0119         60.6     1
4 PMP2   ENSG0000… <chr>           4009         -74.3 1.87e-4  0.0119         52.5     1
5 NEAT1  ENSG0000… <chr>           3530          74.8 1.99e-4  0.0119         51.5     1
# ℹ 5 more variables: df2 <dbl>, lfc <dbl>, did_pval <dbl>, did_adj_pval <dbl>,
#   did_lfc <dbl>

Let’s look at SERPINE1, which encodes PAI-1 (Plasminogen Activator Inhibitor-1), an inhibitor of fibrinolysis that regulates blood clotting, cell migration, and tissue remodeling. We see that it is upregulated by upon panobinostat treatement in a subset of cells (blue). We picked this gene because it (1) had a significant change between the panobinostat and negative control conditions (adj_pval column) and (2) showed larger differential expression for the cells inside the neighbourhood than for the cells outside (did_lfc column).

sel_gene <- "ENSG00000106366" # SERPINE1

dplyr::filter(neighborhoods, name == sel_gene) |>
  dplyr::select(name, adj_pval, lfc, did_adj_pval, did_lfc)
             name  adj_pval       lfc did_adj_pval    did_lfc
1 ENSG00000106366 0.1036024 0.8207343    0.9927934 -0.4241759
p <- tibble(umap = umap) |>
  mutate(de = assay(fit, "DE")[sel_gene,]) |>
  ggplot(aes(x = umap[,1], y = umap[,2])) +
    geom_point(aes(colour = de)) +
    scale_colour_gradient2(low = "#FFD800", high= "#0056B9") +
    coord_fixed()
p
Figure 4: Differential expression of SERPINE1 superimposed on the same UMAP plot as in Figure 2.

To get a further overview, we can make a volcano plot of the differential expression results across all genes.

neighborhoods |>
  drop_na() |>
  ggplot(aes(x = lfc, y = -log10(pval))) + geom_point(aes(col = adj_pval < 0.1)) + 
    scale_colour_manual(values= c(`TRUE` = "orangered", `FALSE` = "thistle3"))
Figure 5: Volcano plot, each point corresponds to one of the genes in our dataset, and the neighbourhood that was found for it.

Using cell type annotation

The analyses up to here were conducted without using pre-existing cell type annotation. Often, such additional cell type information is available or can be obtained from the data by other means. For instance, here, we can distinguish the tumour cells from other, non-malignant cells using the fact that the tumour cells had a deletion of Chromosome 10 and a duplication of Chromosome 7. We build a simple classifier to distinguish the cells accordingly. (This is just to illustrate the process; for a real analysis, you would use more sophisticated methods, and a more detailed cell type classification.)

thresh = c(chr7 = 0.8, chr10 = 2.5)
tumour_label_df <- tibble(cell_id = colnames(fit),
       chr7_total_expr  = colMeans(logcounts(fit)[rowData(fit)$chromosome == "7",]),
       chr10_total_expr = colMeans(logcounts(fit)[rowData(fit)$chromosome == "10",])) |>
  mutate(is_tumour = chr7_total_expr > thresh["chr7"] & chr10_total_expr < thresh["chr10"])

ggplot(tumour_label_df, aes(x = chr10_total_expr, y = chr7_total_expr, colour = is_tumour)) +
  geom_point(size = 0.5) +
  geom_hline(yintercept = thresh["chr7"]) +
  geom_vline(xintercept = thresh["chr10"]) +
  scale_colour_manual(values= c(`TRUE` = "orangered", `FALSE` = "#404040"))
Figure 6: A simple gating strategy to putatively annotate the tumour cells in our data 
tibble(umap = umap) |>
  mutate(is_tumour = tumour_label_df$is_tumour) |>
  ggplot(aes(x = umap[,1], y = umap[,2], colour = is_tumour)) +
    geom_point(size = 0.5) +
    facet_wrap(vars(is_tumour)) + coord_fixed() + 
    scale_colour_manual(values= c(`TRUE` = "orangered", `FALSE` = "#404040"))
Figure 7: The tumour cells appear enriched in the larger of the blobs.

We can use this cell annotation to re-run the neighbourhood finding, now only within the tumour cells, to find tumour subpopulations.

tumour_fit <- fit[, tumour_label_df$is_tumour]
tum_nei <- find_de_neighborhoods(tumour_fit, group_by = vars(patient_id, condition), verbose = FALSE)

as_tibble(tum_nei) |>
  left_join(as_tibble(rowData(fit)[,1:2]), by = c("name" = "gene_id")) |>
  dplyr::relocate(symbol, .before = "name") |>
  filter(adj_pval < 0.1) |>
  arrange(did_pval)  |>
  dplyr::select(symbol, name, neighborhood, n_cells, adj_pval, lfc, did_pval, did_lfc) |>
  print(n = 10)
# A tibble: 24 × 8
   symbol    name            neighborhood  n_cells adj_pval    lfc did_pval did_lfc
   <chr>     <chr>           <I<list>>       <int>    <dbl>  <dbl>    <dbl>   <dbl>
 1 HLA-B     ENSG00000234745 <chr [618]>       618  0.0690  -0.819  0.00136   0.776
 2 HMGB1     ENSG00000189403 <chr [2,940]>    2940  0.0541  -0.773  0.0358    0.691
 3 TUBA1A    ENSG00000167552 <chr [1,521]>    1521  0.0265  -0.816  0.0438    0.409
 4 MT1X      ENSG00000187193 <chr [1,201]>    1201  0.0818   1.83   0.0546    1.30 
 5 NEAT1     ENSG00000245532 <chr [2,637]>    2637  0.00822  1.85   0.106    -0.688
 6 TNFRSF12A ENSG00000006327 <chr [1,603]>    1603  0.0369  -0.834  0.106     0.320
 7 NRCAM     ENSG00000091129 <chr [1,621]>    1621  0.0797  -1.49   0.128     0.788
 8 MT2A      ENSG00000125148 <chr [847]>       847  0.0332   1.14   0.155     0.505
 9 CALM1     ENSG00000198668 <chr [2,104]>    2104  0.0369   0.841  0.193    -0.323
10 CD14      ENSG00000170458 <chr [1,559]>    1559  0.0665  -3.69   0.232     2.20 
# ℹ 14 more rows

Let’s, for example, explore the results for CXCL8. We see in Figure 8 that it is lower expressed in our annotated tumour cells than in the others, and there is a subpopulation of our annotated tumour cells in which the panobinostat treatment upregulates it. Caveat: this analysis is conducted on a subset of the data, using a primitive method to annotate the cells into just two types (tumour and non-tumour). For real analyses, you will be more careful. Nevertheless, this is the basic idea for how lemur can be used to find cell type subpopulations that differ in how they respond to a treatment.

sel_gene <- "ENSG00000169429" # CXCL8

as_tibble(colData(fit)) |>
  mutate(expr = assay(fit, "logcounts")[sel_gene, ]) |>
  mutate(is_tumour = tumour_label_df$is_tumour) |>
  mutate(in_neighborhood = id %in% filter(tum_nei, name == sel_gene)$neighborhood[[1]]) |>
  ggplot(aes(x = condition, y = ifelse(expr < 10, expr, 10))) +
    geom_jitter(size = pi/4, stroke = 0) +
    geom_point(data = . %>% summarize(expr = mean(expr), .by = c(condition, patient_id, is_tumour, in_neighborhood)),
               aes(colour = patient_id), size = 2) +
    stat_summary(fun.data = mean_se, geom = "crossbar", colour = "red") +
    facet_wrap(vars(is_tumour, in_neighborhood), labeller = label_both)
Figure 8: CXCL8 expression in different cell subsets

FAQ

I have already integrated my data using Harmony / MNN / Seurat. Can I call lemur directly with the aligned data?

No. You need to call lemur with the unaligned data so that it can learn how much the expression of each gene changes between conditions.

Can I call lemur with sctransformed instead of log-transformed data?

Yes. You can call lemur with any variance stabilized count matrix. Based on a previous project, I recommend to use log-transformation, but other methods can work just as well.

My data appears less integrated after calling lemur than before. What is happening?!

This can happen if the data has large compositional shifts (for example, if one cell type disappears). The problem is that in this case, the initial linear regression step, which centers the conditions relative to each other, overcorrects and introduces a consistent shift in the latent space. You can either use align_by_grouping / align_harmony to correct for this effect or manually fix the regression coefficient to zero:

fit <- lemur(sce, design = ~ patient_id + condition, n_embedding = 15, linear_coefficient_estimator = "zero")
The conditions still separate if I plot the data using UMAP / t-SNE, even after calling align_harmony / align_neighbors. What should I do?

You can try to increase n_embedding. If this still does not help, the experimental design or the data quality may be unsuitable for inferring differential expression neighbourhoods. But as I haven’t encountered such a dataset yet, I would like to try it out myself. If you can share the data, please open an issue.

How do I make lemur faster?

Several parameters influence the compute resource usage of fitting the LEMUR model and finding neighbourhoods. Here some suggestions and hints:

  • Make sure that your data is stored in memory (not a DelayedArray) either as a sparse dgCMatrix or dense matrix.
  • A larger test_fraction means fewer cells are used to fit the model (and more cells are used for the DE test), which speeds up many steps.
  • A smaller n_embedding reduces the latent dimensions of the fit, which makes the model less flexible, but speeds up the lemur call.
  • Providing a pre-calculated set of matching cells and calling align_grouping is faster than align_harmony.
  • Setting selection_procedure = "contrast" in find_de_neighborhoods often produces better neighbourhoods, but is a lot slower than selection_procedure = "zscore".
  • Setting size_factor_method = "ratio" in find_de_neighborhoods makes the DE more powerful, but is a lot slower than size_factor_method = "normed_sum".

Session Info

sessionInfo()
R version 4.6.1 (2026-06-24)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 26.04 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8       
 [4] LC_COLLATE=en_US.UTF-8     LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C              
[10] LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

time zone: Etc/UTC
tzcode source: system (glibc)

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] SingleCellExperiment_1.35.2 SummarizedExperiment_1.43.0 Biobase_2.73.1             
 [4] GenomicRanges_1.65.1        Seqinfo_1.3.0               IRanges_2.47.2             
 [7] S4Vectors_0.51.5            BiocGenerics_0.59.10        generics_0.1.4             
[10] MatrixGenerics_1.25.0       matrixStats_1.5.0           lubridate_1.9.5            
[13] forcats_1.0.1               stringr_1.6.0               dplyr_1.2.1                
[16] purrr_1.2.2                 readr_2.2.0                 tidyr_1.3.2                
[19] tibble_3.3.1                ggplot2_4.0.3               tidyverse_2.0.0            
[22] lemur_1.11.2               

loaded via a namespace (and not attached):
 [1] gtable_0.3.6              xfun_0.60                 lattice_0.22-9           
 [4] tzdb_0.5.0                vctrs_0.7.3               tools_4.6.1              
 [7] pkgconfig_2.0.3           Matrix_1.7-5              RColorBrewer_1.1-3       
[10] S7_0.2.2                  sparseMatrixStats_1.25.0  lifecycle_1.0.5          
[13] compiler_4.6.1            farver_2.1.2              codetools_0.2-20         
[16] glmGamPoi_1.25.1          htmltools_0.5.9           sys_3.4.3                
[19] buildtools_1.0.0          yaml_2.3.12               pillar_1.11.1            
[22] MASS_7.3-65               uwot_0.2.4                DelayedArray_0.39.3      
[25] abind_1.4-8               RSpectra_0.16-2           tidyselect_1.2.1         
[28] digest_0.6.39             stringi_1.8.7             splines_4.6.1            
[31] labeling_0.4.3            maketools_1.3.2           fastmap_1.2.0            
[34] grid_4.6.1                cli_3.6.6                 SparseArray_1.13.2       
[37] magrittr_2.0.5            S4Arrays_1.13.0           utf8_1.2.6               
[40] withr_3.0.3               DelayedMatrixStats_1.35.0 scales_1.4.0             
[43] timechange_0.4.0          rmarkdown_2.31            XVector_0.53.0           
[46] otel_0.2.0                hms_1.1.4                 beachmat_2.29.0          
[49] evaluate_1.0.5            knitr_1.51                RcppAnnoy_0.0.23         
[52] irlba_2.3.7               rlang_1.3.0               isoband_0.3.0            
[55] Rcpp_1.1.2                glue_1.8.1                jsonlite_2.0.0           
[58] R6_2.6.1

  1. The set.seed step, which resets the random number generator to a defined state, is not strictly necessary, and arguably should be omitted when you apply lemur to your own data. Here, we do it to ensure that this package vignette is rendered the same for everyone. The current implementation of lemur includes some random sampling steps.↩︎

  2. There is a balancing act in the choice of test_fraction: smaller values increase the sensitivity to detect more subtle patterns in the latent space, larger values increase the statistical power of the differentially expression testing.↩︎