# Differential abundance testing with Milo

Emma Dann and Mike Morgan

#### 29 October 2024

#### Package

miloR 2.2.0

    library(miloR)
    library(SingleCellExperiment)
    library(scater)
    library(scran)
    library(dplyr)
    library(patchwork)

# 1 Introduction

Milo is a tool for analysis of complex single cell datasets generated from replicated multi-condition experiments, which detects changes in composition between conditions. While differential abundance (DA) is commonly quantified in discrete cell clusters, Milo uses partially overlapping neighbourhoods of cells on a KNN graph. Starting from a graph that faithfully recapitulates the biology of the cell population, Milo analysis consists of 3 steps:

1. Sampling of representative neighbourhoods
2. Testing for differential abundance of conditions in all neighbourhoods
3. Accounting for multiple hypothesis testing using a weighted FDR procedure that accounts for the overlap of neighbourhoods

In this vignette we will elaborate on how these steps are implemented in the `miloR` package.

# 2 Load data

For this demo we will use a synthetic dataset simulating a developmental trajectory, generated using dyntoy.

    data("sim_trajectory", package = "miloR")

    ## Extract SingleCellExperiment object
    traj_sce <- sim_trajectory[['SCE']]

    ## Extract sample metadata to use for testing
    traj_meta <- sim_trajectory[["meta"]]

    ## Add metadata to colData slot
    colData(traj_sce) <- DataFrame(traj_meta)
    colnames(traj_sce) <- colData(traj_sce)$cell_id

    redim <- reducedDim(traj_sce, "PCA")
    dimnames(redim) <- list(colnames(traj_sce), paste0("PC", c(1:50)))
    reducedDim(traj_sce, "PCA") <- redim

# 3 Pre-processing

For DA analysis we need to construct an undirected KNN graph of single-cells. Standard single-cell analysis pipelines usually do this from distances in PCA. We normalize and calculate principal components using `scater`. I also run UMAP for visualization purposes.

    logcounts(traj_sce) <- log(counts(traj_sce) + 1)
    traj_sce <- runPCA(traj_sce, ncomponents=30)
    traj_sce <- runUMAP(traj_sce)

    plotUMAP(traj_sce)

# 4 Create a Milo object

For differential abundance analysis on graph neighbourhoods we first construct a `Milo` object. This extends the `SingleCellExperiment` class to store information about neighbourhoods on the KNN graph.

## 4.1 From SingleCellExperiment object

The `Milo` constructor takes as input a `SingleCellExperiment` object.

    traj_milo <- Milo(traj_sce)
    reducedDim(traj_milo, "UMAP") <- reducedDim(traj_sce, "UMAP")

    traj_milo

    ## class: Milo
    ## dim: 500 500
    ## metadata(0):
    ## assays(2): counts logcounts
    ## rownames(500): G1 G2 ... G499 G500
    ## rowData names(0):
    ## colnames(500): C1 C2 ... C499 C500
    ## colData names(5): cell_id group_id Condition Replicate Sample
    ## reducedDimNames(2): PCA UMAP
    ## mainExpName: NULL
    ## altExpNames(0):
    ## nhoods dimensions(2): 1 1
    ## nhoodCounts dimensions(2): 1 1
    ## nhoodDistances dimension(1): 0
    ## graph names(0):
    ## nhoodIndex names(1): 0
    ## nhoodExpression dimension(2): 1 1
    ## nhoodReducedDim names(0):
    ## nhoodGraph names(0):
    ## nhoodAdjacency dimension(2): 1 1

## 4.2 From AnnData object (.h5ad)

We can use the `zellkonverter` package to make a `SingleCellExperiment` object from an `AnnData` object stored as `h5ad` file.

    library(zellkonverter)

    # Obtaining an example H5AD file.
    example_h5ad <- system.file("extdata", "krumsiek11.h5ad",
                                package = "zellkonverter")

    example_h5ad_sce <- readH5AD(example_h5ad)
    example_h5ad_milo <- Milo(example_h5ad_sce)

## 4.3 From Seurat object

