hello，大家好，今天我们来学习一个很实用的分析方法，scDiffCom，主要研究组间通讯差异，比如正常和疾病样本之间，不同发育样本之间的通讯差异等等，比较实用。大家多多学习，文章在scAgeCom: a murine atlas of age-related changes in intercellular communication inferred with the package scDiffCom。
我们来看看示例分析，帮助我们理解。
scDiffCom 的目标是研究在 scRNA-seq 数据中检测到的细胞间相互作用如何在两种生物条件之间发生变化。
Some important terminology includes:
- LRI: a known ligand-receptor interaction of the form Apoe:Ldlr (simple) or Col3a1:Itgb1_Itga2 (complex/heteromeric)
- CCI (cell-cell interaction): a communication signal of the form (B cell, T cell; Apoe:Ldlr) where B cell is the emitter cell type expressing the ligand Apoe and T cell is the receiver cell type expressing the receptor Ldlr.
Standard differential analysis
scDiffCom 使用 scRNA-seq Seurat 对象作为输入。 在这里我们使用已下采样的模型对象（R包附带有）。 它包含分为 5 种细胞类型以及将对其进行差异分析的两个年龄组（“年轻”和“老年”）的细胞。
library(Seurat)
## Attaching SeuratObject
library(scDiffCom)
library(data.table)
# Load the toy-model object
seurat_object <- scDiffCom::seurat_sample_tms_liver
seurat_object
## An object of class Seurat
## 726 features across 468 samples within 1 assay
## Active assay: RNA (726 features, 0 variable features)
# Cell-type annotations must be present beforehand in the object meta.data
head(seurat_object[["cell_type"]])
## cell_type
## A16_B000767_B009945_S16.mm10-plus-4-0 T cell
## A16_D045313_B009942_S16.mm10-plus-4-0 endothelial cell of hepatic sinusoid
## A21_B001247_B009941_S21.mm10-plus-4-0 B cell
## A22_B001247_B009941_S22.mm10-plus-4-0 B cell
## A2_D045313_B009942_S2.mm10-plus-4-0 hepatocyte
## A5_B001247_B009941_S5.mm10-plus-4-0 T cell
# The two conditions of interest must also be specified in the object meta.data
table(seurat_sample_tms_liver[["age_group"]])
##
## OLD YOUNG
## 250 218
(可选)研究配体-受体相互作用
scDiffcom 从已知参与配体-受体相互作用 (LRI) 的基因的表达推断细胞类型到细胞类型的通信模式。 软件包包含它自己的内部精选 LRI（用于人和小鼠）数据库，这些数据库是从以前的研究中检索到的。
执行分析时无需显式调用 LRI 数据库。 但可以按如下方式访问和探索它们：
# Load, e.g., the mouse database
data(LRI_mouse)
# Display the data.table of LRIs (more information available in other columns)
LRI_mouse$LRI_curated[, c("LRI")]
## LRI
## 1: 1700013F07Rik:Plscr4
## 2: 2300002M23Rik:Ddr1
## 3: 9530003J23Rik:Itgal
## 4: A2m:Lrp1
## 5: Aanat:Mtnr1a
## ---
## 4737: a:Mc1r
## 4738: a:Mc2r
## 4739: a:Mc3r
## 4740: a:Mc4r
## 4741: a:Mc5r
# Display the data.table of GO Terms attached to each LRI (more information available in other columns)
LRI_mouse$LRI_curated_GO[, c("LRI", "GO_NAME")]
## LRI GO_NAME
## 1: 1700013F07Rik:Plscr4 cellular_component
## 2: 1700013F07Rik:Plscr4 cell part
## 3: 1700013F07Rik:Plscr4 molecular_function
## 4: 1700013F07Rik:Plscr4 binding
## 5: 1700013F07Rik:Plscr4 biological_process
## ---
## 267431: a:Mc5r biological regulation
## 267432: a:Mc5r regulation of biological process
## 267433: a:Mc5r regulation of cellular process
## 267434: a:Mc5r signal transduction
## 267435: a:Mc5r regulation of metabolic process
Run default differential analys
所有结果都是通过简单地调用 run_interaction_analysis 获得的。 检测和差异分析依赖于置换检验。 默认情况下，执行 1000 次迭代，因为这足以进行快速的解释性分析。 In general, we would recommend using 10’000 permutations. Parallel computing can be easily enabled by loading the future package and setting the plan accordingly.
# Load the future package (optional)
library(future)
plan(sequential) # sequentially in the current R process, equivalent to do nothing
#plan(multisession, workers = 4) # background R sessions
#plan(multicore, workers = 4) # forked R processes, not Windows/not RStudio
# Run differential analysis with default parameters
scdiffcom_object <- run_interaction_analysis(
seurat_object = seurat_object,
LRI_species = "mouse",
seurat_celltype_id = "cell_type",
seurat_condition_id = list(
column_name = "age_group",
cond1_name = "YOUNG",
cond2_name = "OLD"
)
)
## Extracting data from assay 'RNA' and slot 'data' (assuming normalized log1p-transformed data).
## Converting normalized data from log1p-transformed to non-log1p-transformed.
## Input data: 726 genes, 468 cells and 5 cell-types.
## Input ligand-receptor database: 4741 mouse interactions.
## Number of LRIs that match to genes present in the dataset: 1219.
## Type of analysis to be performed: differential analysis between YOUNG and OLD cells.
## Total number of potential cell-cell interactions (CCIs): 30475 (5 * 5 * 1219).
## Performing permutation analysis (1000 iterations by batches of 1000) on 9366 potential CCIs.
## Performing batch 1 of 1.
## Filtering and cleaning 'raw' CCIs.
## Returning 1986 detected CCIs.
## Performing over-representation analysis on the categories: LRI, LIGAND_COMPLEX, RECEPTOR_COMPLEX, ER_CELLTYPES, EMITTER_CELLTYPE, RECEIVER_CELLTYPE, GO_TERMS, KEGG_PWS.
## Successfully returning final scDiffCom object.
The output of run_interaction_analysis is an S4 object of class scDiffCom:
scdiffcom_object
## An object of class scDiffCom with name scDiffCom_object
## Analysis performed: differential analysis between YOUNG and OLD cells
## 1986 detected CCIs across 5 cell types
## Over-representation results for LRI, LIGAND_COMPLEX, RECEPTOR_COMPLEX, ER_CELLTYPES, EMITTER_CELLTYPE, RECEIVER_CELLTYPE, GO_TERMS, KEGG_PWS
Explore the detected CCIs
The slot cci_table_detected 是一个 data.table，它包含所有与生物学相关的 CCI（每行一个），其中的列提供有用的信息，例如它们如何在两个条件（这里是年轻和年老）之间进行调节：
# Retrieve and display all detected CCIs
CCI_detected <- GetTableCCI(scdiffcom_object, "detected", simplified = TRUE)
CCI_detected[, c("CCI", "REGULATION")]
## CCI REGULATION
## 1: B cell_B cell_Adam10:Fcer2a FLAT
## 2: B cell_B cell_Btla:Cd79a NSC
## 3: B cell_B cell_Cd48:Cd2 NSC
## 4: B cell_B cell_Cd72:Sema4d NSC
## 5: B cell_B cell_Copa:Cd74 FLAT
## ---
## 1982: myeloid leukocyte_myeloid leukocyte_Vasp:Itga4 FLAT
## 1983: myeloid leukocyte_myeloid leukocyte_Vcam1:Itgb1_Itga4 FLAT
## 1984: myeloid leukocyte_myeloid leukocyte_Vcam1:Itgb1_Itga9 FLAT
## 1985: myeloid leukocyte_myeloid leukocyte_Vcam1:Itgb7_Itga4 FLAT
## 1986: myeloid leukocyte_myeloid leukocyte_Vcam1:Msn FLAT
# Number of CCIs per regulation type
table(CCI_detected$REGULATION)
##
## DOWN FLAT NSC UP
## 120 973 858 35
检测到的 CCI 可以通过查看 data.table 来一个一个地探索，更广泛地基于的over-representation analysis或像这个火山图这样的图表：
if (!require("ggplot2")) install.packages("ggplot2")
## Loading required package: ggplot2
library(ggplot2)
ggplot(
CCI_detected,
aes(
x = LOGFC,
y = -log10(BH_P_VALUE_DE + 1E-2),
colour = REGULATION
)
) + geom_point(
) + scale_colour_manual(
values = c("UP" = "red", "DOWN" = "blue", "FLAT" = "green", "NSC" = "grey")
) + xlab(
"log(FC)"
) + ylab(
"-log10(Adj. p-value)"
)
图片.png
Explore ORA results
The slot ora_table contains a list of data.tables giving the results of the over-representation analysis performed on predefined categories:
# Retrieve the ORA results
ORA_results <- GetTableORA(scdiffcom_object, categories = "all", simplified = TRUE)
# Categories available
names(ORA_results)
## [1] "LRI" "LIGAND_COMPLEX" "RECEPTOR_COMPLEX"
## [4] "ER_CELLTYPES" "EMITTER_CELLTYPE" "RECEIVER_CELLTYPE"
## [7] "GO_TERMS" "KEGG_PWS"
# Explore, e.g., the top 10 LRIs up-regulated based on their ORA score
ORA_results$LRI[, c("VALUE", "ORA_SCORE_UP", "OR_UP", "BH_P_VALUE_UP")][order(-ORA_SCORE_UP)][1:10]
## VALUE ORA_SCORE_UP OR_UP BH_P_VALUE_UP
## 1: H2-D1:Cd8b1 Inf Inf 7.292044e-07
## 2: H2-K1:Cd8b1 Inf Inf 7.292044e-07
## 3: H2-Q6:Cd8b1 Inf Inf 3.012525e-05
## 4: Lck:Cd8a_Cd8b1 Inf Inf 2.883823e-02
## 5: H2-T23:Cd8b1 30.28268 178.35580 8.928196e-05
## 6: Lgals1:Lgals1 20.42120 59.95999 3.485409e-04
## 7: H2-D1:Cd8a 14.36497 58.14560 3.542486e-03
## 8: H2-K1:Cd8a 14.36497 58.14560 3.542486e-03
## 9: H2-T22:Cd8b1 14.36497 58.14560 3.542486e-03
## 10: Ltb:Tnfrsf1a 14.36497 58.14560 3.542486e-03
# Plot the most over-represented up-regulated LRIs
# note that PlotORA returns a ggplot object that you can further optimize (e.g. here to place the legend)
PlotORA(
object = scdiffcom_object,
category = "LRI",
regulation = "UP"
) + theme(
legend.position = c(0.85, 0.4),
legend.key.size = unit(0.4, "cm")
图片.png
A summary of the changes in intercellular communication is given by the following graph that displays the over-represented cell types and cell-type pairs:
if (!require("visNetwork")) install.packages("visNetwork")
## Loading required package: visNetwork
if (!require("igraph")) install.packages("igraph")
## Loading required package: igraph
##
## Attaching package: 'igraph'
## The following objects are masked from 'package:future':
##
## %->%, %<-%
## The following objects are masked from 'package:stats':
##
## decompose, spectrum
## The following object is masked from 'package:base':
##
## union
if (!require("kableExtra")) install.packages("kableExtra")
## Loading required package: kableExtra
if (!require("RColorBrewer")) install.packages("RColorBrewer")
## Loading required package: RColorBrewer
library(visNetwork)
library(igraph)
library(kableExtra)
library(RColorBrewer)
BuildNetwork(
object = scdiffcom_object
)
图片.png
Additional functionalities
Modifying default parameters
# display the first three parameters
head(GetParameters(scdiffcom_object), 3)
## $object_name
## [1] "scDiffCom_object"
##
## $LRI_species
## [1] "mouse"
##
## $seurat_celltype_id
## [1] "cell_type"
比较实用的一个工具，没什么新意，但是比较方便
生活很好，等你超越
