单细胞数据分析实战:用CellChat解析肿瘤微环境中的细胞通讯差异(附完整代码)

肿瘤微环境从来都不是一个静态的战场。当我们通过单细胞测序技术,将成千上万个细胞的身份与状态逐一解码后,一个更深层的问题便浮出水面:这些细胞之间,究竟在“交谈”什么?是谁在发出指令,谁又在默默响应?这些跨越细胞边界的分子对话,如何塑造了肿瘤的免疫抑制、血管新生或转移前生态位的形成?理解细胞间通讯,就是理解肿瘤微环境作为一个“社会系统”如何运作的关键。

CellChat,作为一个基于R语言的开源工具包,为我们提供了一套强大的数学框架和优雅的可视化方案,来系统性地推断、量化和比较细胞间的信号网络。它不仅仅能告诉你“细胞A和细胞B有交流”,更能揭示交流的强度、使用的“语言”(信号通路)、以及在不同生理或病理状态下,这场对话发生了怎样的戏剧性转变。对于致力于肿瘤免疫、微环境异质性或发育生物学的研究者而言,掌握CellChat意味着你手中的单细胞数据,能从一份细胞“花名册”,升级为一张动态的细胞“社交关系图谱”。

本文旨在为生物信息学分析人员和肿瘤研究者提供一份从零开始的实战指南。我们将绕过冗长的理论铺垫,直接切入核心操作,手把手带你完成从数据准备、CellChat对象创建、到差异通讯分析及高级可视化的全流程。文中每一段代码都经过实战检验,并附上了我本人在分析过程中踩过的“坑”及解决方案。无论你是想比较癌与原发灶的通讯差异,还是探索治疗前后微环境对话的改变,这篇文章都将为你提供一个坚实、可复现的起点。

1. 环境准备与数据质控:为CellChat分析奠定基石

在启动任何细胞通讯分析之前,确保你的单细胞数据已经过严格的质控与标准的预处理,是成功的第一步。CellChat虽然强大,但它对输入数据的质量非常敏感。一个常见的误区是,直接将原始的Seurat对象丢给CellChat,结果往往会在后续的概率计算或可视化步骤中遇到各种报错。

首先,你需要一个稳定的R环境。我强烈建议使用R 4.0以上的版本,并在分析前设置好工作目录和包管理。

# 设置工作目录,确保所有数据文件路径正确
setwd("~/your_project_path/")
# 清空环境变量,避免旧对象干扰
rm(list = ls())
# 设置随机种子,保证结果可重复
set.seed(12345)

接下来是核心R包的安装与加载。CellChat依赖于一系列生物信息学包,确保它们被正确安装是关键。

# 安装必要R包(如果尚未安装)
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install(c("Seurat", "Matrix", "dplyr", "ggplot2", "patchwork"))
install.packages("devtools")
devtools::install_github("sqjin/CellChat")
# 加载所有需要的库
library(Seurat)
library(CellChat)
library(patchwork)
library(ggplot2)
library(dplyr)
library(Matrix)

注意CellChat包通过GitHub安装,有时会因网络问题失败。如果遇到devtools::install_github报错,可以尝试先安装remotes包,或从镜像源下载源码包进行本地安装。

现在,让我们加载经过注释的单细胞数据。假设你已经有了一个Seurat对象 seurat_obj,其中包含了细胞类型注释信息(存储在seurat_obj$celltype中)和样本分组信息(例如seurat_obj$group,包含“Tumor”和“Normal”等)。数据质控的重点在于:

  • 基因表达矩阵:确保使用的是归一化后的数据(如RNA assay下的data槽),而非原始计数。
  • 细胞注释:用于分组的细胞类型标签必须是因子(factor)类型,且类别不宜过多过细,否则会导致信号过于稀疏。通常建议将稀有细胞群(细胞数<10)合并或过滤。
  • 样本完整性:如果你要比较不同组别(如原发灶vs转移灶),请确保每个组别内都包含足够多的细胞和主要的细胞类型,否则差异分析会缺乏统计效力。

