R tips：空间最近邻分析

生信菜鸟团

发布于 2025-05-12 10:48:44

2170

Xenium或者VisiumHD等空间组学数据可以进行很多不同于单细胞的空间分析，比如本文提到的最近邻分析，它可以根据坐标判断一个细胞类型周围的最近邻细胞的分布趋势和在不同组别中的对比。

示例数据如下的meta_data所示，正式分析时需要先将数据整合成meta_data所示的格式。

library(tidyverse)
library(Seurat)

meta_data %>% head()
# cell_id   Group Slide_ID        X        Y CellType
# 1 3_aaaacdgo-1 Group_A  Slide_1 2228.621 860.3357   Cell_I
# 2 3_aaaanalm-1 Group_A  Slide_1 2248.922 850.1516   Cell_K
# 3 3_aaaaonpl-1 Group_A  Slide_1 2270.447 856.8892   Cell_K
# 4 3_aaabciio-1 Group_A  Slide_1 2279.562 838.6568   Cell_B
# 5 3_aaabdkea-1 Group_A  Slide_1 2292.319 858.3909   Cell_I
# 6 3_aaacdgip-1 Group_A  Slide_1 2288.924 868.3226   Cell_B

cell_id是数据中的唯一细胞id，尽量保证唯一即可，在最近邻分析时不太重要。

Group是分析样本的组别。

Slide_ID是代表多个芯片的标识，每个芯片上可以有多个样本。由于空间数据的坐标是以芯片或者FOV为一个单元进行标识的，而且从生物学意义上来说，最近邻的空间坐标也只在每个样本内进行计算才有意义（毕竟就算人为把多个样本贴的再近，也不代表任何生物学意义），因此这里的Slide_ID一般是指代芯片ID或者是样本ID，一般优选是样本ID。

X,Y是每个细胞的空间坐标。

CellType是注释号的细胞类型。

使用RANN::nn2进行最近邻计算

由于多个芯片的数据整合在一起了，如上面解释的，每个芯片或样本需要先单独进行最近邻细胞的识别。

# 先拆分成多个芯片的数据
meta_data_ls <-
  meta_data %>%
  group_by(Slide_ID) %>%
  nest() %>%
  deframe()

使用RANN::nn2进行最近邻计算：

res_nn_ls <-
  meta_data_ls %>%
  map(function(dat){
    res <-
      RANN::nn2(
        data  = dat[c("X","Y")],# 去那些地方寻找最近邻细胞
        query = dat[c("X","Y")],# 这些细胞需要寻找最近邻细胞
        k =10# 寻找多个最近邻细胞
      )
    res
  })

注意，RANN::nn2的返回值是一个list，包括nn.idx和nn.dist两个元素：

1. nn.idx是一个矩阵，行是query细胞的数量，列是最近邻细胞的数量，根据距离的远近由近到远依次排列，矩阵的值是最近邻细胞的index。比如Slide_1的nn.idx：

res_nn_ls$Slide_1$nn.idx %>% head
# [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10]
# [1,]    1     2 16356 30348   101 16364 30355 16354     3 16345
# [2,]    2 16364     3 16356     1 16354 30348 30355 16345 16390
# [3,]    3 30355     4 16364     6     5     2 30360 30348     7
# [4,]    4 16357 16364     3 16390     7     5 16354 15735     6
# [5,]    5     6     7     3     4 30360 30378 30355 30365 15735
# [6,]    6     5 30360 30365     3 30355     7 30378     4 30361

第一行就是第一个query细胞，他的最近邻细胞的index是1 2 16356 30348 101 16364 30355 16354 3 16345，那么这10个最近邻细胞就是：

meta_data_ls$Slide_1[res_nn_ls$Slide_1$nn.idx[1, ], ]
# A tibble: 10 × 6
# cell_id      Classification1  Group       X     Y CellType
# <chr>        <chr>            <chr>   <dbl> <dbl> <chr>   
#   1 3_aaaacdgo-1 Fibroblasts      Group_A 2229.  860. Cell_I  
# 2 3_aaaanalm-1 Myeloid cells    Group_A 2249.  850. Cell_K  
# 3 3_bdbjmdoo-1 Fibroblasts      Group_A 2235.  832. Cell_I  
# 4 3_cdbibgio-1 Granulocytes     Group_A 2254.  875. Cell_J  
# 5 3_aaccanjd-1 Fibroblasts      Group_A 2201.  843. Cell_I  
# 6 3_bdbmlcml-1 Fibroblasts      Group_A 2260.  838. Cell_I  
# 7 3_cdbjipci-1 CD4Ts            Group_A 2267.  871. Cell_C  
# 8 3_bdbjhjdp-1 Fibroblasts      Group_A 2253.  827. Cell_I  
# 9 3_aaaaonpl-1 Myeloid cells    Group_A 2270.  857. Cell_K  
# 10 3_bdbhehpa-1 Epithelial cells Group_A 2246.  822. Cell_H

