开工了，到了我们上班的时间了，刚开始，我们来分享一些新出现的方法，看看方法上有哪些更新，我们虽然做不了潮流的引领者，但是要跟随潮流，保证自己不被时代抛弃，今天分享的文献在Comparative analysis of cell-cell communication at single-cell resolution,里面对当前细胞通讯的方法做了一定的总结，我们借鉴一下，并看看作者的方法有哪些可取之处。

目前的单细胞通讯的分析特点

• at the level of the cell type or cluster, discarding single-cell level information.（这个大家应该都非常清楚了，通讯强度就是基因表达的平均值相乘）

然而真实的细胞通讯

• 细胞通讯 does not operate at the level of the group; rather, such interactions take place between individual cells（真正的细胞通讯是细胞之间，而不是细胞类型群体之间，这也是空间转录组为什么分析空间临近通讯的原因）。

如果我们要真正的研究细胞之间的通讯，那么

• analyze interactions at the level of the single cell
• leverage the full information content contained within scRNA-seq data by looking at up- and down-stream cellular activity
• enable comparative analysis between conditions;
• robust to multiple experimental designs

不知道大家对单细胞做细胞通讯的方式有什么更好的建议呢？？？，我们来看看作者开发的方法，Scriabin，an adaptable and computationally-efficient method for CCC analysis.Scriabin 通过结合收集的配体-受体相互作用数据库、下游细胞内信号传导模型、基于anchors的数据集集成和基因网络分析，以单细胞分辨率剖析复杂的交流途径，以在单细胞分辨率下恢复具有生物学意义的 CCC edges。

Scriabin的分析框架

• the cell-cell interaction matrix workflow, optimal for smaller datasets, analyzes communication methods used for each cell-cell pair in the dataset（细胞对分析通讯），CCC 的基本单位是表达配体的发送细胞 N_i，这些配体被接收细胞 N_j 表达的同源受体接收。 Scriabin 通过计算数据集中每对细胞对每个配体-受体对表达的几何平均值，在细胞-细胞相互作用矩阵 M 中编码此信息.由于配体-受体相互作用是定向的，Scriabin将每个细胞分别视为“发送者”（配体表达）和“接收者”（受体表达），从而保持 CCC 网络的定向性质。 M 可以类似于基因表达矩阵进行处理，并用于降维、聚类和差异分析
• the summarized interaction graph workflow, designed for large comparative analyses, identifies cell-cell pairs with different total communicative potential between samples(识别有意义的配受体对)
• the interaction program discovery workflow, suitable for any dataset size, finds modules of co-expressed ligand-receptor pairs（第三点尤为重要）。

图片.png

图片.png

具体的方法

Generation of cell-cell interaction matrix

将一对细胞之间的细胞-细胞相互作用向量定义为每个同源配体-受体对的表达值的几何平均值。 形式上，sender cell N_i 和receiver cell N_j 之间的交互向量 V 由下式给出:

图片.png

矩阵的下游分析（就是Seurat的基本操作）

图片.png

通讯programs的分析，Interaction program discovery workflow

Iterative approximation of a ligand-receptor pair topological overlap matrix (TOM)

图片.png

Identification and significance testing of interaction programs（层次聚类,WGCNA）

图片.png

看一看示例代码

安装

devtools::install_github("BlishLab/scriabin", ref = "main")

single-dataset(单一数据样本)

library(Seurat)
library(SeuratData)
library(scriabin)
library(tidyverse)
#####示例数据
install.packages("https://seurat.nygenome.org/src/contrib/panc8.SeuratData_3.0.2.tar.gz", repos = NULL, type = "source") 
library(panc8.SeuratData)
panc8 <- LoadData("panc8")
####Analyze the indrop dataset for the interaction programs tutorial
panc_fl <- subset(panc8, cells = colnames(panc8)[panc8$tech=="fluidigmc1"])
panc_fl <- SCTransform(panc_fl, verbose = F) %>%
  RunPCA(verbose = F) %>%
  RunUMAP(dims = 1:30, verbose = F)
DimPlot(panc_fl, group.by = "celltype", label = T, repel = T) + NoLegend()

图片.png