一个稳健的数据准备流程可以这样进行:

# 假设 seurat_obj 是你的Seurat对象
# 1. 检查并确保细胞类型注释是因子
if (!is.factor(seurat_obj$celltype)) {
  seurat_obj$celltype <- as.factor(seurat_obj$celltype)
}
# 2. 过滤掉细胞数过少的细胞类型(例如少于20个细胞)
celltype_counts <- table(seurat_obj$celltype)
keep_celltypes <- names(celltype_counts)[celltype_counts >= 20]
seurat_obj <- subset(seurat_obj, subset = celltype %in% keep_celltypes)
# 3. 重新设置因子水平,去除被过滤的空水平
seurat_obj$celltype <- droplevels(seurat_obj$celltype)
# 4. 查看最终的细胞类型分布
print(table(seurat_obj$celltype, seurat_obj$group))

完成这些步骤后,你的数据就为CellChat分析做好了准备。记住,清晰的细胞注释和合理的细胞类型合并策略,是后续获得生物学可解释结果的前提。

2. 构建CellChat对象与基础通讯网络推断

有了准备好的Seurat对象,创建CellChat对象是一个直接的过程。但这一步中隐藏着几个决定分析深度与方向的关键选择。

创建CellChat对象的核心函数是createCellChat。你需要指定用于细胞分组的元数据列。这里,我们使用celltype列。

# 创建CellChat对象,以‘celltype’作为细胞分组依据
cellchat <- createCellChat(object = seurat_obj,
                           group.by = "celltype",
                           assay = "RNA") # 默认是‘RNA’,如果你的数据assay名不同,需指定
# 查看对象基本信息
cellchat

创建对象后,我们需要为其“装配”一个配体-受体数据库。CellChat内置了CellChatDB.mouse(小鼠)和CellChatDB.human(人)两个数据库。选择正确的物种至关重要。

# 使用人类配体-受体数据库
CellChatDB <- CellChatDB.human
# 查看数据库的结构和分类
showDatabaseCategory(CellChatDB)

showDatabaseCategory函数会输出一个表格,展示数据库包含的信号通路分类,如:

分类描述相互作用对数量
Secreted Signaling分泌信号约 2000 对
ECM-Receptor细胞外基质-受体约 500 对
Cell-Cell Contact细胞-细胞接触约 300 对

你可以选择使用全部数据库,也可以根据研究背景使用子集。例如,如果你只关注分泌信号,可以这样做:

# 使用全部数据库(默认,最全面)
CellChatDB.use <- CellChatDB
# 或者,仅使用‘Secreted Signaling’分类
# CellChatDB.use <- subsetDB(CellChatDB, search = "Secreted Signaling")
# 将选定的数据库赋给CellChat对象
cellchat@DB <- CellChatDB.use

接下来是核心的三步曲:识别过表达基因、识别过表达相互作用、并(可选)利用蛋白质互作网络进行数据投射。这里最容易出问题的是计算资源。单细胞数据量巨大时,identifyOverExpressedInteractions可能会消耗大量内存和时间。

# 1. 识别在每个细胞群中过表达的基因
cellchat <- identifyOverExpressedGenes(cellchat)
# 2. 识别过表达的配体-受体相互作用(此步计算量较大)
# 建议在服务器或配置较高的电脑上运行,可使用多线程加速
future::plan("multiprocess", workers = 4) # 设置并行,workers数根据CPU核心数调整
cellchat <- identifyOverExpressedInteractions(cellchat)
# 3. (可选)将基因表达数据投射到蛋白质互作网络,以考虑信号通路层面的关联
cellchat <- projectData(cellchat, PPI.human) # 使用人类PPI网络

完成以上步骤后,就可以开始推断细胞通讯概率了。computeCommunProb函数是CellChat的引擎,它基于表达丰度和统计学检验,计算每对细胞类型之间通过每个配体-受体对进行通讯的概率。