The `Seurat` package includes a converter to `SingleCellExperiment`.

    library(Seurat)
    data("pbmc_small")
    pbmc_small_sce <- as.SingleCellExperiment(pbmc_small)
    pbmc_small_milo <- Milo(pbmc_small_sce)

# 5 Construct KNN graph

We need to add the KNN graph to the Milo object. This is stored in the `graph` slot, in `igraph` format. The `miloR` package includes functionality to build and store the graph from the PCA dimensions stored in the `reducedDim` slot.

    traj_milo <- buildGraph(traj_milo, k = 10, d = 30)

    ## Constructing kNN graph with k:10

In progress: we are perfecting the functionality to add a precomputed KNN graph (for example constructed with Seurat or scanpy) to the `graph` slot using the adjacency matrix.

# 6 1. Defining representative neighbourhoods

We define the neighbourhood of a cell, the index, as the group of cells connected by an edge in the KNN graph to the index cell. For efficiency, we don’t test for DA in the neighbourhood of every cell, but we sample as indices a subset of representative cells, using a KNN sampling algorithm used by Gut et al. 2015.

For sampling you need to define a few parameters:

- `prop`: the proportion of cells to randomly sample to start with (usually 0.1 - 0.2 is sufficient)
- `k`: the k to use for KNN refinement (we recommend using the same k used for KNN graph building)
- `d`: the number of reduced dimensions to use for KNN refinement (we recommend using the same d used for KNN graph building)
- `refined` indicated whether you want to use the sampling refinement algorithm, or just pick cells at random. The default and recommended way to go is to use refinement. The only situation in which you might consider using random instead, is if you have batch corrected your data with a graph based correction algorithm, such as BBKNN, but the results of DA testing will be suboptimal.

    traj_milo <- makeNhoods(traj_milo, prop = 0.1, k = 10, d=30, refined = TRUE)

    ## Checking valid object

    ## Running refined sampling with reduced_dim

Once we have defined neighbourhoods, it’s good to take a look at how big the neighbourhoods are (i.e. how many cells form each neighbourhood). This affects the power of DA testing. We can check this out using the `plotNhoodSizeHist` function. Empirically, we found it’s best to have a distribution peaking between 50 and 100. Otherwise you might consider rerunning `makeNhoods` increasing `k` and/or `prop` (here the distribution looks ludicrous because it’s a small dataset).

    plotNhoodSizeHist(traj_milo)

# 7 Counting cells in neighbourhoods

Now we have to count how many cells from each sample are in each neighbourhood. We need to use the cell metadata and specify which column contains the sample information.

    traj_milo <- countCells(traj_milo, meta.data = data.frame(colData(traj_milo)), samples="Sample")

    ## Checking meta.data validity

    ## Counting cells in neighbourhoods

This adds to the `Milo` object a `n \times m` matrix, where n is the number of neighbourhoods and \(m\) is the number of experimental samples. Values indicate the number of cells from each sample counted in a neighbourhood. This count matrix will be used for DA testing.

    head(nhoodCounts(traj_milo))

    ## 6 x 6 sparse Matrix of class "dgCMatrix"
    ##   B_R1 A_R1 A_R2 B_R2 B_R3 A_R3
    ## 1    6    6    4    5    9    8
    ## 2    7    5    2    5    5    3
    ## 3    9    3    3   11   11    8
    ## 4    5    2    5    8    8    8
    ## 5    3    5    4    4    4    3
    ## 6    5    6    4    8    6    7

# 8 Differential abundance testing

Now we are all set to test for differential abundance in neighbourhoods. We implement this hypothesis testing in a generalized linear model (GLM) framework, specifically using the Negative Binomial GLM implementation in `edgeR`.

We first need to think about our experimental design. The design matrix should match samples to a condition of interest. In this case the `Condition` is the covariate we are going to test for.

    traj_design <- data.frame(colData(traj_milo))[,c("Sample", "Condition")]
    traj_design <- distinct(traj_design)
    rownames(traj_design) <- traj_design$Sample
    ## Reorder rownames to match columns of nhoodCounts(milo)
    traj_design <- traj_design[colnames(nhoodCounts(traj_milo)), , drop=FALSE]

    traj_design

    ##      Sample Condition
    ## B_R1   B_R1         B
    ## A_R1   A_R1         A
    ## A_R2   A_R2         A
    ## B_R2   B_R2         B
    ## B_R3   B_R3         B
    ## A_R3   A_R3         A

