本期内容已同步分享至单细胞天地
细胞通讯是指细胞与细胞之间的联系。细胞和人类一样，多细胞生物的很多细胞会相互作用，形成“细胞社会”，在这个社会里，细胞与细胞之间会发生相互作用和信息的传递，细胞建立通讯联络是必需的。如生物体的生长发育、分化、各种组织器官的形成、组织的维持以及它们各种生理活动的协调。经典的例子莫过于 神经细胞之间的神经递质的传递与接收。
细胞与细胞之间的通讯有三种形式：
- 细胞之间的直接接触；
- 细胞之间通过其间的胞外基质相互联系；
- 细胞之间通过分子间的相互作用产生联系。
在这个过程中，“受体-配体”的概念十分重要，什么是受体，什么是配体？受体与配体之间结合的结果是受体被激活，并产生受体激活后续信号传递的基本步骤。在单细胞分析当中，不同细胞类型之间的通讯可能会对某些生物学过程具有重要意义，因此利用单细胞数据进行细胞通讯分析是单细胞高级分析的一大重点。
CellChat
CellChat是一款2021年发表于Nature Communications的单细胞细胞通讯分析工具。
CellChat上游分析
1. 安装 CellChat
devtools::install_github("sqjin/CellChat")
2. 细胞通讯分析
这里我们使用 pbmc3k 的公共数据集为例来一起学习如何利用 CellChat 进行细胞通讯分析。
2.1 数据输入
CellChat 需要的输入文件包括：
-
细胞的基因表达矩阵（已经经过normalize的）
不同的基因作为行名，细胞ID作为列名。这一点和Seurat对象的结构是保持一致的。 -
细胞的metadata信息
细胞层面的信息，这里可以是我们上游分析的注释信息，可以直接把SeuratObject的metadata提取出来使用。
所以我们的输入文件可以这样准备：
library(CellChat)
library(Seurat)
library(SeuratData)
data("pbmc3k.final")
data <- GetAssayData(object = pbmc3k.final, slot = 'data')
meta <- pbmc3k.final@meta.data
这里涉及到一个小知识点：SeuratObject中RNA assay中不同slot所存储的是什么值：
-
counts：原始的基因counts数，也就是简单reads计数的结果，对于10X Genomics平台数据来说，这是cellranger运行的结果； -
data：这是原始表达矩阵经过质控之后进行NormalizeData()的数据，NormalizeData去除了不同细胞之间测序深度的差异，同时对结果进行了对数化，这个数据是CellChat想要的数据。 -
scale.data：这是函数ScaleData()运行的结果，主要是将每个基因的表达量转换成了符合标准正态分布的数据，从而降低部分细胞异常表达值的影响。
因为实际上我们在使用CellChat时是都已经默认上游的处理已经完成了，所以我在这里不打算介绍CellChat本身自带的函数去normalize原始counts的分析。
2.2 创建CellChat对象
CellChat对象和Seurat对象很像，我们可以这样来进行创建：
cellchat <- createCellChat(object = data,
meta = meta,
group.by = 'seurat_annotations')
这里的group.by参数是不可少的，来自于我们的metadata的列名，我们一般指定为细胞类型注释的结果，这是有利于我们的分析结果的。
除此之外，不得不介绍的还有这个函数不仅仅支持将你自己提取的表达矩阵作为输入，还支持直接用SeuratObject作为输入，不过需要注意的是，如果是多组样本整合的SeuratObject，我们不能利用integrated assay作为输入，因为其含有负值，所以我们可用下面的这段代码实现和前面代码一样的目的：
cellchat <- createCellChat(object = pbmc3k.final,
meta = meta,
group.by = 'seurat_annotations',
assay = 'RNA')
此外还有一个参数do.sparse不要去改动它（默认为TRUE），用稀疏矩阵处理单细胞数据能够节省更多的空间和时间。
简单介绍一下CellChat数据的结构：
str(cellchat)
Formal class 'CellChat' [package "CellChat"] with 14 slots
..@ data.raw : num[0 , 0 ]
..@ data :Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
.. .. ..@ i : int [1:2238732] 29 73 80 148 163 184 186 227 229 230 ...
.. .. ..@ p : int [1:2639] 0 779 2131 3260 4220 4741 5522 6304 7094 7626 ...
.. .. ..@ Dim : int [1:2] 13714 2638
.. .. ..@ Dimnames:List of 2
.. .. .. ..$ : chr [1:13714] "AL627309.1" "AP006222.2" "RP11-206L10.2" "RP11-206L10.9" ...
.. .. .. ..$ : chr [1:2638] "AAACATACAACCAC" "AAACATTGAGCTAC" "AAACATTGATCAGC" "AAACCGTGCTTCCG" ...
.. .. ..@ x : num [1:2238732] 1.64 1.64 2.23 1.64 1.64 ...
.. .. ..@ factors : list()
..@ data.signaling: num[0 , 0 ]
..@ data.scale : num[0 , 0 ]
..@ data.project : num[0 , 0 ]
..@ net : list()
..@ netP : list()
..@ meta :'data.frame': 2638 obs. of 7 variables:
.. ..$ orig.ident : Factor w/ 1 level "pbmc3k": 1 1 1 1 1 1 1 1 1 1 ...
.. ..$ nCount_RNA : num [1:2638] 2419 4903 3147 2639 980 ...
.. ..$ nFeature_RNA : int [1:2638] 779 1352 1129 960 521 781 782 790 532 550 ...
.. ..$ seurat_annotations: Factor w/ 9 levels "Naive CD4 T",..: 2 4 2 3 7 2 5 5 1 6 ...
.. ..$ percent.mt : num [1:2638] 3.02 3.79 0.89 1.74 1.22 ...
.. ..$ RNA_snn_res.0.5 : Factor w/ 9 levels "0","1","2","3",..: 2 4 2 3 7 2 5 5 1 6 ...
.. ..$ seurat_clusters : Factor w/ 9 levels "0","1","2","3",..: 2 4 2 3 7 2 5 5 1 6 ...
..@ idents : Factor w/ 9 levels "Naive CD4 T",..: 2 4 2 3 7 2 5 5 1 6 ...
..@ DB : list()
..@ LR : list()
..@ var.features : list()
..@ dr : list()
..@ options :List of 1
.. ..$ mode: chr "single"
截止到目前为止，我们的CellChat对象结构如上所示，可以看到：
- 我们输入的数据存在了
data这个稀疏矩阵里面，而后面的data.signaling等对象还有待后面的分析进行填充； - 此外，我们应该重点关注的就是
meta部分，这个部分是一个数据框，所以我们还是可以在创建CellChat对象之后来更改细胞的信息，这是比较方便的，可以看到目前CellChat对象对细胞的分析是基于seurat_annotations的，这是我们前面设置好的，如果我们想要改变只需要setIdent(cellchat, ident.use = "labels")就可以了，同时我们还可以通过levels(cellchat@idents)来查看当前的细胞分组信息。
2.3 选择受体配体数据库
CellChat 有一个专门的数据库，叫做CellChatDB，这个数据库是 CellChat 的作者们通过阅读大量文献，手动整理出来的“受体-配体”对，目前有人、鼠以及斑马鱼的版本。其中人的叫做 CellChatDB.human，鼠的叫做 CellChatDB.mouse，斑马鱼的叫做 CellChatDB.zebrafish。关于这两个数据库中具体的“受体-配体”对信息可以通过showDatabaseCategory()函数获得，在这里就不赘述了。
-
mouse
2,021validated molecular interactions, including60% of secrete autocrine/paracrine signaling interactions,21% of extracellular matrix (ECM)-receptor interactionsand19% of cell-cell contact interactions. -
human
1,939validated molecular interactions, including61.8% of paracrine/autocrine signaling interactions,21.7% of extracellular matrix (ECM)-receptor interactionsand16.5% of cell-cell contact interactions.
加载数据库：
CellChatDB <- CellChatDB.human
来看一下数据库的结构（option）：
head(CellChatDB$interaction)
最后，千万不要忘了把数据库嵌入到 CellChat 对象中，否则后面的分析会报错：
cellchat@DB <- CellChatDB
有的时候我们并不想分析所有的细胞通讯，例如我们只关心Secreted Signaling，或者我们只相信经过KEGG分析验证过的细胞通讯，为了节约运行空间和时间，我们可以使用subsetDB()函数来取数据库的子集，简单查看一下数据库的数据结构：
dplyr::glimpse(CellChatDB$interaction)
Rows: 1,939
Columns: 11
$ interaction_name <chr> "TGFB1_TGFBR1_TGFBR2", "TGFB2_TGFBR1_TGFBR2", "TGFB3_TGFBR1_TGFBR2", "TGFB1_ACVR1B_TGFB…
$ pathway_name <chr> "TGFb", "TGFb", "TGFb", "TGFb", "TGFb", "TGFb", "TGFb", "TGFb", "TGFb", "TGFb", "TGFb",…
$ ligand <chr> "TGFB1", "TGFB2", "TGFB3", "TGFB1", "TGFB1", "TGFB2", "TGFB2", "TGFB3", "TGFB3", "TGFB1…
$ receptor <chr> "TGFbR1_R2", "TGFbR1_R2", "TGFbR1_R2", "ACVR1B_TGFbR2", "ACVR1C_TGFbR2", "ACVR1B_TGFbR2…
$ agonist <chr> "TGFb agonist", "TGFb agonist", "TGFb agonist", "TGFb agonist", "TGFb agonist", "TGFb a…
$ antagonist <chr> "TGFb antagonist", "TGFb antagonist", "TGFb antagonist", "TGFb antagonist", "TGFb antag…
$ co_A_receptor <chr> "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "",…
$ co_I_receptor <chr> "TGFb inhibition receptor", "TGFb inhibition receptor", "TGFb inhibition receptor", "TG…
$ evidence <chr> "KEGG: hsa04350", "KEGG: hsa04350", "KEGG: hsa04350", "PMID: 27449815", "PMID: 27449815…
$ annotation <chr> "Secreted Signaling", "Secreted Signaling", "Secreted Signaling", "Secreted Signaling",…
$ interaction_name_2 <chr> "TGFB1 - (TGFBR1+TGFBR2)", "TGFB2 - (TGFBR1+TGFBR2)", "TGFB3 - (TGFBR1+TGFBR2)", "TGFB1…
我们只选择Secreted Signaling相关的细胞间通讯信息：
subsetDB(CellChatDB.human, search = 'Secreted Signaling')
当然我们能做的绝不仅这些，我们可以通过更改subsetDB()函数的key参数来实现任何标准的筛选，例如我们想筛选出evidence为KEGG的细胞间通讯：
subsetDB(CellChatDB.human, search = '^KEGG', key = 'evidence') #regular expression here!!!
2.4 细胞通讯鉴定
在这个部分我们为了降低运行的内存和时间消耗，我们会抽取部分细胞进行分析，分析的流程也很好理解：首先，鉴定出细胞中高表达的受体和配体编码基因，然后再依据这个表达谱构建细胞之间的受体与配体作用对。和单细胞分析一样，这里我们也可以使用 future 包来进行并行计算，加速我们的分析过程。
library(future)
plan(strategy = 'multiprocess', workers = 4)
#subset
cellchat <- subsetData(cellchat)
为什么这里要做subset，因为为了节省运算时间和空间，通过这一步，我们后面的分析将只关注于与细胞通讯有关的基因（从数据框中提取出来的），所以在这里针对于表达矩阵的基因取了子集。
cellchat <- identifyOverExpressedGenes(cellchat)
cellchat <- identifyOverExpressedInteractions(cellchat)
上面两个函数分别鉴定除了每个细胞群中高表达的细胞通讯相关基因和细胞群之间根据细胞通讯相关基因的表达情况（identifyOverExpressedGenes()）最终鉴定出的细胞间通讯关系（identifyOverExpressedInteractions）。
cellchat <- computeCommunProb(cellchat)
cellchat <- filterCommunication(cellchat, min.cells = 10)
注意到为了获得更可信的细胞间通讯，方便后续的验证，这里首先对每个通讯计算了相应的概率值，同时进行了基于每个细胞类群中支持该通讯的细胞数量的过滤。
需要注意的是，预测出的细胞间通讯的数量是和每个细胞群的基因表达量有关的，所以对于基因平均表达量的计算就显得很重要，在默认情况下，computeCommunProb()函数使用的是trimean方法，该方法默认如果一群细胞里面不足25%的细胞表达某个基因的话，这个基因在该群细胞里面的平均表达量就是0。不过你可以通过设置type = "truncatedMean"和trim=来自己指定这个阈值，比如trim=0.1就表示阈值为10%。显然默认方法trimean能帮我们筛选出更少的、更可信的细胞间通讯。此外，考虑到细胞数更多的细胞群之间的细胞通讯信号往往会更强，为了去除不同细胞群之间的细胞数量差异，我们可以尝试使用population.size = TRUE来辅助我们发现稀有细胞类群之间的细胞通讯。你可以通过computeAveExpr(cellchat, features = c("CXCL12","CXCR4"), type = "truncatedMean", trim = 0.1)来提取你所感兴趣的细胞通讯相关基因的平均表达量信息。
2.5 细胞通讯与信号通路关系的构建
细胞之间的通讯往往会是信号通路的重要组成部分，把细胞通讯放到信号通路中进行理解可能会更有利于我们理解生物学过程。
cellchat <- computeCommunProbPathway(cellchat)
每对受配体的细胞间通讯网络会被存储到net的slot中，而每个信号通过的细胞间通讯网络信息将会被存储到netP网络中。
2.6 细胞通讯网络的构建
进一步的，我们将细胞间的通讯进行整合，就能构建出细胞间的通讯网络。
cellchat <- aggregateNet(cellchat)
当然，你可能只关心部分细胞之间的通讯，你可以通过指定 信号发出细胞 和 信号作用细胞 来进行个性化的分析：
?aggregateNet
#output
aggregateNet(
object,
sources.use = NULL,
targets.use = NULL,
signaling = NULL,
pairLR.use = NULL,
remove.isolate = TRUE,
thresh = 0.05,
return.object = TRUE
)
也就是这里的 sources.use 和 targets.use。那么指定成什么呢？这个时候metadata信息的作用就来了，指定的就是前面 group.by 所包含的细胞类群信息。
至此，上游分析已经算是结束，下面就是如何对信息进行展示，即 可视化。
