实验课5 单细胞分析 R代码整理

对应材料：

• 实验/课件/实验课5-单细胞分析-改后版.pptx
• 回放/生物信息学（协和班）第11周星期4第8,9节_笔记.txt
• 老师板书整理

本节考试要求

1. 本节要求

• 这节实验课要求掌握单细胞分析的标准主流程，从 10X 数据输入一直到聚类、注释、marker 提取。
• 要理解每一步是在干什么，不只是记函数名。
• 重点包括：
• 为什么要做细胞质控
• 为什么要筛高变基因
• PCA 后为什么还要选 PC 数量
• 聚类和 resolution 的关系
• 细胞注释为什么主要靠 marker

2. 考试掌握程度

• 老师在回放里说得很明确：实验课考试更可能考流程图，不太会让你现场完整编码。
• 这节课真正要掌握到的程度是：
• 能把单细胞分析流程用文字完整说出来
• 能画出标准流程图
• 能说明每一步的输入、处理、输出
• 知道关键参数大概是什么意思
• 对命令名的要求是“认识并会对应”，不是逐字背会所有代码。

3. 老师重点强调

• 单细胞分析里很多步骤都“默认要做”，尤其是：
• 质控
• 归一化
• 高变基因筛选
• PCA
• 聚类
• 注释
• 线粒体比例过滤是重点，老师多次提到：
• 普通情况下常按 < 10%
• 有些组织更严格时可能按 < 5%
• PC 选择方法要知道有哪几种，但老师也说实际中很多时候直接用前 20 或前 30 个 PC。
• 细胞注释时，主要看 marker 表达，不是完全靠自动工具。
• 课上提到过一些辅助工具，但老师明确说这类工具“不是主线”，不用把注意力放在工具细节上。

4. 与期末重点重合

• 高度重合，优先复习：
• 这是老师最后两节课明确点名的单细胞重点流程之一。
• 非常适合出应用题，尤其是“单细胞分析流程图 + 每步文字说明”。
• 同时也和期末要求中的这些项直接重合：
• 各类分析流程图
• 工具名与分析步骤对应关系
• 实验课代码思路
• 应用题的“文字描述 + 流程图”写法

说明

• 这份整理以老师板书为主，结合课件和第11周实验课回放做了必要补全。
• 老师这节课板书基本从 Read10X 一直写到 FindAllMarkers，主线非常完整。
• 但回放里老师也明确说了：这节实验课考试更重要的是“流程图 + 关键参数 + 每步作用”，不一定要求你把每个函数名一字不差写出来。

一、结合板书整理后的主干 R 代码

library(Seurat)
library(dplyr)
library(clustree)
library(ggplot2)
library(ggraph)

# 1. 读入 10X 数据并构建 Seurat 对象
pbmc.data <- Read10X(data.dir = "./filtered_feature_bc_matrix/")

pbmc <- CreateSeuratObject(
  counts = pbmc.data,
  project = "pbmc3k",
  min.cells = 3,
  min.features = 200
)

# 2.1 计算线粒体比例，先看过滤前 QC
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

p1 <- VlnPlot(
  pbmc,
  features = c("nFeature_RNA", "nCount_RNA", "percent.mt"),
  ncol = 3
)

p1

# 2.2 过滤低质量/异常细胞，再看过滤后 QC
pbmc <- subset(
  pbmc,
  subset = nFeature_RNA > 200 &
    nFeature_RNA < 3000 &
    percent.mt < 10
)

p2 <- VlnPlot(
  pbmc,
  features = c("nFeature_RNA", "nCount_RNA", "percent.mt"),
  ncol = 3
)

p2

# 3. 数据归一化 + 高变基因筛选
pbmc <- NormalizeData(
  pbmc,
  normalization.method = "LogNormalize",
  scale.factor = 10000
)

pbmc <- FindVariableFeatures(
  pbmc,
  selection.method = "vst",
  nfeatures = 2000
)

top10 <- head(VariableFeatures(pbmc), 10)
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)

plot1
plot2

# 4. 标准化 + PCA
all.genes <- rownames(pbmc)

pbmc <- ScaleData(
  pbmc,
  features = all.genes
)

pbmc <- RunPCA(
  pbmc,
  features = VariableFeatures(object = pbmc)
)

# 5.1 PCA 可视化与 ElbowPlot
VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")

DimHeatmap(
  pbmc,
  dims = 1:15,
  cells = 500,
  balanced = TRUE
)

ElbowPlot(pbmc, ndims = 30)

# 5.2 JackStraw 评估主成分
pbmc <- JackStraw(
  pbmc,
  num.replicate = 100,
  dims = 30
)

pbmc <- ScoreJackStraw(
  pbmc,
  dims = 1:30
)

JackStrawPlot(pbmc, dims = 1:30)