# 计算细胞间通讯概率
# raw.use = TRUE 表示使用原始数据(未进行投射的数据),计算更精确但稍慢
# 如果数据量极大,可设置为 FALSE 以使用投射后的数据加速
cellchat <- computeCommunProb(cellchat, raw.use = TRUE)
# 过滤掉可能不可靠的通讯(例如,仅在某个细胞群中极少数细胞表达的相互作用)
cellchat <- filterCommunication(cellchat, min.cells = 10)
# 将相互作用聚合到信号通路水平
cellchat <- computeCommunProbPathway(cellchat)
# 最后,聚合细胞类型级别的网络
cellchat <- aggregateNet(cellchat)

运行完aggregateNet,一个基础的细胞通讯网络就构建完成了。你可以立即进行一些整体可视化,对全局的通讯模式有一个直观认识。

# 可视化整体相互作用数量
par(mfrow = c(1,2))
netVisual_circle(cellchat@net$count, vertex.weight = table(cellchat@idents),
                 weight.scale = TRUE, label.edge = FALSE,
                 title.name = "Number of interactions")
# 可视化整体相互作用强度(权重)
netVisual_circle(cellchat@net$weight, vertex.weight = table(cellchat@idents),
                 weight.scale = TRUE, label.edge = FALSE,
                 title.name = "Interaction strength/weights")

这两张图能迅速告诉你哪些细胞类型是“社交达人”(发出或接收大量信号),以及信号流的总体强度分布。如果发现某个细胞类型完全孤立,可能需要回头检查其注释是否正确,或者是否因细胞数太少而被过度过滤。

3. 深度解析:可视化特定信号通路与模式识别

基础网络图给了我们一个宏观视角,但真正的生物学洞见往往藏在特定的信号通路中。CellChat提供了极其丰富的可视化函数,允许我们从多个维度“解剖”通讯网络。

假设我们对CXCL趋化因子信号通路在肿瘤微环境中的作用感兴趣。首先,我们可以用层次图来展示该通路中信号的流向。

# 定义感兴趣的通路
pathways.show <- c("CXCL")
# 选择接收信号的细胞群(这里以索引号指定,也可用细胞类型名)
# 假设我们想关注该信号如何流向T细胞(假设其索引为1,2,3)
vertex.receiver <- c(1, 2, 3)
# 绘制层次图
netVisual_aggregate(cellchat, signaling = pathways.show,
                    vertex.receiver = vertex.receiver,
                    layout = "hierarchy")

层次图清晰地展示了信号从发送者到接收者的层级关系。而和弦图则能更优美地展示所有细胞类型在该通路下的互作关系,特别适合展示复杂的多对多通讯。

# 绘制和弦图
netVisual_aggregate(cellchat, signaling = pathways.show,
                    layout = "chord")

提示:和弦图在细胞类型较多时可能会显得拥挤。可以通过signaling.name参数修改标题,或使用small.gapbig.gap参数调整布局。

除了看单个通路,我们更关心细胞在整体网络中的角色。CellChat的网络中心性分析可以量化每个细胞类型作为信号发送者、接收者或中介者的重要性。

# 计算网络中心性指标
cellchat <- netAnalysis_computeCentrality(cellchat, slot.name = "netP")
# 可视化特定信号通路中细胞角色的网络图
netAnalysis_signalingRole_network(cellchat, signaling = pathways.show,
                                  width = 20, height = 5, font.size = 8)

更进一步,我们可以使用模式识别算法,将具有相似通讯模式的细胞类型或信号通路聚类在一起。这能帮助我们发现潜在的、功能协同的细胞群体或信号模块。

# 识别传入信号模式
# 首先通过肘部法则选择最佳模式数K
selectK(cellchat, pattern = "incoming")
# 假设根据上图选择 K=3
nPatterns <- 3
cellchat <- identifyCommunicationPatterns(cellchat, pattern = "incoming",
                                          k = nPatterns, height = 15)
