实验课6 拟时序和细胞通讯分析 R代码整理¶

对应材料：

实验/课件/实验课6-拟时序和细胞通讯分析.pptx
回放/生物信息学（协和班）第13周星期4第8,9节_笔记.txt
老师板书整理

说明¶

这份整理以老师板书为主，结合课件和第13周实验课回放做了必要补全。
这节实验课分成两条主线：
Monocle3 拟时序分析
CellChat 细胞通讯分析
课件本身比较简，真正有用的代码主线主要在板书和老师口头说明里。

一、老师对这节实验课的总体要求¶

这一节和前面的单细胞基础分析不太一样，老师在回放里一开始就强调了几件事：

这两类分析原本都不太打算给大家单独做，因为“比较虚”、主观性更强。
后来是因为同学想学，所以额外增加了一节实验课。
这两类分析在实际科研中使用时要非常谨慎，不能轻易把结果当成很硬的结论。
这节课更像“扩展分析”和“了解常用工具包思路”，不是那种最稳妥、最核心的基础流程。

老师口头最强调的风险点有两个：

1. 拟时序的风险点¶

最大的主观步骤是“起点/root 节点”的选择。
你选择的起点不同，整个伪时序方向就可能变化。
所以拟时序虽然图很好看，但“为什么这个点是起点”必须有生物学逻辑支撑。

2. 细胞通讯的风险点¶

细胞通讯依赖配体-受体数据库。
数据库本身的大小、收录范围、阈值设置，都会影响最终结果。
结果里可能有很多配体-受体对，但真正重要的常常只有其中一小部分。
如果没有后续实验验证，细胞通讯分析的结论也要谨慎看待。

一句话总结老师态度：

会用、会看图、会说流程可以；
但不要把它们当成“只要跑出来就一定可信”的分析。

二、Monocle3 拟时序分析：板书整理后的主干 R 代码¶

library(Seurat)
library(monocle3)
library(SingleCellExperiment)

pbmc <- readRDS("pbmc_annotated_for_downstream.rds")

expr_mat <- GetAssayData(
  pbmc,
  assay = "RNA",
  layer = "counts"
)

cell_metadata <- pbmc@meta.data
cell_metadata$cell_id <- rownames(cell_metadata)

gene_metadata <- data.frame(
  gene_short_name = rownames(expr_mat),
  row.names = rownames(expr_mat)
)

cds <- new_cell_data_set(
  expression_data = expr_mat,
  cell_metadata = cell_metadata,
  gene_metadata = gene_metadata
)

cds <- preprocess_cds(cds, num_dim = 20)
cds <- reduce_dimension(cds, reduction_method = "UMAP")
cds <- cluster_cells(cds)
cds <- learn_graph(cds)

cds <- order_cells(cds)

plot_cells(
  cds,
  color_cells_by = "new.cluster.ids",
  label_groups_by_cluster = FALSE,
  label_leaves = TRUE,
  label_branch_points = TRUE
)

plot_cells(
  cds,
  color_cells_by = "pseudotime",
  label_groups_by_cluster = FALSE,
  label_leaves = TRUE,
  label_branch_points = TRUE
)

三、Monocle3 代码对应老师板书的哪几块¶

1. 数据准备¶

对应第一张板书：

library(Seurat)
library(monocle3)
library(SingleCellExperiment)

pbmc <- readRDS("pbmc_annotated_for_downstream.rds")

老师这节课没有让大家从 Read10X() 和 Seurat 预处理开始重新做，而是直接给了一个已经预处理好的对象：

pbmc_annotated_for_downstream.rds

回放里老师也明确说了：

这次发给大家的数据已经是处理好的
直接从这里开始做拟时序和细胞通讯
当然你也可以用自己之前处理好的数据再走一遍
但课堂上建议先用统一给的数据

2. 从 Seurat 对象里取表达矩阵和元数据¶

对应第一、第二张板书：

expr_mat <- GetAssayData(pbmc, assay = "RNA", layer = "counts")

cell_metadata <- pbmc@meta.data
cell_metadata$cell_id <- rownames(cell_metadata)

