miloR发表在2021年Nat Biotechnol上，题为《Differential abundance testing on single-cell data using k-nearest neighbor graphs》。简单来说，miloR是一个用于单细胞数据分析的 R 包，专注于细胞间的差异群体丰度（Differential Abundance, DA）分析。它通过构建细胞邻域（Neighborhood）来检测实验条件或生物状态对特定细胞群体分布的影响。miloR提供了一种统计框架，可以从单细胞测序数据中提取新的生物学见解，特别适合研究复杂的生物学体系，如发育、疾病或治疗响应。
一. miloR算法的工作流程
尽管差异丰度（Differential abundance, DA）通常在离散细胞簇中进行量化，但 MiloR 使用 KNN 图中的部分重叠细胞邻域。从一个真实反映细胞群体生物学的图开始，Milo 分析包括以下三个步骤：
采样代表性邻域； 测试所有邻域中条件的差异丰度； 使用加权假发现率（FDR）程序考虑多个假设检验，该程序考虑邻域的重叠。
MiloR算法的工作流程示意图：
a. 通过一种图采样算法在索引细胞上定义邻域。根据实验设计对细胞进行定量生成计数表。每个邻域的细胞计数使用负二项广义线性模型（GLM）进行建模，并通过假设检验确定差异丰度的邻域。
b. KNN 图的力导向布局，展示从两种实验条件中采样的细胞的模拟连续轨迹。
c. 使用 Milo 进行假设检验能够准确且特异性地检测出差异丰度的邻域（假发现率 FDR 1%）。红点表示差异丰度的邻域。
d. Milo 差异丰度检验结果的图形表示。节点代表邻域，颜色根据其对数倍数变化（log fold change）着色。非差异丰度的邻域（FDR 1%）以白色表示，节点大小与邻域中的细胞数量成比例。图中的边缘表示相邻邻域之间共享的细胞数量。节点布局由力导向嵌入中的邻域索引细胞的位置决定。Nhood 代表邻域。
二. miloR代码实现
miloR官方github在https://github.com/MarioniLab/miloR
教程在：https://bioconductor.org/packages/release/bioc/vignettes/miloR/inst/doc/milo_demo.html
R包安装：
## Milo is available from Bioconductor (preferred stable installation)
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("miloR")
## Install development version
devtools::install_github("MarioniLab/miloR", ref="devel")
Basic Milo example on simulated dataset Milo example on mouse gastrulation dataset: this includes a demo for downstream analysis functions. Integrating miloR in scanpy/anndata workflow (see also milopyfor a full workflow in python)Specifying contrasts of interest for differential abundance testing with Milo Using a mixed effect model for dependendent samples
github上有五个教程。这里我们先看一下第一个基础教程。
Step1. 加载示例数据
#BiocManager::install("miloR")
#BiocManager::install("zellkonverter")
library(miloR)
library(SingleCellExperiment)
library(scater)
library(scran)
library(dplyr)
library(patchwork)
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)
Step2. 数据预处理
为了进行数据分析，我们需要构造一个单细胞的无向KNN图。下面是miloR提供的细胞标准分析代码：
logcounts(traj_sce) <- log(counts(traj_sce) + 1)
traj_sce <- runPCA(traj_sce, ncomponents=30)
traj_sce <- runUMAP(traj_sce)
plotUMAP(traj_sce)
Step3. 构建MiloR对象
对于图邻域的差分丰度分析，我们首先构造一个Milo对象。这扩展了singlecelleexperiment类来存储KNN图上的邻域信息。下面提供了不同单细胞数据格式的构建代码：
3.1 From SingleCellExperiment object
traj_milo <- Milo(traj_sce)
reducedDim(traj_milo, "UMAP") <- reducedDim(traj_sce, "UMAP")
traj_milo
3.2 From AnnData object (.h5ad)
我们可以使用zellkonverter包从保存为h5ad文件的AnnData对象中生成singlecellexexperiment对象。
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)
3.3 From Seurat object
library(Seurat)
data("pbmc_small")
pbmc_small_sce <- as.SingleCellExperiment(pbmc_small)
pbmc_small_milo <- Milo(pbmc_small_sce)
Step4. 构建KNN图
我们需要将KNN图添加到Milo对象中。它以图的格式存储在graph插槽中
traj_milo <- buildGraph(traj_milo, k = 10, d = 30)
Step5. 定义代表性邻域（neighbourhoods）
作者将一个细胞的邻域定义为与其在 KNN 图中通过边相连的细胞群，这个细胞被称为索引细胞（index cell）。为了提高效率，我们不会对每个细胞的邻域进行差异丰度（DA）检验，而是通过一种 KNN 采样算法（参见 Gut et al. 2015）从中选取一部分代表性细胞作为索引进行检验。
在采样过程中，需要定义以下参数：
prop: 随机采样的细胞比例（通常 0.1 - 0.2 即可）。对于少于 30,000 个细胞的数据集，建议使用prop = 0.1；对于更大的数据集，可以选择prop = 0.05，这样计算速度更快。k: KNN 细化所使用的 k 值（建议与 KNN 图构建时使用的 k 值保持一致）。d: KNN 细化所使用的降维数（建议与 KNN 图构建时使用的降维数保持一致）。refined: 指定是否使用采样细化算法，或者仅随机选择细胞。默认推荐使用细化算法。只有在数据经过基于图的批次校正算法（如 BBKNN）处理后，可能会考虑改用随机采样，但这种方法会导致 DA 检验结果次优。refinement_scheme和reduced_dims: 指定使用 PCA 还是图（graph）作为细化方案，默认选择 PCA。
这些参数可以帮助优化采样过程，提高代表性细胞选择的准确性，同时平衡计算效率和生物学意义。
traj_milo <- makeNhoods(traj_milo, prop = 0.1, k = 10, d=30, refined = TRUE)
一旦我们定义了neighbourhood，最好看看neighbourhood有多大（即每个neighbourhood有多少个细胞）。这会影响数据分析测试的能力。我们可以使用plotNhoodSizeHist函数检查这一点。根据经验，我们希望拥有超过5 x N_samples的平均邻域大小。否则，您可能会考虑重新运行makeNhoods，增加k和/或prop。
plotNhoodSizeHist(traj_milo)
Step7. 计算领域中的细胞数量
现在我们要计算每个样本在每个邻域中有多少个细胞。我们需要使用细胞元数据并指定哪一列包含示例信息。
traj_milo <- countCells(traj_milo, meta.data = data.frame(colData(traj_milo)), samples="Sample")
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 7 5 12 8 11
## 2 13 2 1 19 21 1
## 3 10 1 1 4 2 2
## 4 8 . . 14 28 1
## 5 4 6 6 6 7 13
## 6 4 4 4 8 9 6
这将给Milo对象添加一个n \ m矩阵，其中n是邻域的数量，m为实验样本个数。数值表示在一个邻域中计算的每个样本的细胞数。该计数矩阵将用于DA测试。
Step8. 差异丰度检验
现在我们都准备好测试不同neighbourhoods的丰度差异。我们在广义线性模型（GLM）框架中实现了这个假设检验，特别是在edgeR中使用负二项GLM实现。
我们首先需要考虑我们的实验设计。设计矩阵应使样品与感兴趣的条件相匹配。在这种情况下，"Condition"是我们要测试的协变量。
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使用了cydar引入的空间FDR校正的适应性，它解释了邻里之间的重叠。具体来说，每个假设检验的p值由第k个最近邻距离的倒数加权。要使用这个统计数据，我们首先需要在Milo对象中存储最近邻居之间的距离。
traj_milo <- calcNhoodDistance(traj_milo, d=30)
现在我们可以做检验分析，明确定义我们的实验设计。
rownames(traj_design) <- traj_design$Sample
da_results <- testNhoods(traj_milo, design = ~ Condition, design.df = traj_design)
这计算了每个邻域的Fold-change和修正的p值，这表明条件之间是否存在显著的丰度差异。
da_results %>%
arrange(- SpatialFDR) %>%
head()
## logFC logCPM F PValue FDR Nhood SpatialFDR
## 21 -0.007155451 15.06606 0.0001967759 0.9897851 0.9897851 21 0.9897851
## 30 -0.069989274 14.82502 0.0111343047 0.9161290 0.9456816 30 0.9455256
## 6 -0.184051563 15.33523 0.1301442022 0.7188772 0.7668024 6 0.7637508
## 23 -0.301800195 14.73485 0.1838471716 0.6688084 0.7379955 23 0.7341794
## 10 -0.281676336 15.36561 0.3816259156 0.5580706 0.6377950 10 0.6323192
## 28 -0.348771823 15.10259 0.4158152581 0.5201866 0.6165175 28 0.6106588
Step9. 可视化
为了可视化差异丰度（DA）结果，我们构建了一个邻域的抽象图，并将其叠加到单细胞的嵌入图上：
traj_milo <- buildNhoodGraph(traj_milo)
plotUMAP(traj_milo) + plotNhoodGraphDA(traj_milo, da_results, alpha=0.05) +
plot_layout(guides="collect")
官网的这个教程到这里就结束了。其实第二个教程里还有一些实用的分析及可视化代码，可以将这些结果映射到我们感兴趣的细胞类型上，详见https://rawcdn.githack.com/MarioniLab/miloR/7c7f906b94a73e62e36e095ddb3e3567b414144e/vignettes/milo_gastrulation.html#5_Finding_markers_of_DA_populations：
traj_milo <- buildNhoodGraph(traj_milo)
## Plot single-cell UMAP
umap_pl <- plotReducedDim(traj_milo, dimred = "umap", colour_by="stage", text_by = "celltype",
text_size = 3, point_size=0.5) +
guides(fill="none")
## Plot neighbourhood graph
nh_graph_pl <- plotNhoodGraphDA(traj_milo, da_results, layout="umap",alpha=0.1)
umap_pl + nh_graph_pl +
plot_layout(guides="collect")
da_results <- annotateNhoods(traj_milo, da_results, coldata_col = "celltype")
head(da_results)
## logFC logCPM F PValue FDR Nhood SpatialFDR
## 1 -1.8019342 20.07089 1.5747061 0.209611035 0.29435477 1 0.29385683
## 2 4.7439141 19.12468 7.0132433 0.008128239 0.03458566 2 0.03452842
## 3 5.4795871 19.87780 8.6993418 0.003204932 0.02281712 3 0.02248097
## 4 -1.6502089 20.24398 1.2210504 0.269232708 0.36142976 4 0.36054817
## 5 3.2520711 20.61980 3.6364192 0.056612494 0.11608972 5 0.11588354
## 6 0.8336509 19.82476 0.2813193 0.595872418 0.68340624 6 0.68348253
## celltype celltype_fraction
## 1 ExE endoderm 1.0000000
## 2 Rostral neurectoderm 0.5581395
## 3 ExE endoderm 1.0000000
## 4 ExE ectoderm 1.0000000
## 5 ExE ectoderm 1.0000000
## 6 ExE ectoderm 1.0000000
da_results$celltype <- ifelse(da_results$celltype_fraction < 0.7, "Mixed", da_results$celltype)
plotDAbeeswarm(da_results, group.by = "celltype")
