---
title: "Introduction"
author: Constantin Ahlmann-Eltze
date: "`r format(Sys.time(), '%d %B, %Y')`"
format:
  html:
    embed-resources: true
    page-layout: full
    toc: true
highlight-style: tango
vignette: >
  %\VignetteIndexEntry{Introduction}
  %\VignetteEngine{quarto::html}}
  %\VignetteEncoding{UTF-8}
lightbox: true
---

<!-- 'Stangle' analogon: knitr::purl("Introduction.qmd")
-->

```{r}
#| label: setup_alt_IS_THIS_REALLY_NEEDED
#| echo: false
#| eval: false
if(! exists("..options_set") || isFALSE(..options_set)){
  knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>",
    dpi = 40,
    fig_retina = 1,
    dev = "jpeg"
    # dev.args = list(quality = 20)
  )
  ..options_set <- TRUE
}
```
```{r}
#| label: initialize
#| echo: false
#| cache: false
knitr::opts_chunk$set(cache = TRUE, autodep = TRUE)
options(width = 88)
```


# Latent Embedding Multivariate Regression (LEMUR)

<!-- badges: start -->

<!-- badges: end -->

![](../man/figures/lemur-art.jpg){fig-align="center" width="80%"}

[`lemur`](https://doi.org/10.1038/s41588-024-01996-0) 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](https://bioconductor.org/packages/glmGamPoi/)
or [edgeR](https://bioconductor.org/packages/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](https://doi.org/10.1038/s41588-024-01996-0).

![Schematic of the matrix decomposition at the core of LEMUR](man/figures/equation_schematic.png){width="100%"}

## Installation

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

```{r}
#| label: install1
#| eval: false
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`:

```{r}
#| label: install2
#| eval: false
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](https://support.bioconductor.org/) or 
open an issue in the [GitHub repository](https://github.com/const-ae/lemur/) 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.

```{r} 
#| label: preparation
#| echo: false
library("lemur")
set.seed(42)
data("glioblastoma_example_data")
sce <- glioblastoma_example_data[1:50, sample.int(5000, size = 500)]
```

```{r}
#| label: quick_start
#| message: false
#| warning: false
# ... 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)](https://doi.org/10.1186/s13073-021-00894-y). 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 packages^[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.].

```{r}
#| label: load_packages
#| message: false
#| warning: false
library("tidyverse")
library("SingleCellExperiment")
library("lemur")
set.seed(42)
```

```{r}
#| label: load_data
data("glioblastoma_example_data", package = "lemur")
glioblastoma_example_data
```

As you can see above, the `lemur` package ships with a reduced-size version of the glioblastoma data, with `r nrow(glioblastoma_example_data)` genes and `r ncol(glioblastoma_example_data)` cells. A UMAP visualisation indicates that the cells' data separate them by the covariates `patient_id` and `condition`. 

```{r}
#| label: fig-raw_umap
#| fig-cap: "UMAP of `glioblastoma_example_data` prior to between-condition alignment."
#| fig-width: 6
#| fig-height: 4
#| out-width: "80%"
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()
```

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 testing^[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.].

```{r}
#| label: fit_lemur
#| message: false
fit <- lemur(glioblastoma_example_data, design = ~ patient_id + condition, 
             n_embedding = 15, test_fraction = 0.5)
fit
```

The `lemur` function returns an object of class `r class(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`:

```{r}
#| label: fitembedding
fit$embedding |> str()
```

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.

```{r}
#| label: align_lemur
#| message: false
fit <- align_harmony(fit)
```

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

```{r}
#| label: fig-lemur_umap
#| fig-cap: "UMAPs of `fit$embedding`. The points are shown separately for the two conditions, but reside in the same latent space. In contrast to @fig-raw_umap, 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."
#| fig-width: 7
#| fig-height: 3
#| out-width: "100%"
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()
```

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.

```{r}
#| label: lemur_compute_contrasts
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:

```{r}
#| label: fig-umap_de
#| fig-subcap: 
#|   - "Superimposed on the same UMAP plot as in @fig-lemur_umap"
#|   - "The histogram shows that the values are negative for most cells."
#| fig-cap: "Differential expression (log fold changes) of _GAP43_."
#| fig-width: 6
#| fig-height: 4
#| out-width: "100%"
#| layout: "[65, 35]"
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 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.
<!-- plus the sample id (in this case `patient_id`) 
 # WH 2026-03-20: I removed this subclause as it seems unneeded and confusing. -->

```{r}
#| label: fig-de_neighbourhoods
#| message: false 
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.

```{r}
#| label: topgenes
as_tibble(neighborhoods) |>
  left_join(as_tibble(rowData(fit)[,1:2]), by = c("name" = "gene_id")) |>
  relocate(symbol, .before = "name") |>
  arrange(pval) |>
  head(5)
```

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

```{r}
#| label: fig-umap_de2
#| fig-cap: "Differential expression of _SERPINE1_ superimposed on the same UMAP plot as in @fig-lemur_umap."
#| fig-width: 6
#| fig-height: 4
#| out-width: "80%"
sel_gene <- "ENSG00000106366" # SERPINE1

dplyr::filter(neighborhoods, name == sel_gene) |>
  dplyr::select(name, adj_pval, lfc, did_adj_pval, did_lfc)

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
```

<!--
# WH 2026-03-20: I commented the following out since it seems too flimsy and unconvincing.

Next, we are going to try something ambitious: in the LEMUR model, the cells in a neighbourhood are separated from the rest of the cells by a $(k-1)$-dimensional hyperplane in the $k$-dimensional latent space ($k$ being the same as `n_embedding` from above, i.e., $k=$ `r fit$n_embedding`). We can try to approximate this separation as a line in the two-dimensional UMAP plot.

To this end, we create a helper dataframe and use the `geom_density2d` function from `ggplot2`. To avoid the cutting of the boundary to the extremes of the cell coordinates, add `lims` to the plot with an appropriately large limit.

```{r}
#| label: fig-umap_de3
#| fig-cap: "Same as @fig-umap_de2, with an attempt to draw a neighbourhood boundary."
#| fig-width: 6
#| fig-height: 4
#| out-width: "80%"
neighborhood_coordinates <- neighborhoods |>
  dplyr::filter(name == sel_gene) |>
  unnest(c(neighborhood)) |>
  dplyr::rename(cell_id = neighborhood) |>
  left_join(tibble(cell_id = rownames(umap), umap), by = "cell_id") |>
  dplyr::select(name, cell_id, umap)

p + geom_density2d(data = neighborhood_coordinates, breaks = seq(0.2, 0.6, by = 0.1), 
                   contour_var = "ndensity", colour = "#808080") 
```
-->

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

```{r}
#| label: fig-volcano_plot
#| fig-cap: "Volcano plot, each point corresponds to one of the genes in our dataset, and the neighbourhood that was found for it."
#| fig-width: 4
#| fig-height: 2.5
#| out-width: "50%"
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"))
```

<!-- WH 2026-03-20: I deactivated the following chunk since it seems tangential. -->

```{r}
#| label: fig-neighbourhood_size_vs_significance
#| fig-cap: "neighbourhood size vs neighbourhood significance."
#| fig-width: 6
#| fig-height: 4
#| out-width: "50%"
#| eval: false
#| echo: false
neighborhoods |>
  drop_na() |>
  ggplot(aes(x = n_cells, y = -log10(pval))) +
    geom_point(aes(colour = adj_pval < 0.1)) 
```

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

```{r}
#| label: fig-tumour_cell_annotation1
#| fig-cap: "A simple gating strategy to putatively annotate the tumour cells in our data"
#| fig-width: 6
#| fig-height: 4
#| out-width: "60%"
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"))
```
```{r}
#| label: fig-tumour_cell_annotation2
#| fig-cap: "The tumour cells appear enriched in the larger of the blobs."
#| fig-width: 6
#| fig-height: 4
#| out-width: "60%"
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"))
```

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

```{r}
#| label: tumour_de_neighbourhood
#| paged.print: false
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)
```

Let's, for example, explore the results for CXCL8. We see in @fig-tumour_de_neighbourhood_plot 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.

```{r}
#| label: fig-tumour_de_neighbourhood_plot
#| paged.print: false
#| fig-cap: "CXCL8 expression in different cell subsets"
#| fig-width: 7
#| fig-height: 4
#| out-width: "80%"
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) 

```


# 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](https://github.com/satijalab/sctransform) instead of log-transformed data?

Yes. You can call `lemur` with any variance stabilized count matrix. Based on a [previous project](https://www.biorxiv.org/content/10.1101/2021.06.24.449781v4), 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:

```{r}
#| label: fix_linear_coef
#| message: false
#| warning: false
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

```{r}
#| label: session_info
sessionInfo()
```