In arc85/celltalker: Analysis of cell-cell communication from single-cell RNAseq data

This vignette demonstrates an extension of the core functionality of celltalker to replicate samples. The idea is that this functionality can be used to identify consistently differentially expressed ligand and receptor interactions. To demonstrate this functionality, we will leverage our previously published dataset of healthy donor PBMC and healthy donor tonsils.

The data preprocessing steps are documented in the data-raw directory. We will start with a low-dimensional visualization with annotated cell types.

Load packages and setup analysis

suppressMessages({
library(celltalker)
library(Seurat)
library(dplyr)
library(magrittr)
})

## Load data provided with package
data("filtered_lig_rec")
data("overall_metadata")
data("overall_umap")

## Create a Seurat object based on 5 healthy donor PBMC and 5 healhy donor tonsils
ser_obj <- CreateSeuratObject(filtered_lig_rec,meta.data=overall_metadata)
ser_obj[["umap"]] <- CreateDimReducObject(embeddings = overall_umap, key = "UMAP_", assay = DefaultAssay(ser_obj))

DimPlot(ser_obj,group.by="cell_types",split.by="sample_type")

Run celltalker

Next, we will split the data into PBMC and Tonsil Seurat objects and will identify the top ligand and receptor interactions in each of these datasets

## Split dataset
ser_split <- SplitObject(ser_obj,split="sample_type")

## Check out the split data
ser_split

## Run celltalker - PBMC
pbmc_interactions <- celltalk(input_object=ser_split[["PBMC"]],
  metadata_grouping="cell_types",
  ligand_receptor_pairs=ramilowski_pairs,
  number_cells_required=50,
  min_expression=50,
  max_expression=20000,
  scramble_times=10)

## Check out the interactions - PBMC
pbmc_interactions %>%
  mutate(p_val_adj=p.adjust(p_val,method="fdr")) %>%
  filter(p_val_adj<0.05)

## Run celltalker - Tonsil
tonsil_interactions <- celltalk(input_object=ser_split[["Tonsil"]],
  metadata_grouping="cell_types",
  ligand_receptor_pairs=ramilowski_pairs,
  number_cells_required=50,
  min_expression=50,
  max_expression=20000,
  scramble_times=10)

## Check out the interactions - tonsil
tonsil_interactions %>%
  mutate(p_val_adj=p.adjust(p_val,method="fdr")) %>%
  filter(p_val_adj<0.05)

Create circos plots

## Plot top 3 interactions - PBMC
# Identify the top 3 interactions for cell_type1
top_stats_pbmc <- pbmc_interactions %>%
  mutate(p_val_adj=p.adjust(p_val,method="fdr")) %>%
  filter(p_val_adj<0.05) %>%
  group_by(cell_type1) %>%
  top_n(10,interact_ratio) %>%
  ungroup()

# Assign colors to cell types
all_cell_types <- unique(ser_obj[["cell_types"]][,1])

# Define colors
colors_use <- colorRampPalette(RColorBrewer::brewer.pal(n=8,"Set2"))(length(all_cell_types))
names(colors_use) <- all_cell_types

# Suppress messages to silence the circlize functions
suppressMessages(
circos_plot(ligand_receptor_frame=top_stats_pbmc,
  cell_group_colors=colors_use,
  ligand_color="blue",
  receptor_color="red",
  cex_outer=0.5,
  cex_inner=0.4)
)

## Plot top 3 interactions - Tonsil
# Identify the top 3 interactions for cell_type1
top_stats_tonsil <- tonsil_interactions %>%
  mutate(p_val_adj=p.adjust(p_val,method="fdr")) %>%
  filter(p_val_adj<0.05) %>%
  group_by(cell_type1) %>%
  top_n(10,interact_ratio) %>%
  ungroup()

# Suppress messages to silence the circlize functions
suppressMessages(
circos_plot(ligand_receptor_frame=top_stats_tonsil,
  cell_group_colors=colors_use,
  ligand_color="blue",
  receptor_color="red",
  cex_outer=0.5,
  cex_inner=0.4)
)

Identify significantly different interactions across patient samples

## Comparison of group interactions
# Use top 10 interactions from tonsils as input
group_stats <- compare_group_interactions(seurat_object=ser_obj,
  interaction_stats=top_stats_tonsil,
  sample_replicates="sample_id",
  sample_groups="sample_type",
  metadata_grouping="cell_types")

# Extract p values into data.frame and add FDR
mod_p_vals <- do.call(rbind,
  lapply(group_stats,function(x) {
    if (class(x) != "lm") stop("Not an object of class 'lm' ")
    f <- summary(x)$fstatistic
    p <- pf(f[1],f[2],f[3],lower.tail=F)
    attributes(p) <- NULL
    return(p)
  })
)

mod_p_vals <- data.frame(mod_p_vals,p_val_adj=p.adjust(mod_p_vals[,1],method="fdr"))

mod_p_vals

# Boxplot of joint weight interaction across replicates
boxplot_group_interaction(seurat_object=ser_obj,
  interaction_stats=top_stats_tonsil,
  sample_replicates="sample_id",
  sample_groups="sample_type",
  metadata_grouping="cell_types",
  ligand="CD40LG",
  receptor="CD40",
  cell_type1="CD4 T cells",
  cell_type2="B cells")

arc85/celltalker documentation built on July 2, 2023, 2:07 p.m.