单细胞转录组学的记录（part1） | 鲁珀的部落格 lupo's blog

今天第一次把单细胞转录组分析（ssRNA-seq）啃下来了。

在这之前我第一次了解到这东西还是在免疫学的课上，当老师讲授组学的时候，我第一次就意识到了这东西的前景无比巨大，从而着迷上了这东西。正巧最近申请到了硕士的机会，同时博士申请也要准备上马了，借着这个东风于是就开始学了起来。

今天也就简短的记录下使用公开数据PBMC（Peripheral Blood Mononuclear Cells）进行单细胞转录组分析的过程。数据可以从这里获取

本篇文章的所有思路和代码都来自于[Seurat - Guided Clustering Tutorial](Analysis, visualization, and integration of Visium HD spatial datasets with Seurat • Seurat)欢迎大家访问原始网站查看更多的内容。

整个分析的过程使用R和R studio完成，并且使用了Seurat。

首先是对环境的准备，在R studio里面需要加载这三个包

1
2
3

library(dplyr)
library(Seurat)
library(patchwork)

在准备好分析的环境以后，将我们下载的数据文件解压，并且使用Read10X()命令将数据库导入到R studio内

# 加载 PBMC 数据集
pbmc.data <- Read10X(data.dir = "数据文件夹路径")
# 初始化 Seurat 对象
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
#查看 Seurat 对象
pbmc

在完成这三个步骤以后，输入pbmc即可看到下面一系列的输出，这里告诉了我们Seurat 对象的内容

An object of class Seurat 
13714 features across 2638 samples within 1 assay 
Active assay: RNA (13714 features, 2000 variable features)
 3 layers present: counts, data, scale.data
 2 dimensional reductions calculated: pca, umap

随后我们就可以开始进行QC的工作

QC的目的是筛选出低质量细胞

一般的筛选指标包括这三个

每个细胞检测到的基因数量（nFeature_RNA）：
- 过低：可能表示空液滴
- 过高：可能表示双细胞或多细胞。
每个细胞的总 UMI 计数（nCount_RNA）：
- 过低：可能表示细胞质量差。
- 过高：可能表示双细胞或多细胞。
线粒体基因比例（percent.mt）：
- 过高：可能表示细胞死亡或破裂。

通过这下面的代码，对数据集进行QC

#计算线粒体基因比例
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
# 使用提琴图对 QC 指标可视化
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

随后我们要对数据进行过滤

常用的过滤指标有:200 < nFeature_RNA < 2500, 500 < nCount_RNA < 20000,percent.mt < 5

但是这些指标并不唯一，还是要根据提琴图进行判断，在本例中将使用以上的过滤指标对数据进行QC

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

这样，我们就完成了数据的QC。

接下来我们对数据进行归一化，这一步的目的是测序深度差异和技术噪声

1 2	#使用 Log 归一化方法对数据进行归一化 pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)

当完成了这些步骤以后，对数据的前期处理就完成了。我们继续来了解如何进行特征选择

所谓的特征选择，就是筛选具有分辨意义的特征基因（即，它们在某些细胞中高表达，而在其他细胞中低表达），例如CD4基因，我们利用这些基因来筛选特定的细胞群体。我们可以通过以下代码进行

#筛选高变基因
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)

# 选择前10个表达量差异最大的基因
top10 <- head(VariableFeatures(pbmc), 10)

# 绘制高变基因的可视化图
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1 + plot2

在完成了这些步骤以后，我们对数据的前期处理就已经完成了，在第二部分我们继续来了解数据缩放，PCA分析，数据降维。