单细胞RNA测序（scRNA-seq）分析完整教程

目录

1. 简介
2. 环境准备
3. 数据加载
4. 质量控制
5. 数据标准化
6. 特征选择
7. 数据缩放
8. 降维分析
9. 细胞聚类
10. 聚类可视化
11. 标记基因鉴定
12. 细胞类型注释

简介

单细胞RNA测序（scRNA-seq）能够在单个细胞水平上测量基因表达，揭示细胞异质性和发现新的细胞类型。本教程使用Seurat包分析PBMC（外周血单核细胞）数据。

核心概念

• 表达矩阵：行是基因，列是细胞，值是基因在细胞中的表达量
• 质量控制：过滤低质量细胞和低表达基因
• 降维：从数万个基因降到几十个主成分，再到2D可视化
• 聚类：将相似的细胞分组，通常对应不同的细胞类型
• 标记基因：每种细胞类型特异性表达的基因

环境准备

安装和加载必要的包

1
2
3
4
5
6
7
8
9
10
11
12
13
14

# 检查并安装multtest包（用于多重检验）
if(!require(multtest))install.packages("multtest")

# Seurat是单细胞分析的核心包
if(!require(Seurat))install.packages("Seurat")

# dplyr用于数据操作和管道符%>%
if(!require(dplyr))install.packages("dplyr")

# patchwork用于组合多个图形
if(!require(patchwork))install.packages("patchwork")

# R.utils提供文件操作工具（如解压缩）
if(!require(R.utils))install.packages("R.utils")

为什么需要这些包？

• Seurat：提供完整的单细胞分析流程
• dplyr：简化数据操作，使代码更易读
• patchwork：方便地组合多个图形进行比较
• R.utils：处理压缩文件

清空工作环境

1

rm(list = ls())

作用：清除所有已存在的变量，确保分析环境干净，避免变量冲突。

数据加载

下载示例数据

1
2

download.file('https://cf.10xgenomics.com/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz','my.pbmc.gz')
untar(gunzip("my.pbmc.gz"))

数据说明：

• 10X Genomics的PBMC 3k数据集
• 包含约3000个外周血单核细胞
• 已经过初步过滤的数据

读取10X数据

1

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

函数功能：

• 读取10X Genomics格式的数据
• 返回稀疏矩阵（节省内存）
• 矩阵维度：基因 × 细胞

创建Seurat对象

1

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

参数解释：

• counts：原始的基因表达计数矩阵
• project：项目名称，用于标识数据来源
• min.cells = 3：基因至少在3个细胞中表达才保留
• min.features = 200：细胞至少表达200个基因才保留

生物学意义：

• 过滤掉极低表达的基因（可能是技术噪音）
• 过滤掉基因数过少的细胞（可能是空液滴或死细胞）

查看数据基本信息

1
2
3
4
5
6
7
8

pbmc
# 输出：13714 features across 2700 samples within 1 assay

ncol(pbmc)  # 细胞数
# 输出：2700

ncol(pbmc.data)  # 原始数据的细胞数
# 输出：2700

导出表达矩阵（可选）

1
2
3

# Seurat v5的正确写法（v4使用@counts）
lalala <- as.data.frame(GetAssayData(pbmc, slot = "counts"))
write.table(lalala,'mycount.txt',sep = '\t')

注意：Seurat v5改变了数据访问方式，使用$代替@

质量控制

计算线粒体基因百分比

1

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

代码解析：

• pbmc[["percent.mt"]]：在metadata中创建新列
• PercentageFeatureSet()：计算特定基因集的表达百分比
• pattern = "^MT-"：正则表达式，匹配以"MT-"开头的基因

生物学原理：

• 健康细胞：线粒体基因表达比例低（<5-10%）
• 受损/死亡细胞：线粒体基因表达比例高
• 原因：细胞膜破损时，胞质mRNA流失，但线粒体mRNA相对保留

注意事项：

• 小鼠数据使用pattern = "^mt-"（小写）
• 不同组织的阈值可能不同

查看细胞元数据

1

head(pbmc@meta.data, 5)

显示每个细胞的：

• orig.ident：样本来源
• nCount_RNA：总RNA分子数
• nFeature_RNA：检测到的基因数
• percent.mt：线粒体基因百分比

质控指标可视化

1

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

小提琴图解读：

