课前准备--Stereo-seq（HD）细胞类型的分布梯度

copykat.test <- copykat(rawmat=exp.rawdata, 
                        id.type="S", # S是symbol的含义
                        ngene.chr=5, # 规定每条染色体中至少5个基因来计算DNA拷贝数(可调整)
                        win.size=25, # 每个片段至少25个基因
                        KS.cut=0.1, # 增加KS.cut会降低敏感度，通常范围在0.05-0.15
                        sam.name="241016",  # 自行命名
                        distance="euclidean", 
                        norm.cell.names="",
                        output.seg="FLASE", 
                        plot.genes="TRUE", 
                        genome="hg20",
                        n.cores=8)

library(RcppAnnoy)
library(future.apply)
library(Seurat)
library(ggplot2)
library(dplyr)

BuildTree = function(data, n_trees = 50){
  vector_size = ncol(data)
  a_tree <- new(AnnoyEuclidean, vector_size)
  n_samples = nrow(data)
  for (i in 1:n_samples){
    v = as.vector(as.matrix(data[i,]))
    a_tree$addItem(i - 1, v)
  }
  a_tree$build(n_trees)
  return(a_tree)
}

# ========== 2. 获取最近邻函数 ==========
GetNearestNeighbor = function(seurat_obj, n_trees = 50, n_neigh = 10) {
  data = seurat_obj@meta.data[, c("x", "y")]
  a_tree = BuildTree(data = data, n_trees = n_trees)
  neigh_idx = matrix(nrow = nrow(data), ncol = n_neigh)
  neigh_dist = matrix(nrow = nrow(data), ncol = n_neigh)
  res = future_lapply(1:nrow(data), function(x) {
    nei_list = a_tree$getNNsByItemList(x-1, n_neigh, -1, TRUE)
    nei_list$item = nei_list$item + 1
    return(nei_list)
  })
  for (i in 1:nrow(data)) {
    neigh_idx[i, ] <- res[[i]][[1]]
    neigh_dist[i, ] <- res[[i]][[2]]
  }
  return(list(nn.idx = neigh_idx, nn.dists = neigh_dist))
}


seurat_obj_cell@neighbors$Spatial <- GetNearestNeighbor(seurat_obj = seurat_obj_cell, n_neigh = 50)

# 初始定义：D4区域为Tumor core，其余为Stroma(示例)
tumor_domain_cell_idx <- which(seurat_obj_cell@meta.data[, 'domain'] == "D4")
seurat_obj_cell$Tumor_boundary <- "Stroma"
seurat_obj_cell$Tumor_boundary[tumor_domain_cell_idx] <- "Tumor core"

# ========== 从肿瘤核心向外扩张，定义Outer层 ==========
target_cell_idx <- tumor_domain_cell_idx
target_cell_domain <- "Tumor core"
min_neigh <- 10

for (i in 1:3){
  print(paste("Defining Outer", i))
  target_cell_dist <- seurat_obj_cell@neighbors$Spatial$nn.dists[target_cell_idx, ]
  target_cell_neigh <- seurat_obj_cell@neighbors$Spatial$nn.idx[target_cell_idx, ]
  
  # 找到距离≤150μm且满足条件的Stroma细胞
  target_cell_dist_200 <- sapply(1:nrow(target_cell_neigh), function(x){
    idx <- which(target_cell_dist[x, ] <= 150)
    neigh <- unlist(target_cell_neigh[x, idx])
    neigh_domain <- seurat_obj_cell@meta.data[neigh, "Tumor_boundary"]
    if (length(which(neigh_domain == "Stroma")) > 0 && 
        length(which(neigh_domain == target_cell_domain)) > min_neigh){
      return(neigh[neigh_domain == "Stroma"])
    } else{
      return(NA)
    }
  })
  
  target_cell_dist_200_neigh <- unique(unlist(target_cell_dist_200))
  target_cell_dist_200_neigh <- target_cell_dist_200_neigh[!is.na(target_cell_dist_200_neigh)]
  
  seurat_obj_cell$Tumor_boundary[target_cell_dist_200_neigh] <- paste0("Outer ", i)
  target_cell_idx <- target_cell_dist_200_neigh
  target_cell_domain <- paste0("Outer ", i)
}

# ========== 从Outer 1向内扩张，定义Inner层 ==========
o1_domain_cell_idx <- which(seurat_obj_cell@meta.data[, 'Tumor_boundary'] == "Outer 1")
target_cell_idx <- o1_domain_cell_idx
target_cell_domain <- "Outer 1"

for (i in 1:3){
  print(paste("Defining Inner", i))
  target_cell_dist <- seurat_obj_cell@neighbors$Spatial$nn.dists[target_cell_idx, ]
  target_cell_neigh <- seurat_obj_cell@neighbors$Spatial$nn.idx[target_cell_idx, ]
  
  # 找到距离≤150μm且满足条件的Tumor core细胞
  target_cell_dist_200 <- sapply(1:nrow(target_cell_neigh), function(x){
    idx <- which(target_cell_dist[x, ] <= 150)
    neigh <- unlist(target_cell_neigh[x, idx])
    neigh_domain <- seurat_obj_cell@meta.data[neigh, "Tumor_boundary"]
    if (length(which(neigh_domain == "Tumor core")) > 0 && 
        length(which(neigh_domain == target_cell_domain)) > min_neigh){
      return(neigh[neigh_domain == "Tumor core"])
    } else{
      return(NA)
    }
  })
  
  target_cell_dist_200_neigh <- unique(unlist(target_cell_dist_200))
  target_cell_dist_200_neigh <- target_cell_dist_200_neigh[!is.na(target_cell_dist_200_neigh)]
  
  seurat_obj_cell$Tumor_boundary[target_cell_dist_200_neigh] <- paste0("Inner ", i)
  target_cell_idx <- target_cell_dist_200_neigh
  target_cell_domain <- paste0("Inner ", i)
}

# ========== 设置因子水平 ==========
seurat_obj_cell$Tumor_boundary <- factor(seurat_obj_cell$Tumor_boundary, 
  levels = c("Tumor core", "Inner 3", "Inner 2", "Inner 1", 
             "Outer 1", "Outer 2", "Outer 3", "Stroma"))

color_boundary <- c("#A71B4BFF", "#E96F02FF", "#F9C25CFF", "#FEFDBEFF", 
                    "#81DEADFF", "#00A3B6FF", "#584B9F", "#e9e9e9")

# 绘制带图例的PDF
p <- ggplot(data = seurat_obj_cell@meta.data, aes(x = y, y = x, color = Tumor_boundary)) + 
  geom_point(size = 0.01) +
  scale_y_reverse() + coord_equal() + 
  theme_bw() + 
  theme(panel.grid = element_blank(), 
        axis.ticks = element_blank(),
        axis.text = element_blank(), 
        axis.title = element_blank(),
        panel.border = element_blank(),
        legend.title = element_blank(),
        legend.key.spacing.y = unit(0.8, "lines"), 
        legend.spacing.x = unit(1, "lines"), 
        legend.text = element_text(size = 20)) + 
  guides(color = guide_legend(override.aes = list(size = 7), byrow = TRUE)) + 
  scale_color_manual(values = color_boundary)

ggsave(file.path(out_dir, paste0(seurat_obj_cell@project.name, "_tumor_boundary_spatial_colored.pdf")), 
       p, width = 12, height = 10)