set-up

library(SingleCellExperiment) # manipulate sce
library(scater) # plots
library(here) # reproducible paths
library(miloR) # DA abundance package
library(dplyr) # distinct
library(patchwork) #plots
library(dplyr)

Milo vignette followed for this analyisis: miloR

source(here("src/colours.R"))
project <- "fire-mice"
fig_path <- here("outs", project, "plots","DA_miloR", "/")
sce <- readRDS(here("processed", project, "sce_anno_02.RDS"))

Differential abundance testing for KOvsWT

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.

genotypekowt <- sce$genotype %in% c("KO","WT")
milo_ko <- Milo(sce[,genotypekowt])
colData(milo_ko)[["genotype"]] <- droplevels(colData(milo_ko)[["genotype"]])

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. In this case, we specify that we want to build the graph from the MNN corrected PCA dimensions.

For graph building you need to define a few parameters:

d: the number of reduced dimensions to use for KNN refinement. We recommend using the same \(d\) used for KNN graph building.In our corrected dataset we have 23 dimensions.
k: this affects the power of DA testing, since we need to have enough cells from each sample represented in a neighbourhood to estimate the variance between replicates. On the other side, increasing \(k\) too much might lead to over-smoothing. We suggest to start by using the same value for \(k\) used for KNN graph building for clustering and UMAP visualization. We will later use some heuristics to evaluate whether the value of \(k\) should be increased.

milo_ko <- buildGraph(milo_ko, k = 30, d = 23, reduced.dim = "corrected")

## Constructing kNN graph with k:30

Defining representative neighbourhoods on the KNN graph

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.

As well as \(d\) and \(k\), for sampling we need to define a few additional parameters:

prop: the proportion of cells to randomly sample to start with. We suggest using prop=0.1 for datasets of less than 30k cells.
refined: indicates 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.

milo_ko <- makeNhoods(milo_ko, prop = 0.1, k = 30, d=23, refined = TRUE, reduced_dims = "corrected")

## Checking valid object

Once we have defined neighbourhoods, we plot the distribution of neighbourhood sizes (i.e. how many cells form each neighbourhood) to evaluate whether the value of \(k\) used for graph building was appropriate. We can check this out using the plotNhoodSizeHist function.

As a rule of thumb we want to have an average neighbourhood size over 5 x N_samples. (5x12=60)

plotNhoodSizeHist(milo_ko)

Counting cells in neighbourhoods

Milo leverages the variation in cell numbers between replicates for the same experimental condition to test for differential abundance. Therefore 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.