First we generate a cell-cell interaction matrix object on the alpha and beta cells in the dataset. Then we scale, find variable features, run PCA and UMAP, find neighbor graphs and clusters, treating this matrix just like we would a gene expression matrix

panc_fl_ccim <- GenerateCCIM(subset(panc_fl, cells = colnames(panc_fl)[panc_fl$celltype %in% c("alpha","beta")]))
panc_fl_ccim <- ScaleData(panc_fl_ccim) %>% 
  FindVariableFeatures() %>%
  RunPCA() %>% 
  RunUMAP(dims = 1:10) %>%
  FindNeighbors(dims = 1:10) %>%
  FindClusters(resolution = 0.2)
p1 <- DimPlot(panc_fl_ccim, group.by = "sender_celltype")
p2 <- DimPlot(panc_fl_ccim, group.by = "receiver_celltype")
p1+p2
DimPlot(panc_fl_ccim, label = T, repel = T) + NoLegend()

cluster4_5.edges <- FindMarkers(panc_fl_ccim, ident.1 = "4", ident.2 = "5")
cluster4_5.edges %>% top_n(20,wt = abs(avg_log2FC))

CCIMFeaturePlot(panc_fl_ccim, seu = panc_fl, 
            features = c("EFNA1","NPNT"), type_plot = "sender")

图片.png

multi-dataset

library(Seurat)
library(SeuratData)
library(scriabin)

scriabin::load_nichenet_database()

library(SeuratData)
InstallData("ifnb")
ifnb <- UpdateSeuratObject(LoadData("ifnb"))
ifnb <- PercentageFeatureSet(ifnb, pattern = "MT-", col.name = "percent.mt") %>%
  SCTransform(ifnb, vars.to.regress = "percent.mt", verbose = F) %>%
  RunPCA(ifnb, verbose = F) %>% 
  FindNeighbors(ifnb, dims = 1:30, verbose = F) %>%
  FindClusters(ifnb, verbose = F) %>%
  RunUMAP(ifnb, dims = 1:30, verbose = F)
DimPlot(ifnb, label = T, repel = T) + NoLegend()
DimPlot(ifnb, label = T, repel = T, group.by = "seurat_annotations") + NoLegend()
DimPlot(ifnb, group.by = "stim")

图片.png

This dataset is composed of two sub-datasets: PBMCs stimulated with IFNB (STIM) or control (CTRL).

If we assume there is no batch effect between these two datasets, we can observe from these dimensionality reduction projections the intense transcriptional perturbation caused by this stimulation. This is one situation in which it is easy to see why clustering/subclustering for high resolution CCC analysis can be a problematic task: the degree of perturbation is so high that the monocytes from these two datasets will never cluster together. We need a way to align these subspaces together. Given a monocyte in CTRL, what is its molecular counterpart in STIM?

We take recent progress in dataset integration methodology to develop a high-resolution alignment process we call "binning", where we assign each cell a bin identity so that we maximize the similarity of the cells within a bin, and maximize the representation of all the samples we wish to compare in each bin. After bins are identified, they are tested for the significance of their connectivity.

Dataset binning

There are two workflow options for binning in terms of bin identification:

1. Bin by coarse cell identities (recommended). If the user specifies coarse cell identities (eg. in a PBMC dataset, T/NK cells, B cells, myeloid cells), this can prevent anchors from being identified between cells that the user believes to be in non-overlapping cell categories. This can result in cleaner bins, but there must be enough cells in each coarse identity from each dataset to proceed properly. A downside is that there may be some small groups of cells that aren't related to major cell populations (eg. a small population of platelets in a PBMC dataset).
2. Bin all cells. Without specifying coarse cell identities for binning, potentially spurious associations may form between cells (especially in lower quality samples).

Significance testing proceeds by a permutation test: comparing connectivity of the identified bin against randomly generated bins. There are two workflow options for binning in terms of generating random bins:

1. Pick cells from the same cell type in each random bin (recommended). If the user supplies granular cell type identities, a random bin will be constructed with cells from the same cell type identity. The more granular the cell type calls supplied, the more rigorous the significance testing.
2. Without supplying granular cell type identities, the dataset will be clustered and the cluster identities used for significance testing. Generating random bins across the entire dataset would result in very few non-significant bins identified.