Milo uses an adaptation of the Spatial FDR correction introduced by cydar, which accounts for the overlap between neighbourhoods. Specifically, each hypothesis test P-value is weighted by the reciprocal of the kth nearest neighbour distance. To use this statistic we first need to store the distances between nearest neighbors in the Milo object.

    traj_milo <- calcNhoodDistance(traj_milo, d=30)

    ## 'as(<dgTMatrix>, "dgCMatrix")' is deprecated.
    ## Use 'as(., "CsparseMatrix")' instead.
    ## See help("Deprecated") and help("Matrix-deprecated").

Now we can do the test, explicitly defining our experimental design.

    rownames(traj_design) <- traj_design$Sample
    da_results <- testNhoods(traj_milo, design = ~ Condition, design.df = traj_design)

    ## Using TMM normalisation

    ## Running with model contrasts

    ## Performing spatial FDR correction with k-distance weighting

This calculates a Fold-change and corrected P-value for each neighbourhood, which indicates whether there is significant differential abundance between conditions.

    da_results %>%
      arrange(- SpatialFDR) %>%
      head()

    ##          logFC   logCPM          F    PValue       FDR Nhood SpatialFDR
    ## 2  -0.06738152 15.15548 0.01692111 0.9043499 0.9043499     2  0.9043499
    ## 16  0.28867335 14.67741 0.11923759 0.7304479 0.7547962    16  0.7548513
    ## 15 -0.24468397 15.42466 0.34524386 0.6180073 0.6606285    15  0.6570579
    ## 3   0.30679935 15.70291 0.58601131 0.4962218 0.5493884     3  0.5458527
    ## 4  -0.34236939 15.45504 0.85455289 0.4772886 0.5479981     4  0.5439883
    ## 21  0.43799876 15.77691 1.00478270 0.3213754 0.3831783    21  0.3801322

# 9 Visualize neighbourhoods displaying DA

To visualize DA results, we build an abstracted graph of neighbourhoods that we can superimpose on the single-cell embedding.

    traj_milo <- buildNhoodGraph(traj_milo)

    plotUMAP(traj_milo) + plotNhoodGraphDA(traj_milo, da_results, alpha=0.05) +
      plot_layout(guides="collect")

    ## Adding nhood effect sizes to neighbourhood graph attributes

