作者，Evil Genius
最近在看一些文章的時候，發(fā)現(xiàn)在計算細(xì)胞豐度差異的時候采用了miloR的方法，今天我們來關(guān)注一下。
moilR的原文如下，2022年就發(fā)表了。
與之相對的還有一個方法，文章在
其核心原理是通過構(gòu)建細(xì)胞鄰域（Neighborhood）來檢測實(shí)驗(yàn)條件或生物狀態(tài)對特定細(xì)胞群體分布的影響。算法核心原理是基于混合模型框架 和 K-最近鄰（KNN）圖的統(tǒng)計方法。
關(guān)鍵步驟
(1) 構(gòu)建 KNN 圖
基于降維后的細(xì)胞嵌入（如 PCA、UMAP），計算細(xì)胞間的歐氏距離，構(gòu)建 K-最近鄰圖（KNN graph）。
每個細(xì)胞的鄰域定義為與其最近的 k 個鄰居（默認(rèn) k=30），形成重疊的局部細(xì)胞群體。
(2) 定義鄰域（Neighborhoods）
對每個細(xì)胞，統(tǒng)計其所屬的鄰域（即作為其他細(xì)胞的鄰居的次數(shù)），形成鄰域-細(xì)胞矩陣。
通過 圖抽樣（graph sampling） 減少冗余，保留代表性鄰域。（ The first step in our method is to define a set of representative neighborhoods on the KNN graph, where a neighborhood is defined as the group of cells that are connected to an index cell by an edge in the graph. ）
(3) 差異豐度檢驗(yàn)。對每個鄰域中存在的細(xì)胞數(shù)量進(jìn)行計數(shù)（每個實(shí)驗(yàn)樣品），并將其用于條件之間的差異豐度測試。
廣義線性混合模型（GLMM）：
對每個鄰域，擬合一個混合模型（如負(fù)二項(xiàng)分布或泊松分布），以實(shí)驗(yàn)條件（如分組）作為固定效應(yīng)，樣本（個體）作為隨機(jī)效應(yīng)，校正批次效應(yīng)或個體間差異。
counts ~ condition + (1 | sample)
其中 counts 是鄰域中的細(xì)胞數(shù)，condition 是分組變量，sample 是樣本來源。
假設(shè)檢驗(yàn)：
通過似然比檢驗(yàn)（LRT）或 Wald 檢驗(yàn)比較模型（有/無實(shí)驗(yàn)條件效應(yīng)），得到差異豐度的 p 值。
(4) 多重檢驗(yàn)校正
使用 FDR（如 Benjamini-Hochberg）校正鄰域水平的 p 值，控制假陽性率。
(5) 結(jié)果可視化
通過 UMAP 或降維圖標(biāo)記顯著差異的鄰域，結(jié)合火山圖展示 logFC 和顯著性。
把空間的鄰域分析思想帶到了單細(xì)胞組學(xué)。
最后來看看分析代碼
library(miloR)
library(SingleCellExperiment)
library(scater)
library(scran)
library(dplyr)
library(patchwork)
library(Seurat)
pbmc_small = readRDS('test.rds')
pbmc_small_sce <- as.SingleCellExperiment(pbmc_small)
pbmc_small_milo <- Milo(pbmc_small_sce)
####Construct KNN graph
traj_milo <- buildGraph(traj_milo, k = 10, d = 30)
####Defining representative neighbourhoods
traj_milo <- makeNhoods(traj_milo, prop = 0.1, k = 10, d=30, refined = TRUE)
####Counting cells in neighbourhoods
traj_milo <- countCells(traj_milo, meta.data = data.frame(colData(traj_milo))
####Differential abundance testing
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_milo <- calcNhoodDistance(traj_milo, d=30)
rownames(traj_design) <- traj_design$Sample
da_results <- testNhoods(traj_milo, design = ~ Condition, design.df = traj_design)
da_results %>%
arrange(- SpatialFDR) %>%
head()
##Visualize neighbourhoods displaying DA
traj_milo <- buildNhoodGraph(traj_milo)
plotUMAP(traj_milo) + plotNhoodGraphDA(traj_milo, da_results, alpha=0.05) +
plot_layout(guides="collect")