# 5.3 也可以粗略看累计方差贡献
pc_sd <- Stdev(pbmc, reduction = "pca")
cumulative_var <- cumsum(pc_sd^2) / sum(pc_sd^2)

cumulative_var[3]
cumulative_var[20]

# 实验里常直接选前 20 个 PC
dims_used <- 1:20

# 6. 聚类 + UMAP/TSNE
res.used <- seq(0.1, 1, by = 0.1)

pbmc <- pbmc %>%
  FindNeighbors(dims = dims_used) %>%
  FindClusters(resolution = res.used) %>%
  RunUMAP(dims = dims_used) %>%
  RunTSNE(dims = dims_used)

clas.tree.out <- clustree::clustree(pbmc) +
  theme(legend.position = "bottom") +
  scale_color_brewer(palette = "Set1") +
  scale_edge_color_continuous(
    low = "grey80",
    high = "red"
  )

print(clas.tree.out)

sel.clust <- "RNA_snn_res.0.3"
pbmc <- SetIdent(pbmc, value = sel.clust)
pbmc$seurat_clusters <- pbmc@meta.data[[sel.clust]]

# 6.3 聚类结果可视化
DimPlot(
  pbmc,
  reduction = "umap",
  label = TRUE,
  label.size = 5
)

DimPlot(
  pbmc,
  reduction = "tsne",
  label = TRUE,
  label.size = 5
)

# 7. 用 marker 基因染色做细胞注释
FeaturePlot(
  object = pbmc,
  features = c("PTPRC"),
  cols = c("gray", "blue"),
  max.cutoff = 2,
  min.cutoff = 0
)

# 8. 找各 cluster 的 marker 基因
pbmc.markers <- FindAllMarkers(
  pbmc,
  only.pos = TRUE
)

pbmc.markers %>%
  group_by(cluster) %>%
  dplyr::filter(avg_log2FC > 1)

top10 <- pbmc.markers %>%
  group_by(cluster) %>%
  slice_head(n = 10)

DoHeatmap(
  pbmc,
  features = top10$gene
) +
  NoLegend()

二、这份代码对应老师板书的哪几块

1. Read10X() 和 CreateSeuratObject()

对应第一张板书左上：

pbmc.data <- Read10X(data.dir = "./filtered_feature_bc_matrix/")

pbmc <- CreateSeuratObject(
  counts = pbmc.data,
  project = "pbmc3k",
  min.cells = 3,
  min.features = 200
)

这一步的意思就是：

• 读取 filtered_feature_bc_matrix 里的 10X 三个文件
• 构建 Seurat 对象
• 初步过滤：
• 少于 3 个细胞表达的基因去掉
• 少于 200 个基因的细胞去掉

课件也明确说了：

• Read10X() 可以直接读 .gz
• 不需要先手动解压

2. percent.mt 和两次 VlnPlot()

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

pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

p1 <- VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 3000 & percent.mt < 10)

p2 <- VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

这部分是单细胞分析里最重要的质控步骤之一：

• nFeature_RNA：每个细胞检测到的基因数
• nCount_RNA：每个细胞总 UMI 数
• percent.mt：线粒体比例

老师在回放里反复强调：

• 这一步非常关键
• 普通组织通常按 < 10% 过滤线粒体比例
• 高耗能组织可适当放宽
• 但如果放得太离谱，后面结论可能不可靠

3. NormalizeData() 和 FindVariableFeatures()

对应第二张板书右上和右下：

pbmc <- NormalizeData(
  pbmc,
  normalization.method = "LogNormalize",
  scale.factor = 10000
)

pbmc <- FindVariableFeatures(
  pbmc,
  selection.method = "vst",
  nfeatures = 2000
)

下面这几行也是老师板书里明确写了的：

top10 <- head(VariableFeatures(pbmc), 10)
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)

这部分就是：

• 先做归一化
• 再找高变基因
• 再把最突出的高变基因标出来看看

4. ScaleData() 和 RunPCA()

对应第一张板书右上和第三张板书中间：

all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))

这里的逻辑是：

• 对所有基因做标准化
• 但做 PCA 时主要用高变基因

5. 三种 PC 选择方法

对应第三、四张板书左侧：

5.1 ElbowPlot / PCA 可视化

VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")
DimHeatmap(pbmc, dims = 1:15, cells = 500, balanced = TRUE)
ElbowPlot(pbmc, ndims = 30)

5.2 JackStraw

pbmc <- JackStraw(pbmc, num.replicate = 100, dims = 30)
pbmc <- ScoreJackStraw(pbmc, dims = 1:30)
JackStrawPlot(pbmc, dims = 1:30)

5.3 累计方差

板书这里写得比较简略，意思是在看前若干个 PC 对总变异的累计贡献。

