脚本更新---高精度空间转录组的空间度量分析

library(tidyverse)
library(phenoptr)

#####示例数据在https://github.com/akoyabio/phenoptr/blob/main/data/sample_cell_seg_data.rda
#####大家可以研究一下，软件自带的

csd <- sample_cell_seg_data

csd <- csd %>% filter(Phenotype!='other')

distances <- find_nearest_distance(csd)
glimpse(distances)

csd_with_distance <- bind_cols(csd, distances)

merged_with_distance <- merged %>%
  group_by(`Sample Name`) %>%
  do(bind_cols(., find_nearest_distance(.)))

####Analyzing nearest neighbor distances
csd_with_distance %>% group_by(Phenotype) %>% 
  select(Phenotype, starts_with('Distance to')) %>% 
  summarize_all(~round(mean(.), 1))

ggplot(csd_with_distance, aes(`Distance to CD8+`, color=Phenotype)) +
  geom_density(size=1)

# Filter to just CD8+ and CK+ cells
cd8_cells = csd_with_distance %>% filter(select_rows(csd_with_distance, 'CD8+'))
ck_cells = csd_with_distance %>% filter(select_rows(csd_with_distance, 'CK+'))

# For each CD8+ cell, join with the data for the nearest CK+ cell
cd8_to_ck = cd8_cells %>% left_join(ck_cells, by=c('Cell ID CK+'='Cell ID'),
                                    suffix=c('', '.CK'))

####Show the CD8+ and CK+ cells with the nearest neighbors connected
# Read a background image and make a base plot
background_path = 
  system.file("extdata/sample/Set4_1-6plex_[16142,55840]_composite_image.jpg", package='phenoptr')
background = jpeg::readJPEG(background_path) %>% as.raster()
xlim = c(0, 934)
ylim = c(0, 700)
base_plot = ggplot(mapping=aes(`Cell X Position`, `Cell Y Position`)) %>% 
  phenoptr:::add_scales_and_background(background, xlim, ylim, scale_color='white') +
  labs(x='Cell X Position', y='Cell Y Position') +
  scale_color_manual('Phenotype', values=c('CD8+'='red', 'CK+'='cyan2'))

# Add lines and points
base_plot + geom_segment(data=cd8_to_ck,
               aes(xend=`Cell X Position.CK`, yend=`Cell Y Position.CK`),
               color='white') +
  geom_point(data=ck_cells, aes(color='CK+'), size=1) +
  geom_point(data=cd8_cells, aes(color='CD8+'), size=1) +
  labs(title='Nearest CK+ to each CD8+')

ck_to_cd8 = ck_cells %>% left_join(cd8_cells, by=c('Cell ID CD8+'='Cell ID'),
                                    suffix=c('', '.CD8'))
base_plot + 
  geom_segment(data=ck_to_cd8,
               aes(xend=`Cell X Position.CD8`, yend=`Cell Y Position.CD8`),
               color='white') +
  geom_point(data=ck_cells, aes(color='CK+'), size=1) +
  geom_point(data=cd8_cells, aes(color='CD8+'), size=1) +
  labs(title='Nearest CD8+ to each CK+')

mutual = ck_to_cd8 %>% filter(`Cell ID`==`Cell ID CK+.CD8`)
base_plot +
  geom_segment(data=mutual,
               aes(xend=`Cell X Position.CD8`, yend=`Cell Y Position.CD8`),
               size=1, color='white') +
  geom_point(data=ck_cells, aes(color='CK+'), size=1) +
  geom_point(data=cd8_cells, aes(color='CD8+'), size=1) +
  labs(title='Mutual nearest neighbors - CD8+ and CK+')

count_within(csd, from='CD68+', to='CK+', radius=25)

count_within(csd, from='CK+', to='CD68+', radius=25)