另，由于进行RANN::nn2计算时的data和query是同一个数据，因此nn.idx的第一列（也就是最近的那个细胞）一定是其自身，也那么nn.idx的第一列的值一定是1,2,3,4,...按顺序排序的值，而真正需要关注的最近邻细胞应该是nn.idx的第二列。

1. nn.dist和nn.idx结构一样，不过其值是最近邻细胞的距离值。同样的道理，其第一列的值一定是0。

res_nn_ls$Slide_1$nn.dist %>% head()
# [,1]     [,2]     [,3]     [,4]     [,5]     [,6]     [,7]     [,8]     [,9]    [,10]
# [1,]    0 22.71249 29.00059 29.14782 32.56646 38.22030 40.13437 41.47847 41.96744 42.17485
# [2,]    0 16.12265 22.55427 22.63088 22.71249 23.77691 25.10910 27.80881 28.50026 29.06299
# [3,]    0 14.47012 20.38398 21.40854 21.72842 21.92357 22.55427 23.55163 24.31147 31.65161
# [4,]    0 19.68257 19.72214 20.38398 20.58600 20.76804 23.49836 29.20852 29.83792 31.10790
# [5,]    0 10.49589 15.74977 21.92357 23.49836 25.56396 27.69624 28.01566 30.04889 34.96368
# [6,]    0 10.49589 15.75802 20.36185 21.72842 21.78226 26.17876 27.97772 31.10790 32.90458

数据清洗和绘图展示

数据清洗，拿到真正的最近邻细胞的类型以及相应的组别及距离。

res_nn_df_ls <-
  res_nn_ls %>%
  imap(function(res_nn, slide_name){
    meta <- meta_data_ls[[slide_name]]

  # celltype
      center_cell <- meta$CellType
      nn_cell <- meta$CellType[res_nn$nn.idx[,2]]# 第二列是最近邻细胞

  # distacne
      dist_val <- res_nn$nn.dists[,2]

  # return
    data.frame(
      center = center_cell,
      others = nn_cell,
      distance = dist_val,
      Group= meta$Group
    )
  })

# 多个芯片的数据可以合并在一起了
res_nn_df <- bind_rows(res_nn_df_ls)

res_nn_df %>% head
# center others distance   Group
# 1 Cell_I Cell_K 22.71249 Group_A
# 2 Cell_K Cell_I 16.12265 Group_A
# 3 Cell_K Cell_C 14.47012 Group_A
# 4 Cell_B Cell_K 19.68257 Group_A
# 5 Cell_I Cell_B 10.49589 Group_A
# 6 Cell_B Cell_I 10.49589 Group_A

用于进一步绘图的数据是res_nn_df，其中center是中心细胞，others是中心细胞的最近邻细胞，distance是两个细胞的距离，Group是样本的组别。

密度图

仅仅以Cell_A为例，查看其最近邻细胞的分布情况：

res_nn_df %>%
  dplyr::filter(center == "Cell_A") %>% 
  ggplot(aes(x = distance, fill = Group)) +
  geom_density(alpha = 0.5) +
  facet_wrap(~ others) +
  scale_x_log10()

res_nn_df %>%
  dplyr::filter(center == "Cell_A") %>% 
  ggplot(aes(x = distance, fill = others)) +
  geom_density(alpha = 0.5) +
  facet_wrap(~ Group) +
  scale_x_log10()

堆叠密度图

仅仅以Cell_A为例，查看其最近邻细胞的分布情况：

res_nn_df %>%
  dplyr::filter(center =="Cell_A")%>%
  ggplot(aes(x = distance, y =Group, fill = after_stat(x)))+
  ggridges::geom_density_ridges_gradient()+
  facet_wrap(~ others)+
  ggridges::theme_ridges()+
  scale_x_log10()+
  scale_fill_viridis_c(option ="C", transform ="log10")

res_nn_df %>%
  dplyr::filter(center =="Cell_A")%>%
  ggplot(aes(x = distance, y = others, fill = after_stat(x)))+
  ggridges::geom_density_ridges_gradient()+
  facet_wrap(~Group)+
  ggridges::theme_ridges()+
  scale_x_log10()+
  scale_fill_viridis_c(option ="C", transform ="log10")