我把它整理成了更稳妥的写法：

pc_sd <- Stdev(pbmc, reduction = "pca")
cumulative_var <- cumsum(pc_sd^2) / sum(pc_sd^2)

回放里老师对这一段讲得很明确：

• 5.1、5.2、5.3 都是在帮助你选 PC 数量
• 但考试不一定要求你把三种图和命令都死背
• 实际上常直接取前 20 个左右 PC
• 经验上前 15 到 30 个 PC 都比较常见

6. 聚类、clustree、UMAP/TSNE

对应第四张板书：

res.used <- seq(0.1, 1, by = 0.1)

pbmc <- pbmc %>%
  FindNeighbors(dims = dims_used) %>%
  FindClusters(resolution = res.used) %>%
  RunUMAP(dims = dims_used) %>%
  RunTSNE(dims = dims_used)

再用：

clas.tree.out <- clustree::clustree(pbmc) + ...

去比较不同 resolution 下 cluster 的分裂情况。

最后板书明确选的是：

sel.clust <- "RNA_snn_res.0.3"

也就是课件里的：

• Resolution = 0.3

然后板书接着写：

pbmc <- SetIdent(pbmc, value = sel.clust)
pbmc$seurat_clusters <- pbmc@meta.data[[sel.clust]]

接下来再画：

DimPlot(..., reduction = "umap")
DimPlot(..., reduction = "tsne")

回放里老师也专门说了：

• UMAP 和 TSNE 都可以画
• 但更常用的是 UMAP
• 因为 UMAP 往往把 cluster 分得更开一些

7. FeaturePlot() 做 marker 染色

对应第五张板书右下：

FeaturePlot(
  object = pbmc,
  features = c("PTPRC"),
  cols = c("gray", "blue"),
  max.cutoff = 2,
  min.cutoff = 0
)

老师在回放里强调了：

• PTPRC 是免疫细胞常用 marker
• 如果要看别的细胞类型，就把 PTPRC 换成别的 marker 基因
• 细胞注释主要还是以 marker 表达为主

8. FindAllMarkers() 和热图

对应第五张板书左侧：

pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE)

再配合：

pbmc.markers %>%
  group_by(cluster) %>%
  dplyr::filter(avg_log2FC > 1)

以及：

top10 <- pbmc.markers %>%
  group_by(cluster) %>%
  slice_head(n = 10)

DoHeatmap(pbmc, features = top10$gene) + NoLegend()

这部分就是：

• 为每个 cluster 找 marker
• 取每类最典型的一批 marker
• 画热图辅助判断 cluster 身份

三、课件补充但板书没完全展开的代码

1. GetAssayData() 读取表达矩阵

课件专门提醒了 Seurat v5：

GetAssayData(pbmc, assay = "RNA", layer = "counts")
GetAssayData(pbmc, assay = "RNA", layer = "data")
GetAssayData(pbmc, assay = "RNA", layer = "scale.data")

也就是说：

• counts：原始 UMI
• data：归一化后矩阵
• scale.data：标准化后矩阵

这部分板书没展开写，但课件点得很明确。

2. SingleR 自动注释

课件后面补了一个 SingleR 方案：

library(SingleR)
library(celldex)

hpca.se <- celldex::HumanPrimaryCellAtlasData()
data.mat <- GetAssayData(pbmc, assay = "RNA", layer = "data")
clusters <- pbmc@meta.data[[sel.clust]]

pred.hesc <- SingleR(
  test = data.mat,
  ref = hpca.se,
  labels = hpca.se$label.main,
  clusters = clusters
)

celltype <- data.frame(
  ClusterID = rownames(pred.hesc),
  celltype = pred.hesc$labels,
  stringsAsFactors = FALSE
)

pbmc@meta.data$singleR <- celltype[
  match(clusters, celltype$ClusterID),
  "celltype"
]

DimPlot(
  pbmc,
  reduction = "umap",
  group.by = "singleR",
  label = TRUE,
  label.size = 4
)

不过回放里老师对这类“工具型注释”说得比较谨慎：

• 可以作为辅助
• 但不是最核心
• 主要还是以 marker 表达和生物学背景判断为主

3. 手工重命名 cluster

课件里还给了一个按 marker 和背景知识手工命名 cluster 的例子。

这一段只适合在你得到和课件一致的 cluster 结果时参考：

new.cluster.ids <- c(
  "Myeloid", "Epithelial", "Myeloid", "Myeloid", "TNK", "Myeloid",
  "Myeloid", "Myeloid", "Endothelial", "Myeloid", "Fibroblast",
  "Bcell", "Epithelial", "Epithelial"
)

names(new.cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, new.cluster.ids)
pbmc@meta.data$new.cluster.ids <- Idents(pbmc)

