Introduction to single-cell RNA-seq analysis
11 Differential Abundance Analysis
In the previous section we discussed how we can perform differential expression analysis using pseudo-bulk samples, when multiple biological replicates are available in a multi-condition experiment. Differential expression analysis addresses the question of whether a gene has a different average expression between the two groups being compared. However, another useful question to address is whether the cell composition also differs between conditions. We could imagine that B-cell precursor populations expand in leukaemia samples, leading to reduced relative abundance of T and NK cells.
This type of analysis is referred to as differential abundance analysis.
We start our analysis, as always, by loading our packages:
# load packages
library(BiocParallel)
library(scran)
library(scater)
library(miloR)
library(tidyverse)
library(patchwork)We will continue working with our PBMMC and ETV6-RUNX1 samples, as we did in the differential expression section. Note that we’ve modified the data used in this section to make sure we have some significant DA results to show (in reality, the original data did not show any differences).
# Load data
# Note: this data has been modified for the purpose of this demonstration
# To make sure we have some significant DA results to show.
sce <- readRDS("RObjects/MiloDemo.sce.rds")One method to perform DA analysis is to simply count how many cells occur in each label + condition combination and then test for differences in the counts between the two conditions.
table(sce$label, sce$SampleName)##
## ETV6RUNX1-1 ETV6RUNX1-2 ETV6RUNX1-3 ETV6RUNX1-4 PBMMC-1 PBMMC-2 PBMMC-3
## B (c3) 548 544 563 550 66 61 48
## B (c6) 316 495 61 264 133 100 120
## B (c5) 115 327 124 910 86 76 85
## Monocytes (c14) 5 6 3 9 244 286 258
## T (c2) 2 28 336 40 1180 1172 1201
## Erythrocytes (c8) 4 17 259 898 404 426 378
## CD20+ B (c10) 7 113 464 141 420 401 392
## NK T (c11) 3 13 51 27 390 407 440
## T (c19) 11 61 24 18 37 29 37
## Unknown (c9) 16 99 58 16 3 1 0
## B (c1) 6 4380 70 16 356 324 331
## T (c13) 7 41 49 28 3 3 0
## B (c4) 7 43 81 1438 97 100 93
## B (c7) 3 464 75 12 1 2 8
## B (c12) 6 76 1580 84 27 27 29
## B (c16) 3 105 7 850 7 10 6
## Monocytes (c21) 1 3 15 36 70 93 91
## T (c18) 0 9 10 2 1 1 0
## B (c23) 0 23 8 5 30 24 28
## Erythrocytes (c15) 0 5 80 909 8 7 13
## B (c22) 0 2 0 1 26 30 21
## Erythrocytes (c17) 0 0 3 6 319 338 319
## T (c20) 0 0 0 0 132 107 109As these are count data, statistical methods used in flow-cytometry have been adapted to test for differences in cell abundance from such a matrix of counts. This is the approach used in the Differential Abundance Analysis chapter in the OSCA book. However, this approach does not take into account the fact that cells are not independent observations, but rather are related to each other in a complex way.
This approach also relies on a fixed set of clusters determined by us beforehand, which can be limiting in cases where the changes in abundance are more continuous (e.g. along a developmental trajectory, or where cell states are not completely discrete). To address this limitation, Dann et al. (2022) developed a method where cell abundance differences are tested along neighbourhoods of cells on a KNN graph. This means that the results aren’t dependent on our clustering results, and are instead treated in a more “continuous” way.
To address this issue, we can use a method called Milo.
This method has been implemented by the authors in the R/Bioconductor package MiloR, which we cover in this section. Starting from a graph that faithfully recapitulates the biology of the cell population, the Milo analysis consists of 3 steps:
- Building a k-nearest neighbours (KNN) graph
- Sampling representative neighbourhoods in the graph (for computational efficiency)
- Testing for differential abundance of conditions in all neighbourhoods
- Accounting for multiple hypothesis testing using a weighted FDR procedure that accounts for the overlap of neighbourhoods
As the package was written by the authors of the scater/scran suite of packages for Bioconductor the structures unpinning it are based on SingleCellExperiment. The previous add-on section demonstrated how we can easily convert our Seurat object to a SingleCellExperiment object. The next step in the analysis is to turn our single cell object into a Milo object. This is very similar to the SingleCellExperiment object, but also includes information about the neighbourhood graphs that are built during the analysis.
# check the contents of the object
sce## class: SingleCellExperiment
## dim: 28343 30167
## metadata(0):
## assays(2): counts logcounts
## rownames(28343): ENSG00000238009 ENSG00000241860 ... FAM41AY1 FAM224B
## rowData names(0):
## colnames(30167): ETV6RUNX1-1_AAACCTGAGACTTTCG-1 ETV6RUNX1-1_AAACCTGGTCTTCAAG-1 ... PBMMC-3_TTTGTCATCAGTTGAC-1 PBMMC-3_TTTGTCATCTCGTTTA-1
## colData names(12): orig.ident nCount_RNA ... ident label
## reducedDimNames(3): PCA HARMONY UMAP
## mainExpName: SCT
## altExpNames(1): RNA# plot UMAP done on the batch-corrected data
plotReducedDim(sce, dimred = "UMAP",
colour_by = "label",
text_by = "label")
# convert to a Milo object
milo <- Milo(sce)
milo## class: Milo
## dim: 28343 30167
## metadata(0):
## assays(2): counts logcounts
## rownames(28343): ENSG00000238009 ENSG00000241860 ... FAM41AY1 FAM224B
## rowData names(0):
## colnames(30167): ETV6RUNX1-1_AAACCTGAGACTTTCG-1 ETV6RUNX1-1_AAACCTGGTCTTCAAG-1 ... PBMMC-3_TTTGTCATCAGTTGAC-1 PBMMC-3_TTTGTCATCTCGTTTA-1
## colData names(12): orig.ident nCount_RNA ... ident label
## reducedDimNames(3): PCA HARMONY UMAP
## mainExpName: SCT
## 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 1Notice how there are now several slots with prefix nhoods, which we will explore as we go through the analysis.
11.1 Construct KNN graph
# Build KNN graph ----The first step in our analysis is to build a KNN graph from our cells. This is very similar to what we did earlier in the clustering session, except now the graph will be stored inside the object:
# add KNN graph to Milo object
milo <- buildGraph(milo,
k = 48,
d = 20,
reduced.dim = "HARMONY",
BPPARAM = MulticoreParam(7))## Constructing kNN graph with k:48A couple of things to note about this:
kis the number of nearest neighbours to build the graph. This can be adjusted depending on how much of a fine resolution you want. If you use too high a number, you will end up loosing resolution as a higher diversity of cells will be connected in the graph. On the other hand, if you use too low a number, you may increase noise in the data and loose statistical power, as only very few cells will be connected to each other in a neighbourhood. The author’s recommendation is to use a value of k as you did for clustering and UMAP visualisation - this is what we’ve done in this case.dis the number of dimensions from our dimensionality reduction matrix to use. In this case we’re using the MNN-corrected matrix, which contains 50 dimensions, and we use all of them (50 is also the default, so we could have left this option out).
The object now has a new object inside the graph slot, which a standard object type from the igraph package:
# the graph is stored as a standard igraph object
graph(milo)## IGRAPH 9d5cb1d U--- 30167 1035669 --
## + edges from 9d5cb1d:
## [1] 1-- 12 1-- 53 1-- 93 1-- 95 1-- 111 1-- 116 1-- 135 1-- 169 1-- 192 1-- 194 1-- 217 1-- 247 1-- 287 1-- 343 1-- 402 1-- 430 1-- 487 1-- 500 1-- 590 1-- 620 1-- 686 1-- 710 1-- 740 1-- 786
## [25] 1-- 912 1-- 962 1-- 1061 1-- 1066 1-- 1145 1-- 1151 1-- 1196 1-- 1213 1-- 1242 1-- 1368 1-- 1377 1-- 1402 1-- 1457 1-- 1487 1-- 1544 1-- 1553 1-- 1557 1-- 1574 1-- 1602 1-- 1641 1-- 1741 1-- 1760 1-- 1818 1-- 1893
## [49] 1-- 1930 1-- 1978 1-- 2002 1-- 2022 1-- 2063 1-- 2080 1-- 2241 1-- 2265 1-- 2347 1-- 2483 1-- 2507 1-- 2578 1-- 2654 1--20571 1--20652 1--25307 1--25366 1--26113 1--26814 1--27026 1--27331 1--27795 1--27928 1--28640
## [73] 1--29677 2-- 3 2-- 10 2-- 30 2-- 48 2-- 109 2-- 154 2-- 177 2-- 197 2-- 233 2-- 312 2-- 373 2-- 399 2-- 400 2-- 460 2-- 551 2-- 577 2-- 591 2-- 596 2-- 719 2-- 731 2-- 742 2-- 762 2-- 851
## [97] 2-- 898 2-- 941 2-- 1014 2-- 1083 2-- 1204 2-- 1347 2-- 1395 2-- 1411 2-- 1482 2-- 1558 2-- 1560 2-- 1687 2-- 1744 2-- 1756 2-- 1759 2-- 1874 2-- 1908 2-- 2035 2-- 2042 2-- 2079 2-- 2107 2-- 2112 2-- 2149 2-- 2191
## [121] 2-- 2196 2-- 2303 2-- 2315 2-- 2318 2-- 2368 2-- 2501 2-- 2552 2-- 2558 2-- 2622 2-- 2700 2-- 3854 2-- 3892 2-- 8760 2--24290 3-- 10 3-- 30 3-- 48 3-- 109 3-- 154 3-- 220 3-- 233 3-- 373 3-- 399 3-- 425
## [145] 3-- 460 3-- 533 3-- 596 3-- 731 3-- 851 3-- 1014 3-- 1026 3-- 1048 3-- 1482 3-- 1518 3-- 1560 3-- 1564 3-- 1687 3-- 1756 3-- 1759 3-- 1772 3-- 1874 3-- 1888 3-- 1895 3-- 1908 3-- 2035 3-- 2082 3-- 2106 3-- 2107
## [169] 3-- 2131 3-- 2191 3-- 2196 3-- 2303 3-- 2315 3-- 2318 3-- 2341 3-- 2393 3-- 2501 3-- 2552 3-- 2558 3-- 2622 3-- 2646 3-- 6268 3-- 9229 3--21021 3--23192 3--23194 3--23729 3--24290 3--24482 3--24594 3--25181 3--25879
## [193] 3--25917 4-- 9 4-- 75 4-- 115 4-- 228 4-- 323 4-- 337 4-- 365 4-- 370 4-- 409 4-- 437 4-- 572 4-- 574 4-- 597 4-- 599 4-- 677 4-- 727 4-- 767 4-- 815 4-- 835 4-- 922 4-- 924 4-- 1050 4-- 1060
## + ... omitted several edges11.2 Define neighbourhoods
The next step in the analysis is to define cell neighbourhoods. This essentially consists of picking a focal cell and counting how many other cells it is connected to in the graph (and how many come from each group under comparison). However, if we did this for every single cell in the data, it could get computationally quite intractable. Instead, Milo implements a sampling step, where so-called “index cells” are sampled from the larger graph, and those will be the neighbourhoods used in the DA analysis.
This is done with the makeNhodds() function:
# sample index cells to define neighbourhoods
milo <- makeNhoods(milo,
prop = 0.1,
k = 48,
d = 20,
reduced_dims = "HARMONY")
# check our object again
milo## class: Milo
## dim: 28343 30167
## metadata(0):
## assays(2): counts logcounts
## rownames(28343): ENSG00000238009 ENSG00000241860 ... FAM41AY1 FAM224B
## rowData names(0):
## colnames(30167): ETV6RUNX1-1_AAACCTGAGACTTTCG-1 ETV6RUNX1-1_AAACCTGGTCTTCAAG-1 ... PBMMC-3_TTTGTCATCAGTTGAC-1 PBMMC-3_TTTGTCATCTCGTTTA-1
## colData names(12): orig.ident nCount_RNA ... ident label
## reducedDimNames(3): PCA HARMONY UMAP
## mainExpName: SCT
## altExpNames(0):
## nhoods dimensions(2): 30167 2324
## nhoodCounts dimensions(2): 1 1
## nhoodDistances dimension(1): 0
## graph names(1): graph
## nhoodIndex names(1): 2324
## nhoodExpression dimension(2): 1 1
## nhoodReducedDim names(0):
## nhoodGraph names(0):
## nhoodAdjacency dimension(2): 1 1Some things to note:
propis the proportion of cells to sample (the authors advise 0.1 - 0.2).kanddshould be set the same as the k value used to build the original KNN graph. These values will be used for the graph-sampling algorithm used behind the scenes.
We can also see that our Milo object now has the nhoods slots populated. In this case, it indicates that the sampling algorithm picked 2324 index cells to form neighbourhoods for DA analysis.
One good QC metric to check at this stage is to create a histogram of cell counts per neighbourhood. Like we said earlier, when we set the value of k to define our KNN graph, we want to make sure the value is not too low, such that most neighbourhoods have very few cells, nor for it to be so high that we have very large (and presumably heterogenous) neighbourhoods.
Conveniently, MiloR has a plotting function just for this purpose:
# distribution of neighbourhood sizes
plotNhoodSizeHist(milo) +
geom_vline(xintercept = 100, col = "salmon")
The authors of Milo have stated several different parameters of deciding if your histogram is correct. Either ‘peaking above 20’, ‘peaking between 50 and 100’ or ‘an average neighbourhood size over 5 x N_samples’. Realistically, all of these statements give numbers in the same ballpark and so we can make a judgement on our data. In our case, we can see that our histogram peaks at ~100, suggesting a good neighbourhood size. If this was not the case, we could re-run the analysis from the start, making a KNN graph with a higher/lower k value.
11.3 Counting cells
After defining our sample of neighbourhoods, the next step consists of counting how many cells there are in each neighbourhood for each sample replicate. In this step we need to define which column of our metadata we want to use to do the counting. We should count cells at the biological replicate leve, which in our case is stored in the SampleName column.
# count cells in each neighbourhood
milo <- countCells(milo,
meta.data = colData(milo),
samples = "SampleName")
# Milo now has a counts matrix
head(nhoodCounts(milo))## 6 x 7 sparse Matrix of class "dgCMatrix"
## ETV6RUNX1-1 ETV6RUNX1-2 ETV6RUNX1-3 ETV6RUNX1-4 PBMMC-1 PBMMC-2 PBMMC-3
## 1 . 1 1 1 27 22 29
## 2 . 111 15 . . . .
## 3 23 19 27 21 . . .
## 4 . . . . 18 28 28
## 5 . 87 . . . . 2
## 6 . 5 . 11 38 21 21The dimensions of the counts matrix corresponds to the number of neighbourhoods (rows) and samples (columns), in our case a 2324 by 7 matrix.
11.4 Neighbourhood connectivity
There is one more step that we need to do before we are ready to run our DA analysis. It consists of calculating the distances between each neighbourhood. As we said, the neighbourhoods on our graph may partially overlap, so when Milo corrects our p-values for multiple testing, it takes into account the spatial dependence of our tests. For example, neighbourhoods that are closer to each other may have similar p-values, and so we should avoid penalising them too much, otherwise we will sacrifice statistical power.
# calculate distances between neighbourhoods - for p-value correction
# This step takes a few minutes
milo <- calcNhoodDistance(milo, d = 20, reduced.dim = "HARMONY")As before, the value of d (number of dimensions to consider from our MNN-corrected matrix) should be the same that was used for building our initial graph.
11.5 Running DA tests
Finally, we are ready to run the actual differential abundance step. Similarly to the pseudo-bulk appproach for differential expression, MiloR takes advantage of the edgeR package and its negative binomial linear model implementation. This provides a suitable statistical model to account for over-dispersed count data, as is common with these kind of datasets.
First, we need to define a simple table with information about our samples - we will use this table to define our model formula (similarly to what we did in the differential expression step). In this case, we want to detect DA between our PBMMC and ETV6-RUNX1 sample groups, so we define a table based on those two groups. We could also include a batch column in this table, but to keep the demonstration simple we will skip this. We would normally do this if we know there is a batch effect that we want to account for it in DA testing.
# define a table for our model design
sample_info <- unique(colData(milo)[, c("SampleName", "SampleGroup")])
rownames(sample_info) <- sample_info$SampleName
# order the table to match the order of samples in our counts matrix
sample_info <- sample_info[colnames(nhoodCounts(milo)), ]
sample_info## DataFrame with 7 rows and 2 columns
## SampleName SampleGroup
## <factor> <factor>
## ETV6RUNX1-1 ETV6RUNX1-1 ETV6RUNX1
## ETV6RUNX1-2 ETV6RUNX1-2 ETV6RUNX1
## ETV6RUNX1-3 ETV6RUNX1-3 ETV6RUNX1
## ETV6RUNX1-4 ETV6RUNX1-4 ETV6RUNX1
## PBMMC-1 PBMMC-1 PBMMC
## PBMMC-2 PBMMC-2 PBMMC
## PBMMC-3 PBMMC-3 PBMMCNow we can do the DA test, explicitly defining our experimental design. In this case as discussed we will test the difference between sample groups. The testNhoods function calculates a Fold-change and corrected P-value for each neighbourhood, which indicates whether there is significant differential abundance between sample groups.
# run DA test
da_results <- testNhoods(milo,
design = ~ SampleGroup,
design.df = sample_info,
reduced.dim = "HARMONY")
# results are returned as a data.frame
da_results %>%
arrange(SpatialFDR) %>%
head()## logFC logCPM F PValue FDR Nhood SpatialFDR
## 1 6.353066 9.682331 61.68553 7.590678e-05 0.004369078 1 0.003461454
## 15 8.784424 9.715095 103.00022 2.965902e-04 0.004369078 15 0.003461454
## 17 7.710021 9.878173 60.70449 9.336862e-05 0.004369078 17 0.003461454
## 25 6.160069 9.837182 53.69824 1.388821e-04 0.004369078 25 0.003461454
## 35 8.822261 9.750387 103.99192 2.905231e-04 0.004369078 35 0.003461454
## 37 7.841412 10.004035 64.41996 7.690225e-05 0.004369078 37 0.00346145411.6 Inspecting DA results
A good diagnostic plot to make after running our analysis is a histogram of p-values. We should expect a distribution with a peak close to 0 (corresponding to differentially abundant neighbourhoods) and tailing off towards 1. This blog article explains the different p-value histogram profiles you may see and what they can mean.
# p-value histogram
ggplot(da_results, aes(PValue)) +
geom_histogram(bins = 50)
Next, we can get an overview of our results as a volcano plot, marking a significance threshold of 10% FDR:
# volcano plot
# each point in this plot corresponds to a neighbourhood (not a cell)
ggplot(da_results, aes(logFC, -log10(SpatialFDR))) +
geom_point(aes(colour = FDR < 0.1)) +
geom_hline(yintercept = 1)
As we can see, several neighbourhoods fall below our FDR threshold, indicating significant differential abundance of cells between PBMMC and ETV6-RUNX1 samples.
There is an unsual discontinuity in our logFC axis, suggesting a sudden change in abundance in some of our neighbourhoods. We can get a better idea of the fold changes by visualising them as a graph for our neighbourhoods (rather than the original single-cell graph, which would be too big to display). We can then super-impose this graph on our original UMAP (or t-SNE), to compare with our original analysis.
# build neighbourhood graph embedding
milo <- buildNhoodGraph(milo)# our original UMAP with our previously annotated cell labels
umap_plot <- plotReducedDim(milo,
dimred = "UMAP",
colour_by = "label",
text_by = "label") +
theme(legend.position = "none")
# the neighbourhood map adjusted to match UMAP embedding
nh_graph_plot <- plotNhoodGraphDA(milo,
da_results,
res_column = "logFC",
layout = "UMAP",
alpha = 0.05)
# the two plots together side-by-side
umap_plot + nh_graph_plot +
plot_layout(guides="collect")
On the left we have our original UMAP, with cell/cluster annotations we did previously (by standard clustering and using known marker genes to manually annotate our cells). On the right we have the Milo neighbourhood graph, where each node represents a neighbourhood, coloured by the log fold-change in abundance between PBMMC and ETV6-RUNX1 samples (positive values represent higher abundance in ETV6-RUNX1, and vice-versa). We can see, for example, a cluster of cells with negative log fold changes around our large B cell cluster, which likely explains the discontinuity in values we saw earlier in our volcano plot.
Although we have annotations for our cells, these annotations at the moment are not present in the differential abundance table:
head(da_results)## logFC logCPM F PValue FDR Nhood SpatialFDR
## 1 6.353066 9.682331 61.685530 7.590678e-05 0.004369078 1 0.003461454
## 2 -7.494346 8.957664 3.951141 1.014005e-01 0.235307859 2 0.231165778
## 3 -7.063460 8.593517 17.084573 7.471359e-03 0.032334150 3 0.029750833
## 4 8.660067 9.599543 84.810279 4.502935e-04 0.004685037 4 0.003746897
## 5 -3.444498 8.604922 1.233463 3.024654e-01 0.464593214 5 0.449783234
## 6 4.200126 9.791865 8.975795 1.941770e-02 0.064191656 6 0.060118360We can transfer our cell labels to the neighbourhood DA results, by simply counting how many cells within that neighbourhood share the same label, and assign the label with the highest counts.
# annotate our neighbourhood DA results with our cell labels
da_results <- annotateNhoods(milo, da_results, coldata_col = "label")
head(da_results)## logFC logCPM F PValue FDR Nhood SpatialFDR label label_fraction
## 1 6.353066 9.682331 61.685530 7.590678e-05 0.004369078 1 0.003461454 NK T (c11) 1.0000000
## 2 -7.494346 8.957664 3.951141 1.014005e-01 0.235307859 2 0.231165778 B (c7) 0.9523810
## 3 -7.063460 8.593517 17.084573 7.471359e-03 0.032334150 3 0.029750833 B (c3) 1.0000000
## 4 8.660067 9.599543 84.810279 4.502935e-04 0.004685037 4 0.003746897 Erythrocytes (c17) 1.0000000
## 5 -3.444498 8.604922 1.233463 3.024654e-01 0.464593214 5 0.449783234 B (c1) 0.9325843
## 6 4.200126 9.791865 8.975795 1.941770e-02 0.064191656 6 0.060118360 CD20+ B (c10) 1.0000000The result includes the fraction of cells in the neighbourhood that shared that label. We can look at the distribution of this fraction as a quality check:
# histogram of fraction of cells in the neighbourhood with the same label
ggplot(da_results, aes(label_fraction)) +
geom_histogram(bins = 50)
We should expect the plot to be very biased towards 1, as is seen. This makes sense, since Milo’s neighbourhoods should be very similar to our previously-defined clusters. After all, both our clustering and Milo used KNN graphs as the starting point, and we’ve used the same settings for the k and d parameters in both steps of our analysis.
Despite most neighbourhoods being homogenous, some seem to have mixed labels. We can highlight these in our results table:
# add "mixed" label to neighbourhoods with less 70% consistency
da_results$label <- ifelse(da_results$label_fraction < 0.7,
"Mixed",
da_results$label)
head(da_results)## logFC logCPM F PValue FDR Nhood SpatialFDR label label_fraction
## 1 6.353066 9.682331 61.685530 7.590678e-05 0.004369078 1 0.003461454 NK T (c11) 1.0000000
## 2 -7.494346 8.957664 3.951141 1.014005e-01 0.235307859 2 0.231165778 B (c7) 0.9523810
## 3 -7.063460 8.593517 17.084573 7.471359e-03 0.032334150 3 0.029750833 B (c3) 1.0000000
## 4 8.660067 9.599543 84.810279 4.502935e-04 0.004685037 4 0.003746897 Erythrocytes (c17) 1.0000000
## 5 -3.444498 8.604922 1.233463 3.024654e-01 0.464593214 5 0.449783234 B (c1) 0.9325843
## 6 4.200126 9.791865 8.975795 1.941770e-02 0.064191656 6 0.060118360 CD20+ B (c10) 1.0000000Finally, we can visualise the distribution of DA fold changes in different labels.
# distribution of logFC across neighbourhood labels
plotDAbeeswarm(da_results, group.by = "label")
We can see that several clusters are enriched and some are depleted between our sample groups. There are not many DA neighbourhoods with a mixed label either. And we can see that B cells cluster 3 has quite a discontinuity in log fold changes - this is what was causing the odd distribution in the volcano plot, and we saw earlier is likely due to a sub-population of B cells with different representation across our samples.
It is also worth noting that this analysis may also reveal uncorrected batch effects, which would lead to an artificial separation of our samples according to their group. As always, we should always proceed in our analysis with care and back our conclusions with further confirmatory experiments.
There are further downstream analysis possible to identify genes that differentiate between our different neighbourhoods. This is similar in spirit to the marker gene analysis we did earlier, but based on the neighbourhood graph constructed by Milo. You can consult the package’s main vignette to see examples of these downstream analysis.
sessionInfo()
## R version 4.5.1 (2025-06-13)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 22.04.5 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0 LAPACK version 3.10.0
##
## locale:
## [1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8 LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8 LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C
## [10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Europe/London
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] patchwork_1.3.2 lubridate_1.9.5 forcats_1.0.1 stringr_1.6.0 dplyr_1.2.0 purrr_1.2.1 readr_2.1.6
## [8] tidyr_1.3.2 tibble_3.3.1 tidyverse_2.0.0 miloR_2.6.0 edgeR_4.8.0 limma_3.66.0 scater_1.38.0
## [15] ggplot2_4.0.2 scran_1.38.0 scuttle_1.20.0 SingleCellExperiment_1.32.0 SummarizedExperiment_1.40.0 Biobase_2.70.0 GenomicRanges_1.62.0
## [22] Seqinfo_1.0.0 IRanges_2.44.0 S4Vectors_0.48.0 BiocGenerics_0.56.0 generics_0.1.4 MatrixGenerics_1.22.0 matrixStats_1.5.0
## [29] BiocParallel_1.44.0
##
## loaded via a namespace (and not attached):
## [1] gridExtra_2.3 rlang_1.1.7 magrittr_2.0.4 compiler_4.5.1 vctrs_0.7.1 pkgconfig_2.0.3 fastmap_1.2.0 XVector_0.50.0 labeling_0.4.3 ggraph_2.2.2 rmarkdown_2.30
## [12] tzdb_0.5.0 pracma_2.4.6 ggbeeswarm_0.7.2 xfun_0.54 bluster_1.20.0 cachem_1.1.0 beachmat_2.26.0 jsonlite_2.0.0 DelayedArray_0.36.0 tweenr_2.0.3 irlba_2.3.5.1
## [23] parallel_4.5.1 cluster_2.1.8.1 R6_2.6.1 stringi_1.8.7 bslib_0.9.0 RColorBrewer_1.1-3 jquerylib_0.1.4 numDeriv_2016.8-1.1 Rcpp_1.1.0 knitr_1.50 splines_4.5.1
## [34] timechange_0.4.0 Matrix_1.7-4 igraph_2.2.1 tidyselect_1.2.1 dichromat_2.0-0.1 abind_1.4-8 yaml_2.3.10 viridis_0.6.5 codetools_0.2-20 lattice_0.22-7 withr_3.0.2
## [45] S7_0.2.1 evaluate_1.0.5 polyclip_1.10-7 pillar_1.11.1 rprojroot_2.1.1 hms_1.1.4 scales_1.4.0 gtools_3.9.5 glue_1.8.0 metapod_1.18.0 tools_4.5.1
## [56] BiocNeighbors_2.4.0 ScaledMatrix_1.18.0 locfit_1.5-9.12 graphlayouts_1.2.2 tidygraph_1.3.1 cowplot_1.2.0 grid_4.5.1 beeswarm_0.4.0 BiocSingular_1.26.0 ggforce_0.5.0 vipor_0.4.7
## [67] cli_3.6.5 rsvd_1.0.5 S4Arrays_1.10.0 viridisLite_0.4.2 gtable_0.3.6 sass_0.4.10 digest_0.6.38 SparseArray_1.10.1 ggrepel_0.9.6 dqrng_0.4.1 farver_2.1.2
## [78] memoise_2.0.1 htmltools_0.5.8.1 lifecycle_1.0.5 here_1.0.2 statmod_1.5.1 MASS_7.3-65