# 用河流图展示模式
netAnalysis_river(cellchat, pattern = "incoming")
# 用点图展示模式与通路的关联
netAnalysis_dot(cellchat, pattern = "incoming")

河流图直观地展示了不同模式下的信号流强度分布,而点图则量化了每个信号通路对每个模式的贡献度。这些分析能帮你回答诸如“哪些细胞类型接收了相似的信号组合?”或“哪些信号通路倾向于协同作用?”等问题。

4. 差异细胞通讯分析:比较不同条件下的对话差异

这才是大多数肿瘤微环境研究的核心:比较不同生理或病理状态下的细胞通讯有何不同。例如,比较原发性肿瘤和转移灶,或者治疗前与治疗后的样本。CellChat提供了系统的方法进行这种比较分析。

首先,你需要为每个条件独立创建CellChat对象。假设我们已经有了代表“Primary_Tumor”和“Metastasis”的两个Seurat对象,并分别创建了cellchat_primarycellchat_meta

# 为每个样本组创建独立的CellChat对象(示例代码,需对每个组别执行)
# cellchat_primary <- createCellChat(seurat_primary, group.by = "celltype")
# cellchat_meta <- createCellChat(seurat_meta, group.by = "celltype")
# ...(对每个对象执行第2章中的数据库设置、过表达识别、概率计算等步骤)

# 将对象存入列表
object.list <- list(Primary = cellchat_primary, Metastasis = cellchat_meta)
# 合并CellChat对象,为比较分析做准备
cellchat.merged <- mergeCellChat(object.list, add.names = names(object.list))

合并后,我们可以从多个层面进行差异比较。

1. 整体相互作用数量与强度的差异: 这能告诉你两个条件下,细胞社交网络的活跃度发生了全局性还是局部性的变化。

# 比较相互作用数量
gg1 <- netVisual_diffInteraction(cellchat.merged, comparison = c(1,2),
                                 weight.scale = TRUE, measure = "count")
# 比较相互作用强度(权重)
gg2 <- netVisual_diffInteraction(cellchat.merged, comparison = c(1,2),
                                 weight.scale = TRUE, measure = "weight")
gg1 + gg2

热图中,红色表示在第二个条件(如转移灶)中增强的相互作用,蓝色则表示减弱。

2. 信号通路信息流的差异: 哪些信号通路在条件间发生了显著上调或下调?

# 识别并可视化差异最大的信号通路
gg <- rankNet(cellchat.merged, mode = "comparison",
              stacked = FALSE, # 设置为TRUE可得到堆叠条形图
              font.size = 4)
print(gg)

这张图按差异程度排序,清晰展示了在转移灶中相对活跃或沉默的信号通路。

3. 保守与特异性的信号通路: 哪些信号通路是两种条件共有的(保守),哪些是某一种条件特有的?

# 计算功能相似性,识别保守通路
rankSimilarity(cellchat.merged, type = "functional")
# 识别并提取保守和条件特异的信号通路
# 这里以保守通路为例,提取在两种条件下都显著活跃的通路
conserved_pathways <- identifyConservedPatterns(cellchat.merged, type = "functional")

4. 特定配体-受体对的差异: 当你有一个明确的假设时,可以深入查看特定细胞类型间,通过特定配体-受体对的通讯变化。

# 气泡图展示特定发送者与接收者之间所有相互作用的差异
netVisual_bubble(cellchat.merged,
                 sources.use = c("Macrophage"), # 发送者细胞类型
                 targets.use = c("CD8_Tcell", "Treg"), # 接收者细胞类型
                 comparison = c(1,2), # 比较第1组和第2组
                 angle.x = 90)

这张气泡图会展示巨噬细胞向CD8 T细胞或Treg细胞发送的所有信号中,哪些在转移灶中显著增强(红色)或减弱(蓝色)。这是将高通量分析与具体生物学假设连接起来的强大工具。