gene_metadata <- data.frame(
  gene_short_name = rownames(expr_mat),
  row.names = rownames(expr_mat)
)

这部分本质上是在做格式转换：

Monocle3 不直接吃 Seurat 对象
所以要把表达矩阵、细胞信息、基因信息分别整理出来

老师在回放里专门说了：

这类代码“看起来很多”
但前面很大一部分其实都只是为了满足对象格式要求
真正的核心算法代码只占后面几行

3. 创建 `cell_data_set`¶

对应第二张板书左侧：

cds <- new_cell_data_set(
  expression_data = expr_mat,
  cell_metadata = cell_metadata,
  gene_metadata = gene_metadata
)

这一步是 Monocle3 的关键入口：

把整理好的表达矩阵、细胞元数据、基因元数据组合成 cds

4. 核心分析步骤¶

对应第二张板书左下和右侧：

cds <- preprocess_cds(cds, num_dim = 20)
cds <- reduce_dimension(cds, reduction_method = "UMAP")
cds <- cluster_cells(cds)
cds <- learn_graph(cds)
cds <- order_cells(cds)

老师在回放里对这部分讲得很明确：

前面都是“创建对象”
真正核心步骤就在这几行
主要逻辑就是：
预处理
降维
聚类
学习轨迹图
排序得到伪时序

也就是说，拟时序这部分最该记住的主线是：

准备Monocle3输入格式
-> 创建cds对象
-> preprocess_cds
-> reduce_dimension
-> cluster_cells
-> learn_graph
-> order_cells
-> plot_cells

5. `plot_cells()` 两种常见画法¶

对应第二张板书右侧：

plot_cells(
  cds,
  color_cells_by = "new.cluster.ids",
  ...
)

plot_cells(
  cds,
  color_cells_by = "pseudotime",
  ...
)

这两张图分别表示：

按 cluster 看轨迹结构
按 pseudotime 看早晚顺序

课件第 4 页也正好给了这两类示意图：

按聚类显示轨迹
按拟时序显示轨迹

四、Monocle3 相关要求：老师特别强调什么¶

1. 真正最主观的是 `order_cells()` 的起点选择¶

板书里写的是最简单形式：

cds <- order_cells(cds)

但老师在回放里专门强调：

这一步最关键、也最主观
根节点选得不一样，最后的伪时序方向会变
所以不能只看图漂亮不漂亮，要能解释为什么这里是“早期细胞”

2. 这部分更适合“会说流程”，不适合死背长代码¶

老师对这节课整体的态度很明确：

代码主线并不复杂
复杂的是前面的格式准备
所以复习时更应该抓住“输入是什么、对象怎么构建、核心步骤是什么、结果怎么看”

五、CellChat 细胞通讯分析：板书整理后的主干 R 代码¶

library(Seurat)
library(CellChat)
library(patchwork)

pbmc <- readRDS("pbmc_annotated_for_downstream.rds")
Idents(pbmc) <- "new.cluster.ids"

data.input <- GetAssayData(
  pbmc,
  assay = "RNA",
  layer = "data"
)

meta <- data.frame(
  labels = Idents(pbmc),
  row.names = colnames(pbmc)
)

cellchat <- createCellChat(
  object = data.input,
  meta = meta,
  group.by = "labels"
)

CellChatDB <- CellChatDB.human
cellchat@DB <- CellChatDB

cellchat <- subsetData(cellchat)
cellchat <- identifyOverExpressedGenes(cellchat, do.fast = FALSE)
cellchat <- identifyOverExpressedInteractions(cellchat)

cellchat <- computeCommunProb(cellchat)
cellchat <- filterCommunication(cellchat, min.cells = 10)
cellchat <- computeCommunProbPathway(cellchat)
cellchat <- aggregateNet(cellchat)

groupSize <- as.numeric(table(cellchat@idents))

netVisual_circle(
  cellchat@net$count,
  vertex.weight = groupSize,
  weight.scale = TRUE,
  label.edge = FALSE,
  title.name = "Number of interactions"
)

