From bf5c71246c8b972a000ad16a7bad9e5a1044c32f Mon Sep 17 00:00:00 2001 From: ixxmu Date: Sun, 24 Nov 2024 01:33:09 +0000 Subject: [PATCH] =?UTF-8?q?miloR=E5=8D=95=E7=BB=86=E8=83=9E=E5=B7=AE?= =?UTF-8?q?=E5=BC=82=E4=B8=B0=E5=BA=A6=E5=88=86=E6=9E=90?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- ...70\260\345\272\246\345\210\206\346\236\220.md" | 15 +++++++++++++++ 1 file changed, 15 insertions(+) create mode 100644 "docs/2024-11/miloR\345\215\225\347\273\206\350\203\236\345\267\256\345\274\202\344\270\260\345\272\246\345\210\206\346\236\220.md" diff --git "a/docs/2024-11/miloR\345\215\225\347\273\206\350\203\236\345\267\256\345\274\202\344\270\260\345\272\246\345\210\206\346\236\220.md" "b/docs/2024-11/miloR\345\215\225\347\273\206\350\203\236\345\267\256\345\274\202\344\270\260\345\272\246\345\210\206\346\236\220.md" new file mode 100644 index 00000000..9982ef09 --- /dev/null +++ "b/docs/2024-11/miloR\345\215\225\347\273\206\350\203\236\345\267\256\345\274\202\344\270\260\345\272\246\345\210\206\346\236\220.md" @@ -0,0 +1,15 @@ +--- +title: "miloR单细胞差异丰度分析" +date: 2024-11-24T01:32:24Z +draft: ["false"] +tags: [ + "fetched", + "生信随笔" +] +categories: ["Acdemic"] +--- +miloR单细胞差异丰度分析 by 生信随笔 +------ +

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算法的工作流程示意图：

image-20241115210539230

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") 
  1. Basic Milo example on simulated dataset
  2. Milo example on mouse gastrulation dataset: this includes a demo for downstream analysis functions.
  3. Integrating miloR in scanpy/anndata workflow (see also milopy for a full workflow in python)
  4. Specifying contrasts of interest for differential abundance testing with Milo
  5. 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)
image-20241016151136405

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)
image-20241016151246258

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")
image-20241115212516767

官网的这个教程到这里就结束了。其实第二个教程里还有一些实用的分析及可视化代码，可以将这些结果映射到我们感兴趣的细胞类型上，详见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")
image-20241115213209259


+
+原文链接