In stemangiola/tidyseurat: Brings Seurat to the Tidyverse

knitr::opts_chunk$set(echo = TRUE)

Case study code comparison

Calculate gamma-delta signature and plot

tidyseurat

seurat_obj <- readRDS("dev/PBMC_tidy_clean_scaled_UMAP_cluster_cell_type.rds")

seurat_obj = 
  seurat_obj %>%
  filter(first.labels == "T_cells") %>%
  RunPCA() %>% 
  RunUMAP(dims=1:30) %>%
  mutate(type=case_when(sample %in% c("GSE115189", "SRR11038995", "SRR7244582") ~ "A", TRUE ~ "B"))

library(tidygate)
library(ggplot2)
library(purrr)
library(patchwork)

# Calculate gamma delta signature
seurat_obj_sig = 

    seurat_obj %>%
    join_features(
      features = c("CD3D", "TRDC", "TRGC1", "TRGC2", "CD8A", "CD8B"), 
      shape = "wide", 
      assay = "SCT"
    ) %>%


    mutate(signature_score = 
             scales::rescale(CD3D + TRDC + TRGC1 + TRGC2, to=c(0,1)) - 
             scales::rescale(CD8A + CD8B, to=c(0,1))
    ) 

p1 = seurat_obj_sig %>%


  # Subsample
  add_count(sample, name = "tot_cells") %>%
  mutate(min_cells = min(tot_cells)) %>%
  group_by(sample) %>%
  sample_n(min_cells) %>%

  # Plot
  pivot_longer(cols=c("CD3D", "TRDC", "TRGC1", "TRGC2", "CD8A", "CD8B", "signature_score")) %>%
  mutate(value = case_when(value>0 ~ value)) %>%
  group_by(name) %>%
  mutate(value = scale(value)) %>%
  ggplot(aes(UMAP_1, UMAP_2, color=value)) +
  geom_point(shape=".") +
  facet_grid(type~name) +
  scale_color_viridis_c() +
  custom_theme

# Test differential abundance
p2 = seurat_obj_sig %>%

  # Gating
  mutate(gamma_delta = gate_chr(
        UMAP_1,
        UMAP_2, 
        .color =  signature_score, 
        .size=0.1
    )) %>%

  # Calculate proportions
  add_count(sample, name = "tot_cells") %>%
  count(sample, type, tot_cells, gamma_delta) %>%
  mutate(frac = n/tot_cells) %>%
  filter(gamma_delta == 1) %>%

  # Plot
  ggplot(aes(type, frac)) + 
  geom_boxplot() + 
  geom_point() +
  custom_theme

p = p1 / (p2  | plot_spacer()) +  plot_layout(guides = "collect")

ggsave("dev/summary_statistics.pdf", p,  device = "pdf", width = 183, height = 150, units = "mm", useDingbats=FALSE)

Seurat

library(Seurat)
library(gatepoints)
library(dplyr)

# Calculate gamma delta signature
signature_score_1 = 
  seurat_obj[c("CD3D", "TRDC", "TRGC1", "TRGC2"),] %>%
  GetAssayData(assay="SCT", slot="data") %>%
  colSums() %>%
  scales::rescale(to=c(0,1))

signature_score_2 = 
  seurat_obj[c("CD8A", "CD8B"),] %>%
  GetAssayData(assay="SCT", slot="data") %>%
  colSums() %>%
  scales::rescale(to=c(0,1))
seurat_obj$signature_score = signature_score_1 - signature_score_2

# Subsample
splits = colnames(seurat_obj) %>% split(seurat_obj$sample) 
min_size = splits %>% sapply(length) %>%  min()
cell_subset = splits %>%  lapply(function(x) sample(x, min_size)) %>%  unlist()
seurat_obj = seurat_obj[,cell_subset]

# Plot
DefaultAssay(seurat_obj) = "SCT"

seurat_obj %>%
FeaturePlot(
  features = c("signature_score", "CD3D", "TRDC", "TRGC1", "TRGC2", "CD8A", "CD8B"),
  split.by = "type",
  min.cutoff = 0.1
) 


# Gating
p = FeaturePlot(seurat_obj, features = "signature_score")
seurat_obj$within_gate = colnames(seurat_obj) %in% CellSelector(plot = p)

# Calculate proportions
seurat_obj[[]] %>%
  add_count(sample, name = "tot_cells") %>%
  count(sample, type, tot_cells, within_gate) %>%
  mutate(frac = n/tot_cells) %>%
  filter(within_gate == T) %>%

  # Plot
  ggplot(aes(type, frac)) + 
  geom_boxplot() + 
  geom_point()