1. 数据准备

library(miloR)
library(SingleCellExperiment)
library(scater)
library(scran)
library(dplyr)
library(patchwork)
 
#Load data
load(file = "./sc.RData")
sce <- as.SingleCellExperiment(sc)
 
#Visualize the data
plotReducedDim(sce, colour_by="celltype", dimred = "UMAP") 

2. 构建领域，对领域中的细胞计数
   d：用于KNN细化的降维数，类似于聚类降维时选择的dims。
   k：分布峰值在 50 到 100 之间，平均邻域大小超过 5 x N_samples。

#Create a Milo object
scmilo <- Milo(sce)
 
#Construct KNN graph
scmilo <- buildGraph(scmilo, k = 30, d = 30,reduced.dim = "PCA")
 
#Defining representative neighbourhoods on the KNN graph
set.seed(10)
scmilo <- makeNhoods(scmilo, prop = 0.1,
                     k = 30, d = 30, 
                     refined = TRUE, reduced_dims = "PCA")
plotNhoodSizeHist(scmilo)
 
#Counting cells in neighbourhoods
scmilo <- countCells(scmilo, 
                     meta.data = as.data.frame(colData(scmilo)), 
                     sample = "sample")
head(nhoodCounts(scmilo))
## 6 x 6 sparse Matrix of class "dgCMatrix"
##   2  3  6 4 10 14
## 1 .  1 14 1 99 20
## 2 7 10 35 4  2  4
## 3 3  4 29 4 26  6
## 4 5 17 37 7  .  1
## 5 .  .  3 . 33 12
## 6 .  1  8 3 69 11

3. 定义分组，差异测试

#Defining experimental design
sc_design <- data.frame(colData(scmilo))[,c("sample", "group")]
 
sc_design <- distinct(sc_design)
rownames(sc_design) <- sc_design$sample
 
#Computing neighbourhood connectivity
scmilo <- calcNhoodDistance(scmilo, d = 30, reduced.dim = "PCA")
 
#Testing
results <- testNhoods(scmilo, design = ~ group, 
                      design.df = sc_design, reduced.dim="PCA")
head(results)
 
results %>%
  arrange(SpatialFDR) %>%
  head() 

4.可视化

#Inspecting DA testing results
ggplot(results, aes(PValue)) + geom_histogram(bins=50)
ggplot(results, aes(logFC, -log10(SpatialFDR))) + 
  geom_point() +
  geom_hline(yintercept = 1) 

scmilo <- buildNhoodGraph(scmilo)
 
## Plot single-cell UMAP
umap_pl <- plotReducedDim(scmilo, dimred = "UMAP", 
                          colour_by="celltype", text_by = "celltype", 
                          text_size = 3, point_size=0.5) +
  guides(fill="none")
 
## Plot neighbourhood graph
nh_graph_pl <- plotNhoodGraphDA(scmilo, results, layout="UMAP",alpha = 1) 
 
umap_pl + nh_graph_pl +
  plot_layout(guides="collect")

results <- annotateNhoods(scmilo, results, coldata_col = "celltype")
head(results)
 
plotDAbeeswarm(results, group.by = "celltype",alpha = 1)
#alpha值调高可以看所有的结果