• nFeature_RNA：基因数分布，识别双细胞（异常高）或低质量细胞（异常低）
• nCount_RNA：总分子数分布，反映测序深度
• percent.mt：线粒体百分比，高值提示细胞损伤

特征相关性分析

1
2
3

plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
CombinePlots(plots = list(plot1, plot2))

散点图意义：

• 图1：分子数vs线粒体百分比，识别异常细胞
• 图2：分子数vs基因数，应呈正相关，偏离者可能是双细胞

过滤低质量细胞

1

pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)

过滤标准：

• nFeature_RNA > 200：去除基因数过少的细胞（空液滴）
• nFeature_RNA < 2500：去除基因数过多的细胞（双细胞）
• percent.mt < 5：去除线粒体含量高的细胞（死细胞）

如何确定阈值？

• 观察分布图，找到明显的异常值
• 参考文献中类似数据的阈值
• 根据具体实验调整

数据标准化

1

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

为什么需要标准化？

不同细胞的RNA总量差异原因：

1. 细胞大小不同
2. 细胞周期阶段不同
3. 技术偏差（捕获效率）

LogNormalize方法

步骤1：CPM标准化

1

标准化值 = (基因计数 / 细胞总计数) × scale.factor

步骤2：对数转换

1

最终值 = log(标准化值 + 1)

计算示例

假设两个细胞：

标准化后表达值相同，消除了细胞大小的影响！

参数说明：

• scale.factor = 10000：相当于CPM（counts per ten thousand）
• log1p转换：稳定方差，使数据更接近正态分布

特征选择

寻找高变基因

1

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

什么是高变基因？

高变基因在不同细胞间表达差异大，它们：

• 能够区分不同细胞类型
• 包含最多的生物学信息
• 不是"管家基因"（在所有细胞中稳定表达）

VST方法原理

1. 计算均值-方差关系

   • 高表达基因往往有高方差（技术噪音）
   • 建立模型预测期望方差

2. 标准化方差

   1	标准化方差 = (观察方差 - 预期方差) / 预期方差

3. 选择top 2000基因

查看高变基因

1
2

# 前10个高变基因
top10 <- head(VariableFeatures(pbmc), 10)

典型的高变基因包括：

• 免疫球蛋白基因（IGLC2, IGKC）
• HLA基因（抗原呈递）
• 细胞类型标记基因

可视化高变基因

1
2
3
4
5
6
7
8

# 基础散点图
plot1 <- VariableFeaturePlot(pbmc)

# 添加基因标签
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)

# 并排显示
plot1 + plot2

图形解读：

• X轴：平均表达量
• Y轴：标准化方差
• 红点：被选中的2000个高变基因
• 标签：前10个最高变基因

数据缩放

1

pbmc <- ScaleData(pbmc, features = rownames(pbmc))

什么是ScaleData？

对每个基因进行Z-score标准化：

1

缩放值 = (表达值 - 基因均值) / 基因标准差

结果：每个基因均值=0，标准差=1

为什么要缩放？

1. 消除量级差异

   • 基因A：表达范围1-10
   • 基因B：表达范围1000-10000
   • 不缩放，基因B会主导分析

2. PCA的要求

   • PCA对方差敏感
   • 需要所有变量在同一尺度

3. 可视化需要

   • 热图等需要标准化数据

注意：features = rownames(pbmc)缩放所有基因，耗时较长。通常只缩放高变基因即可。

降维分析

主成分分析（PCA）

1

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

PCA原理

将高维数据（2000个基因）投射到低维空间（50个主成分）：

• 每个PC是基因的线性组合
• PC1解释最多的变异
• PC之间相互正交

生物学解释

每个PC通常代表特定生物学过程：

• PC1：常区分髓系vs淋巴系
• PC2：常区分T细胞vs B细胞
• PC3-5：细胞状态、细胞周期等

查看PCA结果

1
2
3
4
5
6
7
8
9
10
11

# 简要信息
print(pbmc[["pca"]], dims = 1:5, nfeatures = 5)

# 基因贡献可视化
VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")

# 细胞在PC空间的分布
DimPlot(pbmc, reduction = "pca")

# 主成分热图
DimHeatmap(pbmc, dims = 1:15, cells = 500, balanced = TRUE)

确定使用多少个PC

方法1：ElbowPlot（推荐）

1

ElbowPlot(pbmc)

解读：