5. 高级技巧、结果解读与避坑指南

经过上述步骤,你已经完成了一套完整的CellChat分析。但在将结果转化为生物学故事时,还有一些高级技巧和注意事项。

结果解读的核心在于区分“相关性”与“因果性”。CellChat推断的是潜在的通讯概率,基于的是配体与受体基因的共表达。它强烈暗示了通讯的可能性,但并非直接证明。因此,在阐述结论时,应使用“可能介导”、“提示存在……信号”等谨慎措辞,并结合已知文献进行佐证。

高级可视化定制能让你的图表更专业。例如,自定义颜色方案以匹配你文章中的细胞类型配色。

# 定义细胞类型颜色向量
colors_celltype <- c("B_cell" = "#1f77b4",
                     "CD4_Tcell" = "#ff7f0e",
                     "CD8_Tcell" = "#2ca02c",
                     "Macrophage" = "#d62728",
                     "Fibroblast" = "#9467bd")
# 在可视化函数中通过color.use参数调用
netVisual_aggregate(cellchat, signaling = pathways.show,
                    layout = "circle",
                    color.use = colors_celltype)

常见报错与解决方案:

  1. Error in subsetData(cellchat) 或后续步骤报错:最常见的原因是创建CellChat对象时,指定的group.by列包含NA值或不是因子。务必在创建对象前检查并清理元数据。
  2. identifyOverExpressedInteractions 运行极慢或内存不足:单细胞数据量过大(如超过5万个细胞)。解决方案:
    • 使用future::plan("multisession", workers = 6)设置更多并行workers。
    • computeCommunProb中尝试设置raw.use = FALSE,使用投射后的数据(会损失一些精度但快很多)。
    • 考虑先对细胞进行适度的下采样(如每个细胞类型随机取2000个细胞),进行探索性分析。
  3. 可视化图形中文字重叠或布局混乱:调整图形尺寸和字体大小。几乎所有netVisual_*函数都有width, height, font.size, vertex.label.cex等参数供调整。
  4. 差异分析时出现“对象长度不一致”错误:确保用于比较的多个CellChat对象中的细胞类型基本一致。如果某个细胞类型只在一个条件中出现,需要使用liftCellChat函数进行“对齐”。

最后,分析的可复现性至关重要。养成保存中间RDS文件和记录完整R脚本的习惯。一个推荐的项目目录结构如下:

your_project/
├── data/
│   ├── raw_seurat.rds
│   └── processed_seurat.rds
├── script/
│   ├── 01_data_preprocessing.R
│   ├── 02_cellchat_analysis.R
│   └── 03_figure_generation.R
├── results/
│   ├── cellchat_primary.rds
│   ├── cellchat_metastasis.rds
│   └── figures/ # 存放所有输出图片
└── README.md # 记录分析流程和参数

将CellChat分析整合进你的单细胞研究流程,绝非简单的“跑一个流程”。它要求你对生物学问题有深刻理解,对数据质量有严格把控,并对分析结果的统计本质有清醒认识。我最初几次使用时,曾被那些绚丽的网络图吸引,却忽略了统计效力不足带来的假阳性。后来才明白,在设置min.cells过滤阈值、解释差异分析p值时,必须结合样本的生物学重复和细胞数量来综合判断。现在,我更倾向于将CellChat的结果视为生成假设的工具,它为我指明了哪些细胞对话最值得用实验手段去进一步验证。当你看到巨噬细胞与T细胞之间的某条通路在耐药样本中特异性增强时,那种将计算预测与潜在机制联系起来的瞬间,正是生物信息学分析最令人兴奋的时刻。

Logo

DAMO开发者矩阵,由阿里巴巴达摩院和中国互联网协会联合发起,致力于探讨最前沿的技术趋势与应用成果,搭建高质量的交流与分享平台,推动技术创新与产业应用链接,围绕“人工智能与新型计算”构建开放共享的开发者生态。

更多推荐