您现在的位置是：首页 > 云平台

当前栏目

CellChat学习笔记【一】——通讯网络构建

网络笔记学习构建通讯 CellChat

2023-06-13 09:13:22 时间

细胞通讯是指细胞与细胞之间的联系。细胞和人类一样，多细胞生物的很多细胞会相互作用，形成“细胞社会”，在这个社会里，细胞与细胞之间会发生相互作用和信息的传递，细胞建立通讯联络是必需的。如生物体的生长发育、分化、各种组织器官的形成、组织的维持以及它们各种生理活动的协调。经典的例子莫过于 神经细胞之间的神经递质的传递与接收。

细胞与细胞之间的通讯有三种形式：

细胞之间的直接接触；
细胞之间通过其间的胞外基质相互联系；
细胞之间通过分子间的相互作用产生联系。

在这个过程中，“受体-配体”的概念十分重要，什么是受体，什么是配体？受体与配体之间结合的结果是受体被激活，并产生受体激活后续信号传递的基本步骤。在单细胞分析当中，不同细胞类型之间的通讯可能会对某些生物学过程具有重要意义，因此利用单细胞数据进行细胞通讯分析是单细胞高级分析的一大重点。

CellChat

CellChat是一款2021年发表于Nature Communications的单细胞细胞通讯分析工具。

CellChat上游分析

安装 CellChat

devtools::install_github("sqjin/CellChat")

细胞通讯分析

这里我们使用 pbmc3k 的公共数据集为例来一起学习如何利用 CellChat 进行细胞通讯分析。

数据输入

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的分析。

创建CellChat对象

CellChat对象和Seurat对象很像，我们可以这样来进行创建：

cellchat <- createCellChat(object = data,
                           meta = meta,
                           group.by = 'seurat_annotations')

## The cell groups used for CellChat analysis are  Naive CD4 T Memory CD4 T CD14+ Mono B CD8 T FCGR3A+ Mono NK DC Platelet

这里的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)来查看当前的细胞分组信息。

选择受体配体数据库

CellChat 有一个专门的数据库，叫做CellChatDB，这个数据库是 CellChat 的作者们通过阅读大量文献，手动整理出来的“受体-配体”对，目前有人、鼠以及斑马鱼的版本。其中人的叫做 CellChatDB.human，鼠的叫做 CellChatDB.mouse，斑马鱼的叫做 CellChatDB.zebrafish。关于这两个数据库中具体的“受体-配体”对信息可以通过showDatabaseCategory()函数获得，在这里就不赘述了。

CellChat数据库

mouse

2,021 validated molecular interactions, including 60% of secrete autocrine/paracrine signaling interactions, 21% of extracellular matrix (ECM)-receptor interactions and 19% of cell-cell contact interactions.

human

1,939 validated molecular interactions, including 61.8% of paracrine/autocrine signaling interactions, 21.7% of extracellular matrix (ECM)-receptor interactions and 16.5% of cell-cell contact interactions.

加载数据库：

CellChatDB <- CellChatDB.human

来看一下数据库的结构（option）：

head(CellChatDB$interaction)

##                        interaction_name pathway_name ligand      receptor
## TGFB1_TGFBR1_TGFBR2 TGFB1_TGFBR1_TGFBR2         TGFb  TGFB1     TGFbR1_R2
## TGFB2_TGFBR1_TGFBR2 TGFB2_TGFBR1_TGFBR2         TGFb  TGFB2     TGFbR1_R2
## TGFB3_TGFBR1_TGFBR2 TGFB3_TGFBR1_TGFBR2         TGFb  TGFB3     TGFbR1_R2
## TGFB1_ACVR1B_TGFBR2 TGFB1_ACVR1B_TGFBR2         TGFb  TGFB1 ACVR1B_TGFbR2
## TGFB1_ACVR1C_TGFBR2 TGFB1_ACVR1C_TGFBR2         TGFb  TGFB1 ACVR1C_TGFbR2
## TGFB2_ACVR1B_TGFBR2 TGFB2_ACVR1B_TGFBR2         TGFb  TGFB2 ACVR1B_TGFbR2
##                          agonist      antagonist co_A_receptor
## TGFB1_TGFBR1_TGFBR2 TGFb agonist TGFb antagonist              
## TGFB2_TGFBR1_TGFBR2 TGFb agonist TGFb antagonist              
## TGFB3_TGFBR1_TGFBR2 TGFb agonist TGFb antagonist              
## TGFB1_ACVR1B_TGFBR2 TGFb agonist TGFb antagonist              
## TGFB1_ACVR1C_TGFBR2 TGFb agonist TGFb antagonist              
## TGFB2_ACVR1B_TGFBR2 TGFb agonist TGFb antagonist              
##                                co_I_receptor       evidence         annotation
## TGFB1_TGFBR1_TGFBR2 TGFb inhibition receptor KEGG: hsa04350 Secreted Signaling
## TGFB2_TGFBR1_TGFBR2 TGFb inhibition receptor KEGG: hsa04350 Secreted Signaling
## TGFB3_TGFBR1_TGFBR2 TGFb inhibition receptor KEGG: hsa04350 Secreted Signaling
## TGFB1_ACVR1B_TGFBR2 TGFb inhibition receptor PMID: 27449815 Secreted Signaling
## TGFB1_ACVR1C_TGFBR2 TGFb inhibition receptor PMID: 27449815 Secreted Signaling
## TGFB2_ACVR1B_TGFBR2 TGFb inhibition receptor PMID: 27449815 Secreted Signaling
##                          interaction_name_2
## TGFB1_TGFBR1_TGFBR2 TGFB1 - (TGFBR1+TGFBR2)
## TGFB2_TGFBR1_TGFBR2 TGFB2 - (TGFBR1+TGFBR2)
## TGFB3_TGFBR1_TGFBR2 TGFB3 - (TGFBR1+TGFBR2)
## TGFB1_ACVR1B_TGFBR2 TGFB1 - (ACVR1B+TGFBR2)
## TGFB1_ACVR1C_TGFBR2 TGFB1 - (ACVR1C+TGFBR2)
## TGFB2_ACVR1B_TGFBR2 TGFB2 - (ACVR1B+TGFBR2)