Here we will define a coarse cell identity to use for bin identification, and then use the Seurat-provided annotations for significance testing.

Additionally, Scriabin allows different reductions to be used for neighbor graph construction. For example, the results from a reference-based sPCA can be used for binning if the cell type relationships in the reference are considered more informative.

ifnb$coarse_ident <- mapvalues(ifnb$seurat_annotations, from = unique(ifnb$seurat_annotations),
                               to = c("myeloid","myeloid","T/misc","T/misc",
                                      "T/misc","T/misc","T/misc","B",
                                      "B","myeloid","myeloid","T/misc","T/misc"))
ifnb <- BinDatasets(ifnb, split.by = "stim", dims = 1:30, 
                    coarse_cell_types = "coarse_ident", sigtest_cell_types = "seurat_annotations")

图片.png

ifnb_split <- SplitObject(ifnb, split.by = "stim")
sum_ig <- AssembleInteractionGraphs(ifnb, by = "prior", split.by = "stim")
ifnb_split <- pblapply(ifnb_split, function(x) {BuildPriorInteraction(x, correct.depth = T)})
ogig <- lapply(ifnb_split, function(x) {
  as.matrix(x@graphs$prior_interaction)
})

图片.png

interaction-programs

library(Seurat)
library(SeuratData)
library(scriabin)
library(tidyverse)
library(ComplexHeatmap)
library(cowplot)
install.packages("https://seurat.nygenome.org/src/contrib/panc8.SeuratData_3.0.2.tar.gz", repos = NULL, type = "source") 
library(panc8.SeuratData)
panc8 <- LoadData("panc8")

panc_id <- subset(panc8, cells = colnames(panc8)[panc8$tech=="indrop"])
panc_id <- SCTransform(panc_id, verbose = F) %>%
  RunPCA(verbose = F) %>%
  RunUMAP(dims = 1:30, verbose = F)
DimPlot(panc_id, group.by = "celltype", label = T, repel = T) + NoLegend()

图片.png

panc_ip <- InteractionPrograms(panc_id, iterate.threshold = 300)
#test for interaction program significance
panc_ip_sig <- InteractionProgramSignificance(panc_ip, n.replicate = 500)
#score cells by expression of interaction program
panc_id <- ScoreInteractionPrograms(panc_id, panc_ip_sig)

panc_id_ip_lig <- as.matrix(panc_id@meta.data %>% 
  select("celltype",
         starts_with("ligands")) %>%
  group_by(celltype) %>%
  summarise_if(is.numeric, mean) %>% column_to_rownames("celltype"))
Heatmap(panc_id_ip_lig, show_column_names = F)
panc_id_ip_rec <- as.matrix(panc_id@meta.data %>% 
  select("celltype",
         starts_with("receptors")) %>%
  group_by(celltype) %>%
  summarise_if(is.numeric, mean) %>% column_to_rownames("celltype"))
Heatmap(panc_id_ip_rec, show_column_names = F)

图片.png

act_stellate_ip <- panc_id_ip_lig["activated_stellate",]
poi <- gsub("ligands_","",names(which(act_stellate_ip==max(act_stellate_ip))))
#Seurat's FeaturePlot has a nice option to blend expression of two features together on the same plot
p <- FeaturePlot(panc_id, 
          features = c(paste0("ligands_",poi),
                            paste0("receptors_",poi)), 
          blend = T, combine = F, 
          cols = c("grey90","purple","yellowgreen"), order = T)
p[[3]] + NoLegend()
DimPlot(panc_id, group.by = "celltype", label = T, repel = T) + NoLegend()

图片.png

moi <- reshape2::melt(mod_df %>% dplyr::filter(name==poi) %>%
  select("lr_pair",contains("connectivity"))) %>% arrange(-value)
moi$lr_pair <- factor(moi$lr_pair, levels = unique(moi$lr_pair))
ggplot(moi, aes(x = lr_pair, y = value, color = variable)) + 
  geom_point() + theme_cowplot() + ggpubr::rotate_x_text() + labs(x = NULL, y = "Intramodular\nconnectivity")

图片.png

最后问一句，大家觉得这个方法改进之处合理么？

生活很好，有你更好