netVisual_circle(
  cellchat@net$weight,
  vertex.weight = groupSize,
  weight.scale = TRUE,
  label.edge = FALSE,
  title.name = "Interaction weights"
)

六、CellChat 代码对应老师板书的哪几块¶

1. 数据准备¶

对应第三张板书左侧：

library(Seurat)
library(CellChat)
library(patchwork)

pbmc <- readRDS("pbmc_annotated_for_downstream.rds")
Idents(pbmc) <- "new.cluster.ids"

这里老师直接把 cluster 身份设成：

Idents(pbmc) <- "new.cluster.ids"

也就是说，后面的细胞通讯分析，是以“已经注释好的细胞类型/cluster”作为分组基础的。

2. 取表达矩阵和分组信息¶

对应第三张板书左下：

data.input <- GetAssayData(pbmc, assay = "RNA", layer = "data")

meta <- data.frame(
  labels = Idents(pbmc),
  row.names = colnames(pbmc)
)

这里注意一个和拟时序不一样的点：

拟时序这里用的是 counts
CellChat 这里老师板书写的是 layer = "data"

也就是用归一化后的表达矩阵。

3. 创建 `CellChat` 对象并指定数据库¶

对应第三张板书右侧上半部分：

cellchat <- createCellChat(
  object = data.input,
  meta = meta,
  group.by = "labels"
)

CellChatDB <- CellChatDB.human
cellchat@DB <- CellChatDB

也就是说：

输入表达矩阵
输入细胞分组信息
指定使用人类数据库

4. 前置筛选步骤¶

对应第三张板书右侧中间：

cellchat <- subsetData(cellchat)
cellchat <- identifyOverExpressedGenes(cellchat, do.fast = FALSE)
cellchat <- identifyOverExpressedInteractions(cellchat)

老师在回放里对这一块说得很清楚：

这几步是 CellChat 的核心步骤之一
它本质上是在找高表达的配体和受体，以及高表达的相互作用对

5. 通讯概率、通路、网络聚合¶

对应第四张板书左侧：

cellchat <- computeCommunProb(cellchat)
cellchat <- filterCommunication(cellchat, min.cells = 10)
cellchat <- computeCommunProbPathway(cellchat)
cellchat <- aggregateNet(cellchat)

老师在回放里说得很明确：

CellChat 真正做运算的核心也就是这几步
其他很多代码都只是为了满足对象格式要求

也就是说，这一段最值得背的是：

subsetData
-> identifyOverExpressedGenes
-> identifyOverExpressedInteractions
-> computeCommunProb
-> filterCommunication
-> computeCommunProbPathway
-> aggregateNet

6. 两类最常见可视化¶

对应第四张板书左右两边：

groupSize <- as.numeric(table(cellchat@idents))

netVisual_circle(
  cellchat@net$count,
  ...
  title.name = "Number of interactions"
)

netVisual_circle(
  cellchat@net$weight,
  ...
  title.name = "Interaction weights"
)

这两张图分别对应课件第 7 页：

通讯数量 Number of interactions
通讯强度 Interaction weights

课件也专门解释了：

边越粗，说明预测到的相互作用数量越多，或者整体通讯强度越高

七、CellChat 相关要求：老师特别强调什么¶

1. 这类结果非常依赖数据库和阈值¶

回放里老师对这一块的提醒很强：

细胞通讯分析靠的是配体-受体数据库
数据库的内容和阈值会直接影响结果
所以不能把输出图当成“绝对真相”

2. 最后真正关心的，常常不是整张圈图，而是具体哪一对配体-受体最重要¶

老师回放里举的意思很清楚：

某两类细胞可能有很多互作
但真正值得深究的，往往是其中显著性最高、最合理的那一对或几对 ligand-receptor

所以真正做科研时，往往会进一步：

从中间结果里把具体互作对抓出来
再看是哪条通路、哪对配体-受体最关键
而不是停留在一张总览图上

3. 这部分比拟时序更依赖后续验证¶

老师的总体态度是：