最后，千万不要忘了把数据库嵌入到 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", "TGFB…
## $ pathway_name       <chr> "TGFb", "TGFb", "TGFb", "TGFb", "TGFb", "TGFb", "TG…
## $ ligand             <chr> "TGFB1", "TGFB2", "TGFB3", "TGFB1", "TGFB1", "TGFB2…
## $ receptor           <chr> "TGFbR1_R2", "TGFbR1_R2", "TGFbR1_R2", "ACVR1B_TGFb…
## $ agonist            <chr> "TGFb agonist", "TGFb agonist", "TGFb agonist", "TG…
## $ antagonist         <chr> "TGFb antagonist", "TGFb antagonist", "TGFb antagon…
## $ co_A_receptor      <chr> "", "", "", "", "", "", "", "", "", "", "", "", "",…
## $ co_I_receptor      <chr> "TGFb inhibition receptor", "TGFb inhibition recept…
## $ evidence           <chr> "KEGG: hsa04350", "KEGG: hsa04350", "KEGG: hsa04350…
## $ annotation         <chr> "Secreted Signaling", "Secreted Signaling", "Secret…
## $ interaction_name_2 <chr> "TGFB1 - (TGFBR1+TGFBR2)", "TGFB2 - (TGFBR1+TGFBR2)…

我们只选择Secreted Signaling相关的细胞间通讯信息：

subsetDB(CellChatDB.human, search = 'Secreted Signaling')

当然我们能做的绝不仅这些，我们可以通过更改subsetDB()函数的key参数来实现任何标准的筛选，例如我们想筛选出evidence为KEGG的细胞间通讯：

subsetDB(CellChatDB.human, search = '^KEGG', key = 'evidence') #regular expression here!!!

细胞通讯鉴定

在这个部分我们为了降低运行的内存和时间消耗，我们会抽取部分细胞进行分析，分析的流程也很好理解：首先，鉴定出细胞中高表达的受体和配体编码基因，然后再依据这个表达谱构建细胞之间的受体与配体作用对。和单细胞分析一样，这里我们也可以使用 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)来提取你所感兴趣的细胞通讯相关基因的平均表达量信息。

细胞通讯与信号通路关系的构建

细胞之间的通讯往往会是信号通路的重要组成部分，把细胞通讯放到信号通路中进行理解可能会更有利于我们理解生物学过程。

cellchat <- computeCommunProbPathway(cellchat)

每对受配体的细胞间通讯网络会被存储到net的slot中，而每个信号通过的细胞间通讯网络信息将会被存储到netP网络中。

细胞通讯网络的构建

进一步的，我们将细胞间的通讯进行整合，就能构建出细胞间的通讯网络。

cellchat <- aggregateNet(cellchat)

当然，你可能只关心部分细胞之间的通讯，你可以通过指定 信号发出细胞 和 信号作用细胞 来进行个性化的分析：

?aggregateNet

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 所包含的细胞类群信息。

至此，上游分析已经结束，后面我们将继续分享如何理解这些分析结果，并对其进行可视化。

猜你喜欢

【C 语言】结构体 ( 结构体数组作为函数参数 | 数组在栈内存创建 )
Linux带给你的网吧新体验（linux网吧）
课题介绍-上篇｜2023腾讯犀牛鸟精英人才计划
node.js中的console.error方法使用说明
android图片文字识别器,图片转换文字识别器[通俗易懂]
实验性的测试标记揭示了Win10 Sun Valley的浮动菜单
densenet121网络结构_resnet网络结构详解
MySQL查询如何实现一列字段相加（mysql 一列字段相加）
快速导出Oracle视图数据（oracle导出视图）
权限Linux文件权限掌握关键！（linux什么是文件）
TOMCAT启动报错：严重: Error filterStart
C++文件输入输出流（无师自通）
java锁：第四章：读写锁
Linux 升级 JDK：突破技术壁垒（linux升级jdk）
Linux下如何调整系统语言（linux设置系统语言）
【基于蜂窝网络的物联网开发】你需要了解的若干关于物联网卡的常规知识
【APUE】文件 I/O 操作
Linux下文件快速查找的秘密（文件查找linux）
SQL Server之上求下技术：精准轻松编写（sqlserver上下标）
解锁虚拟空间之旅快速安装Redis（虚拟空间安装 redis）
windows下安装mongodb以及node.js连接mongodb实例
【Kotlin】Kotlin 单例 ( 懒汉式与恶汉式 | Java 单例 | Kotlin 单例 | 对象声明 | 伴生对象 | get 方法 | ? 与 !! 判空 )
ORA-13045: invalid compatibility flag ORACLE 报错故障修复远程处理
火车头采集常用的纯正则过滤

zl程序教程