Session Info

    sessionInfo()

    ## R version 4.4.1 (2024-06-14 ucrt)
    ## Platform: x86_64-w64-mingw32/x64
    ## Running under: Windows Server 2022 x64 (build 20348)
    ##
    ## Matrix products: default
    ##
    ##
    ## locale:
    ## [1] LC_COLLATE=C
    ## [2] LC_CTYPE=English_United States.utf8
    ## [3] LC_MONETARY=English_United States.utf8
    ## [4] LC_NUMERIC=C
    ## [5] LC_TIME=English_United States.utf8
    ##
    ## time zone: America/New_York
    ## tzcode source: internal
    ##
    ## attached base packages:
    ## [1] stats4    stats     graphics  grDevices utils     datasets  methods
    ## [8] base
    ##
    ## other attached packages:
    ##  [1] MouseThymusAgeing_1.13.0    patchwork_1.3.0
    ##  [3] dplyr_1.1.4                 scran_1.34.0
    ##  [5] scater_1.34.0               ggplot2_3.5.1
    ##  [7] scuttle_1.16.0              SingleCellExperiment_1.28.0
    ##  [9] SummarizedExperiment_1.36.0 Biobase_2.66.0
    ## [11] GenomicRanges_1.58.0        GenomeInfoDb_1.42.0
    ## [13] IRanges_2.40.0              S4Vectors_0.44.0
    ## [15] BiocGenerics_0.52.0         MatrixGenerics_1.18.0
    ## [17] matrixStats_1.4.1           miloR_2.2.0
    ## [19] edgeR_4.4.0                 limma_3.62.0
    ## [21] BiocStyle_2.34.0
    ##
    ## loaded via a namespace (and not attached):
    ##   [1] RColorBrewer_1.1-3      jsonlite_1.8.9          magrittr_2.0.3
    ##   [4] magick_2.8.5            ggbeeswarm_0.7.2        farver_2.1.2
    ##   [7] rmarkdown_2.28          zlibbioc_1.52.0         vctrs_0.6.5
    ##  [10] memoise_2.0.1           tinytex_0.53            htmltools_0.5.8.1
    ##  [13] S4Arrays_1.6.0          AnnotationHub_3.14.0    curl_5.2.3
    ##  [16] BiocNeighbors_2.0.0     SparseArray_1.6.0       sass_0.4.9
    ##  [19] pracma_2.4.4            bslib_0.8.0             cachem_1.1.0
    ##  [22] igraph_2.1.1            mime_0.12               lifecycle_1.0.4
    ##  [25] pkgconfig_2.0.3         rsvd_1.0.5              Matrix_1.7-1
    ##  [28] R6_2.5.1                fastmap_1.2.0           GenomeInfoDbData_1.2.13
    ##  [31] digest_0.6.37           numDeriv_2016.8-1.1     colorspace_2.1-1
    ##  [34] AnnotationDbi_1.68.0    dqrng_0.4.1             irlba_2.3.5.1
    ##  [37] ExperimentHub_2.14.0    RSQLite_2.3.7           beachmat_2.22.0
    ##  [40] labeling_0.4.3          filelock_1.0.3          fansi_1.0.6
    ##  [43] httr_1.4.7              polyclip_1.10-7         abind_1.4-8
    ##  [46] compiler_4.4.1          bit64_4.5.2             withr_3.0.2
    ##  [49] BiocParallel_1.40.0     viridis_0.6.5           DBI_1.2.3
    ##  [52] highr_0.11              ggforce_0.4.2           MASS_7.3-61
    ##  [55] rappdirs_0.3.3          DelayedArray_0.32.0     bluster_1.16.0
    ##  [58] gtools_3.9.5            tools_4.4.1             vipor_0.4.7
    ##  [61] beeswarm_0.4.0          glue_1.8.0              grid_4.4.1
    ##  [64] cluster_2.1.6           generics_0.1.3          gtable_0.3.6
    ##  [67] tidyr_1.3.1             BiocSingular_1.22.0     tidygraph_1.3.1
    ##  [70] ScaledMatrix_1.14.0     metapod_1.14.0          utf8_1.2.4
    ##  [73] XVector_0.46.0          ggrepel_0.9.6           BiocVersion_3.20.0
    ##  [76] pillar_1.9.0            stringr_1.5.1           splines_4.4.1
    ##  [79] tweenr_2.0.3            BiocFileCache_2.14.0    lattice_0.22-6
    ##  [82] FNN_1.1.4.1             bit_4.5.0               tidyselect_1.2.1
    ##  [85] locfit_1.5-9.10         Biostrings_2.74.0       knitr_1.48
    ##  [88] gridExtra_2.3           bookdown_0.41           xfun_0.48
    ##  [91] graphlayouts_1.2.0      statmod_1.5.0           stringi_1.8.4
    ##  [94] UCSC.utils_1.2.0        yaml_2.3.10             evaluate_1.0.1
    ##  [97] codetools_0.2-20        ggraph_2.2.1            tibble_3.2.1
    ## [100] BiocManager_1.30.25     cli_3.6.3               uwot_0.2.2
    ## [103] munsell_0.5.1           jquerylib_0.1.4         Rcpp_1.0.13
    ## [106] dbplyr_2.5.0            png_0.1-8               parallel_4.4.1
    ## [109] blob_1.2.4              viridisLite_0.4.2       scales_1.3.0
    ## [112] purrr_1.0.2             crayon_1.5.3            rlang_1.1.4
    ## [115] cowplot_1.1.3           KEGGREST_1.46.0