milo_ko <- countCells(milo_ko, meta.data = as.data.frame(colData(milo_ko)), sample="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(milo_ko))

## 6 x 8 sparse Matrix of class "dgCMatrix"
##   EKDL2651 EKDL2653 EKDL2654 EKDL2656 EKDL2657 EKDL2659 EKDL2660 EKDL2662
## 1       16        3       13        1        2        5        3        2
## 2        5        4        6        3        4        7        6        3
## 3        8        9        9        1        2        7        3        5
## 4       38       16       31        9        9       18        7        4
## 5       14        6       43       27        .        7       31        8
## 6       10        7       11       17        4        4       21        5

Defining experimental design

For this section we refer back to the contrast vignette

The design matrix should match each sample to the experimental condition of interest for DA testing. We also include the sequencing.batch column in the design matrix. This represents a known technical covariate that we want to account for in DA testing.

design <- data.frame(colData(milo_ko))[,c("Sample", "batch", "genotype")]
design$batch <- as.factor(design$batch)
design <- distinct(design)
rownames(design) <- design$Sample
contrastKO <- c("genotypeKO - genotypeWT") # the syntax is

Computing neighbourhood connectivity

Milo uses an adaptation of the Spatial FDR correction introduced by cydar, where we correct p-values accounting for the amount of 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. This is done by the calcNhoodDistance function N.B. this step is the most time consuming of the analysis workflow and might take a “couple” (this takes ages! ) of minutes for large (or not that large) datasets.

milo_ko <- calcNhoodDistance(milo_ko, d=23, reduced.dim = "corrected")

saveRDS(milo_ko, here("processed", project, "milo_ko.rds"))

Testing KO

# we need to use the ~ 0 + Variable expression here so that we have all of the levels of our variable as separate columns in our model matrix
da_results_ko <- testNhoods(milo_ko, design = ~ 0  + genotype + batch, design.df = design, model.contrasts = contrastKO, reduced.dim = "corrected")

## Using TMM normalisation

## Performing spatial FDR correction withk-distance weighting

table(da_results_ko$SpatialFDR < 0.1)

## 
## FALSE  TRUE 
##  1945   500

This calculates a Fold-change and corrected P-value for each neighbourhood, which indicates whether there is significant differential abundance between developmental stages. The main statistics we consider here are:

logFC: indicates the log-Fold change in cell numbers between samples from E7.5 and samples from E7.0
PValue: reports P-values before FDR correction
SpatialFDR: reports P-values corrected for multiple testing accounting for overlap between neighbourhoods

da_results_ko %>%
  arrange(SpatialFDR) %>%
  head()

##          logFC    logCPM        F       PValue          FDR Nhood   SpatialFDR
## 32   -8.692787 10.305361 47.20832 6.906301e-12 1.039088e-08    32 8.353722e-09
## 480  -8.524413 10.109033 45.42693 1.704576e-11 1.041922e-08   480 8.353722e-09
## 707  -8.540875 10.198406 45.99934 1.274954e-11 1.039088e-08   707 8.353722e-09
## 1807 -8.573780 10.208376 46.32634 1.080098e-11 1.039088e-08  1807 8.353722e-09
## 499  -8.450609 10.073669 44.43143 2.825260e-11 1.381552e-08   499 1.108741e-08
## 1304 -8.208373  9.922538 43.94743 3.612398e-11 1.472052e-08  1304 1.201859e-08

Inspecting DA testing results KOvsWT

We can start inspecting the results of our DA analysis from a couple of standard diagnostic plots. We first inspect the distribution of uncorrected P values, to verify that the test was balanced.

ggplot(da_results_ko, aes(PValue)) + geom_histogram(bins=50) + ggtitle(contrastKO)

Then we visualize the test results with a volcano plot (remember that each point here represents a neighbourhood, not a cell).

ggplot(da_results_ko, aes(logFC, -log10(SpatialFDR))) + 
  geom_point() +
  geom_hline(yintercept = 1) ## Mark significance threshold (10% FDR)

To visualize DA results relating them to the embedding of single cells, we can build an abstracted graph of neighbourhoods that we can superimpose on the single-cell embedding. Here each node represents a neighbourhood, while edges indicate how many cells two neighbourhoods have in common. Here the layout of nodes is determined by the position of the index cell in the UMAP embedding of all single-cells. The neighbourhoods displaying significant DA are colored by their log-Fold Change.

milo_ko <- buildNhoodGraph(milo_ko)
## Plot single-cell UMAP
tsne_pl <- plotReducedDim(milo_ko, dimred = "TSNE", colour_by="genotype", text_by = "celltype", 
                          text_size = 3, point_size=0.5) +
  guides(fill="none")  + scale_color_manual(values = c(col_wt_het_ko[1],col_wt_het_ko[3])) + labs(color="genotype")

## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.

## Plot neighbourhood graph
nh_graph_pl <- plotNhoodGraphDA(milo_ko, da_results_ko, layout="TSNE",alpha=0.1) +  scale_fill_gradient2(high = scales::muted("red"), mid = "white", low = scales::muted("blue")) + labs(fill = "logFC")

## Scale for 'fill' is already present. Adding another scale for 'fill', which
## will replace the existing scale.

tsne_pl + nh_graph_pl +
  plot_layout(guides="collect")

We might also be interested in visualizing whether DA is particularly evident in certain cell types. To do this, we assign a cell type label to each neighbourhood by finding the most abundant cell type within cells in each neighbourhood. We can label neighbourhoods in the results data.frame using the function annotateNhoods. This also saves the fraction of cells harbouring the label.

da_results_ko <- annotateNhoods(milo_ko, da_results_ko, coldata_col = "clusters_named")
head(da_results_ko)

##          logFC   logCPM            F     PValue       FDR Nhood SpatialFDR
## 1 -1.416412279 8.426049 2.9132643409 0.08789621 0.3083303     1  0.3010076
## 2 -0.369591374 8.305650 0.2227145766 0.63699332 0.8794310     2  0.8736272
## 3  0.009624541 8.432442 0.0001260405 0.99104282 0.9967243     3  0.9968246
## 4 -0.775972214 9.679941 1.1639268335 0.28068811 0.6178972     4  0.6069391
## 5 -0.526154209 9.812559 0.4848690863 0.48624712 0.7973831     5  0.7887570
## 6 -0.472514196 9.130728 0.4447783609 0.50484609 0.8115376     6  0.8032122
##   clusters_named clusters_named_fraction
## 1       mOligo_1               1.0000000
## 2       mOligo_2               0.6842105
## 3       mOligo_1               1.0000000
## 4       mOligo_1               1.0000000
## 5        Astro_1               1.0000000
## 6        Astro_1               1.0000000

While neighbourhoods tend to be homogeneous, we can define a threshold for clusters_named_fraction to exclude neighbourhoods that are a mix of cell types.

ggplot(da_results_ko, aes(clusters_named_fraction)) + geom_histogram(bins=50)

da_results_ko$clusters_named <- ifelse(da_results_ko$clusters_named_fraction < 0.5, "Mixed", da_results_ko$clusters_named)

Now we can visualize the distribution of DA Fold Changes in different cell types

da_results_ko$clusters_named <- factor(da_results_ko$clusters_named, rev(c("BAMs", "Microglia", "Immune", "Astro_1", "Astro_2", "Astro_Oligo", "OPC", "pOPC", "iOligo", "mOligo_1", "mOligo_2", "mOligo_3", "mOligo_4", "Endothelial", "Mural_cells", "ChP_epithelia", "iNeurons_&_NRPs", "mNeuron_ex", "mNeuron_in", "Mixed")))
plotDAbeeswarm(da_results_ko, group.by = "clusters_named")  +  scale_colour_gradient2(high = scales::muted("red"), mid = "white", low = scales::muted("blue")) + xlab("")

## Scale for 'colour' is already present. Adding another scale for 'colour',
## which will replace the existing scale.

Astro1

da_results_ko <- annotateNhoods(milo_ko, da_results_ko, coldata_col = "Sample")

## Converting Sample to factor...

head(da_results_ko)

##          logFC   logCPM            F     PValue       FDR Nhood SpatialFDR
## 1 -1.416412279 8.426049 2.9132643409 0.08789621 0.3083303     1  0.3010076
## 2 -0.369591374 8.305650 0.2227145766 0.63699332 0.8794310     2  0.8736272
## 3  0.009624541 8.432442 0.0001260405 0.99104282 0.9967243     3  0.9968246
## 4 -0.775972214 9.679941 1.1639268335 0.28068811 0.6178972     4  0.6069391
## 5 -0.526154209 9.812559 0.4848690863 0.48624712 0.7973831     5  0.7887570
## 6 -0.472514196 9.130728 0.4447783609 0.50484609 0.8115376     6  0.8032122
##   clusters_named clusters_named_fraction   Sample Sample_fraction
## 1       mOligo_1               1.0000000 EKDL2651       0.3555556
## 2       mOligo_2               0.6842105 EKDL2659       0.1842105
## 3       mOligo_1               1.0000000 EKDL2653       0.2045455
## 4       mOligo_1               1.0000000 EKDL2651       0.2878788
## 5        Astro_1               1.0000000 EKDL2654       0.3161765
## 6        Astro_1               1.0000000 EKDL2660       0.2658228

da_results_ko %>% filter(clusters_named == "Astro_1") %>% filter(SpatialFDR<0.1) %>% filter(logFC<0)

##          logFC   logCPM        F      PValue        FDR Nhood SpatialFDR
## 28   -1.674167 9.547669 5.377601 0.020424247 0.10017371    28 0.09754041
## 89   -2.185764 8.775484 8.235711 0.004119150 0.02462426    89 0.02406802
## 317  -1.839686 9.277276 6.117850 0.013404828 0.07033220   317 0.06856572
## 471  -1.968074 9.023544 6.324977 0.011926110 0.06325236   471 0.06164098
## 952  -1.824557 8.983894 5.934862 0.014867865 0.07677381   952 0.07486769
## 1476 -2.279488 8.757116 8.444411 0.003672625 0.02233723  1476 0.02182178
## 1605 -2.017823 8.877335 6.664466 0.009854637 0.05390288  1605 0.05258510
##      clusters_named clusters_named_fraction   Sample Sample_fraction
## 28          Astro_1               1.0000000 EKDL2660       0.3888889
## 89          Astro_1               1.0000000 EKDL2651       0.3500000
## 317         Astro_1               1.0000000 EKDL2654       0.3483146
## 471         Astro_1               1.0000000 EKDL2660       0.3380282
## 952         Astro_1               1.0000000 EKDL2654       0.2500000
## 1476        Astro_1               0.9298246 EKDL2660       0.2982456
## 1605        Astro_1               1.0000000 EKDL2660       0.3064516

The nhoods with DA are made of a mix of samples.

Testing HET vs WT

genotypehetwt <- sce$genotype %in% c("HET","WT")
milo_het <- Milo(sce[,genotypehetwt])
colData(milo_het)[["genotype"]] <- droplevels(colData(milo_het)[["genotype"]])
milo_het <- buildGraph(milo_het, k = 30, d = 23, reduced.dim = "corrected")

## Constructing kNN graph with k:30

milo_het <- makeNhoods(milo_het, prop = 0.1, k = 30, d=23, refined = TRUE, reduced_dims = "corrected")

## Checking valid object

plotNhoodSizeHist(milo_het)

milo_het <- countCells(milo_het, meta.data = as.data.frame(colData(milo_het)), sample="Sample")

## Checking meta.data validity

## Counting cells in neighbourhoods

head(nhoodCounts(milo_het))

## 6 x 8 sparse Matrix of class "dgCMatrix"
##   EKDL2651 EKDL2652 EKDL2654 EKDL2655 EKDL2657 EKDL2658 EKDL2660 EKDL2661
## 1       15        4        8        7        2        6        6        5
## 2        7        2       10        8        1        5        2        4
## 3       11        .        8       15        1        5       11        4
## 4       12        1       11       18        4        3        8        6
## 5       44        4        4       22        3        6       18        6
## 6       26        9        4       13        4        2        1        2

design_het <- data.frame(colData(milo_het))[,c("Sample", "batch", "genotype")]
design_het$batch <- as.factor(design_het$batch)
design_het <- distinct(design_het)
rownames(design_het) <- design_het$Sample

contrastHET <- c("genotypeHET - genotypeWT")

milo_het <- calcNhoodDistance(milo_het, d=23, reduced.dim = "corrected")

saveRDS(milo_het, here("processed", project, "milo_het.rds"))

da_results_het <- testNhoods(milo_het, design = ~ 0 + genotype + batch, design.df = design_het, model.contrasts = "genotypeHET - genotypeWT", reduced.dim = "corrected")

## Using TMM normalisation

## Performing spatial FDR correction withk-distance weighting

#source(here("src/subset_milo_het.R"))
# we need to use the ~ 0 + Variable expression here so that we have all of the levels of our variable as separate columns in our model matrix
#da_results_het <- testNhoods(milo_het, design = ~ 0 + genotype + batch, design.df = design_het, model.contrasts = "genotypeHET - genotypeWT", reduced.dim = "corrected")
table(da_results_het$SpatialFDR < 0.1)

## 
## FALSE 
##  2422

ggplot(da_results_het, aes(PValue)) + geom_histogram(bins=50) + ggtitle("genotypeWT - genotypeHET") + scale_y_continuous(limits = c(0,700)) + geom_abline(intercept=dim((da_results_het))[1]/50)

Visualise the results (HETvsWT)