DimPlot(
  pbmc,
  reduction = "umap",
  label = TRUE,
  pt.size = 0.5,
  group.by = "new.cluster.ids"
)

这段不是板书主线，但它体现了课件后半部分“最终把 cluster 命名成细胞类型”的思路。

四、结合回放，老师口头最强调的是什么

1. 流程图比死背命令更重要

回放里老师明确说：

• 实验课考试更可能让你画流程图
• 不一定要求你完整编码
• 写流程图时不需要把每个命令都写出来

所以这节单细胞实验课，你最该会写的是这条主线：

10X数据输入
-> 创建Seurat对象
-> 计算percent.mt
-> 画QC图
-> 设置阈值过滤细胞
-> 归一化
-> 高变基因筛选
-> 标准化
-> PCA降维
-> 选择PC数量
-> FindNeighbors/FindClusters
-> UMAP/TSNE可视化
-> marker染色注释
-> FindAllMarkers
-> 热图展示

2. 这几个数字要记住

这节课最容易被拿来考的关键参数是：

• min.cells = 3
• min.features = 200
• 200 < nFeature_RNA < 3000
• percent.mt < 10
• 高变基因数 2000
• 常用 dims = 1:20
• 最终示例 resolution = 0.3

另外根据回放，老师还补充了两点：

• 线粒体比例在普通组织里常见按 < 10%，有时也会更严格到 < 5%
• PC 数量很多时候经验上在 15 到 30 之间，实际分析里直接取前 20 个非常常见

3. 单细胞注释主要看 marker，不是只靠工具

回放里老师反复强调：

• FeaturePlot 看 marker 表达是主线
• SingleR 这类工具只是辅助
• 最后判断 cluster 属于什么细胞类型，还是要结合 marker 和生物学背景

老师在回放里还专门提醒：

• 理想情况下，一个好 marker 应该对某一团细胞具有比较强的特异性
• 但实际图上常会看到别的 cluster 也有一些发蓝的细胞，这往往和聚类误差、数据噪声等有关
• 所以看 marker 时不要机械地“见蓝就认定”，而要结合整体表达模式判断

4. UMAP 比 TSNE 更常用

老师口头明确提到：

• 两种都可以画
• 但 UMAP 更常用
• 因为它通常把 cluster 分得更开，视觉上更清楚

五、你复习时最该盯住的代码逻辑

如果不要求逐行默写代码，这节课最值得背的是下面这一条主线：

Read10X
-> CreateSeuratObject
-> PercentageFeatureSet
-> VlnPlot看QC
-> subset按阈值过滤
-> NormalizeData
-> FindVariableFeatures
-> ScaleData
-> RunPCA
-> ElbowPlot / JackStraw / 累计方差选PC
-> FindNeighbors
-> FindClusters
-> clustree选resolution
-> RunUMAP / RunTSNE
-> FeaturePlot看marker
-> FindAllMarkers
-> DoHeatmap

其中最容易单独拿出来考的，是这几个点：

1. 为什么要先算 percent.mt
2. 为什么 VlnPlot 要在过滤前后各画一次
3. 为什么要先找高变基因再做 PCA
4. 三种选 PC 方法各自是在干什么
5. resolution 变大时 cluster 数为什么会增多
6. 为什么 marker 表达是细胞注释的主线
7. 为什么流程图里不需要写满函数命令，但必须把分析逻辑写全

六、如果考试要你手写伪代码，可以写成这样

输入：10X 单细胞表达矩阵

1. 读取 filtered_feature_bc_matrix
2. 构建 Seurat 对象
3. 计算线粒体基因比例
4. 绘制 QC 小提琴图
5. 根据 nFeature_RNA、percent.mt 等阈值过滤细胞
6. 对数据做归一化
7. 筛选高变基因
8. 对数据做标准化
9. 进行 PCA 降维
10. 结合 ElbowPlot、JackStraw 等方法选择 PC 数量
11. 计算邻居图并进行聚类
12. 用 clustree 选择合适 resolution
13. 用 UMAP 或 TSNE 可视化聚类结果
14. 用 marker 基因表达对 cluster 进行注释
15. 找出各 cluster 的 marker 基因
16. 用热图展示 top marker genes

七、一句提醒

这节实验课代码量看着大，但真正考试最重要的不是每个函数名，而是你要会解释：

• 输入是什么
• 过滤依据是什么
• PC 怎么选
• resolution 怎么选
• marker 怎么帮助注释细胞类型
• 最终输出了哪些图和哪些结果

再补一句回放里的原意：

• 如果考试考这节实验课，老师更可能让你“画出单细胞分析流程图，并把各种可能的分析写一下”
• 写流程图时不需要逐个写命令名
• 但必须知道每一步在做什么、前后是怎么衔接的