• 显示每个PC解释的变异量
• 找"肘部"：曲线明显变平的点
• 通常选择10-30个PC

方法2：JackStraw分析（可选）

1
2
3

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

原理：统计检验每个PC的显著性
缺点：计算耗时

细胞聚类

构建细胞邻近图

1

pbmc <- FindNeighbors(pbmc, dims = 1:10)

功能：

• 基于前10个PC计算细胞间距离
• 构建KNN（K近邻）图
• 构建SNN（共享近邻）图

聚类分析

1

pbmc <- FindClusters(pbmc, resolution = 0.5)

算法：Louvain算法（社区检测）
参数：

• resolution：控制聚类粒度
  • 低值（0.1-0.5）：少而大的聚类
  • 高值（0.5-2）：多而小的聚类

输出信息解读：

1

Number of communities: 9

表示识别出9个细胞群

聚类可视化

UMAP降维

1
2

pbmc <- RunUMAP(pbmc, dims = 1:10)
DimPlot(pbmc, reduction = "umap")

UMAP特点：

• 非线性降维
• 保留局部和全局结构
• 更好地分离不同细胞群

t-SNE降维

1
2

pbmc <- RunTSNE(pbmc, dims = 1:10)
DimPlot(pbmc, reduction = "tsne")

t-SNE特点：

• 非线性降维
• 强调局部结构
• 聚类边界更清晰

UMAP vs t-SNE：

• UMAP更快，保留更多全局结构
• t-SNE聚类分离更明显
• 两者互补，可同时使用

标记基因鉴定

特定聚类的标记基因

1

cluster5.markers <- FindMarkers(pbmc, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25)

比较聚类5与聚类0、3的差异基因

所有聚类的标记基因

1

pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)

参数说明：

• only.pos = TRUE：只找高表达基因
• min.pct = 0.25：基因至少在25%的细胞中表达
• logfc.threshold = 0.25：log2倍数变化阈值

查看每个聚类的top标记基因

1
2
3

pbmc.markers %>% 
  group_by(cluster) %>% 
  top_n(n = 2, wt = avg_log2FC)

结果解读：

• p_val：原始p值
• avg_log2FC：平均log2倍数变化
• pct.1：聚类内表达比例
• pct.2：其他聚类表达比例
• p_val_adj：校正后p值

标记基因可视化

小提琴图

1

VlnPlot(pbmc, features = c("NKG7", "PF4"), slot = "counts", log = TRUE)

显示基因在各聚类的表达分布

特征图

1

FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ"))

在UMAP上显示基因表达

热图

1
2
3
4
5

top10 <- pbmc.markers %>% 
  group_by(cluster) %>% 
  top_n(n = 10, wt = avg_log2FC)
  
DoHeatmap(pbmc, features = top10$gene) + NoLegend()

显示各聚类的标记基因表达模式

细胞类型注释

基于标记基因注释

1
2
3
4

new.cluster.ids <- c("Naive CD4 T", "CD14+ Mono", "Memory CD4 T", "B", "CD8 T", 
                     "FCGR3A+ Mono", "NK", "DC", "Platelet")
names(new.cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, new.cluster.ids)

常见PBMC细胞类型标记

可视化注释结果

1

DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()

显示带细胞类型标签的UMAP图

保存结果

1

saveRDS(pbmc, 'pbmc.rds')

保存完整的Seurat对象，包含所有分析结果

重新加载

1
2

pbmc <- readRDS('pbmc.rds')
DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()

总结

单细胞RNA测序分析的标准流程：

1. 数据准备：加载数据，创建Seurat对象
2. 质量控制：过滤低质量细胞
3. 数据处理：标准化、特征选择、缩放
4. 降维分析：PCA、UMAP/t-SNE
5. 细胞聚类：识别细胞群
6. 生物学解释：标记基因、细胞类型注释

重要提示

1. 参数选择：没有通用的最佳参数，需根据数据特点调整
2. 生物学验证：结合已知标记基因验证结果
3. 迭代优化：可能需要多次调整参数获得最佳结果
4. 版本兼容：注意Seurat版本差异（v4 vs v5）

扩展学习

• 批次效应校正：多样本整合分析
• 轨迹推断：细胞发育轨迹分析
• 细胞通讯：配受体相互作用分析
• 空间转录组：结合空间信息的分析

参考资源

• Seurat官方教程
• 单细胞最佳实践
• CellMarker数据库

__END__