可以做
可以看
但如果没有后续实验验证，结论一定要保守

八、这节实验课最该怎么复习¶

这节课不适合按“逐行背代码”的方式复习，更适合按下面两条主线背：

1. 拟时序主线¶

读入预处理好的Seurat对象
-> 提取counts矩阵、细胞元数据、基因元数据
-> 创建Monocle3的cds对象
-> preprocess_cds
-> reduce_dimension
-> cluster_cells
-> learn_graph
-> order_cells
-> 按cluster和按pseudotime画轨迹图

2. 细胞通讯主线¶

读入预处理好的Seurat对象
-> 设置分组身份
-> 提取归一化表达矩阵和meta
-> 创建CellChat对象
-> 指定人类数据库
-> 找高表达基因和高表达互作
-> 计算通讯概率
-> 过滤低支持结果
-> 按通路和整体网络聚合
-> 画通讯数量和通讯强度图

九、如果考试只让你写流程图或伪代码，最稳妥的写法¶

1. 拟时序分析¶

输入：已经完成单细胞基础分析和注释的Seurat对象

1. 提取表达矩阵、细胞元数据、基因元数据
2. 创建Monocle3对象
3. 进行预处理和降维
4. 对细胞重新聚类
5. 学习轨迹图结构
6. 选择起点并排序细胞
7. 绘制按cluster和按pseudotime显示的轨迹图

2. 细胞通讯分析¶

输入：已经完成单细胞基础分析和注释的Seurat对象

1. 提取归一化表达矩阵和细胞分组信息
2. 创建CellChat对象
3. 指定物种数据库
4. 筛选高表达基因和高表达互作对
5. 计算细胞间通讯概率
6. 过滤低支持结果
7. 按通路和整体网络聚合
8. 绘制通讯数量图和通讯强度图

十、一句提醒¶

这节实验课最重要的不是“我能不能把所有函数名一字不差背下来”，而是你要会解释：

为什么这两类分析都需要在单细胞基础分析之后再做
拟时序为什么主观，主观在哪一步
细胞通讯为什么依赖数据库和阈值
哪几步是真正核心计算
最终图应该怎么解释，哪些地方不能过度解读

实验课6 拟时序和细胞通讯分析 R代码整理¶

说明¶

一、老师对这节实验课的总体要求¶

1. 拟时序的风险点¶

2. 细胞通讯的风险点¶

二、Monocle3 拟时序分析：板书整理后的主干 R 代码¶

三、Monocle3 代码对应老师板书的哪几块¶

1. 数据准备¶

2. 从 Seurat 对象里取表达矩阵和元数据¶

3. 创建 cell_data_set¶

4. 核心分析步骤¶

5. plot_cells() 两种常见画法¶

四、Monocle3 相关要求：老师特别强调什么¶

1. 真正最主观的是 order_cells() 的起点选择¶

2. 这部分更适合“会说流程”，不适合死背长代码¶

五、CellChat 细胞通讯分析：板书整理后的主干 R 代码¶

六、CellChat 代码对应老师板书的哪几块¶

1. 数据准备¶

2. 取表达矩阵和分组信息¶

3. 创建 CellChat 对象并指定数据库¶

4. 前置筛选步骤¶

5. 通讯概率、通路、网络聚合¶

6. 两类最常见可视化¶

七、CellChat 相关要求：老师特别强调什么¶

1. 这类结果非常依赖数据库和阈值¶

2. 最后真正关心的，常常不是整张圈图，而是具体哪一对配体-受体最重要¶

3. 这部分比拟时序更依赖后续验证¶

八、这节实验课最该怎么复习¶

1. 拟时序主线¶

2. 细胞通讯主线¶

九、如果考试只让你写流程图或伪代码，最稳妥的写法¶

1. 拟时序分析¶

2. 细胞通讯分析¶

十、一句提醒¶

3. 创建 `cell_data_set`¶

5. `plot_cells()` 两种常见画法¶

1. 真正最主观的是 `order_cells()` 的起点选择¶

3. 创建 `CellChat` 对象并指定数据库¶