Compiled: March 18, 2021

来源：https://satijalab.org/seurat/articles/pbmc3k_tutorial.html

这个教程可以说是每个单细胞学习的人的最初的入门教程，里面的东西可以反复琢磨，标准分析基本上就这么多内容。教程中提到了大量的资源介绍，参数含义解读，甚至高级部分的简洁版算法介绍，结果解读，可视化方法，以及数据变换后使用什么函数来提取对应的变换后的数据。

1.数据介绍

此次分析所有数据来源于Peripheral Blood Mononuclear Cells (PBMC)，使用10X Genomics， Illumina NextSeq 500进行测序，共获得2700个单细胞。可以在这里进行下载：https://cf.10xgenomics.com/samples/cell/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz

数据下载下来后进行解压，主要包含以下三个文件：

image-20210404002351434.png

2.加载包和读取数据

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

# Load the PBMC dataset, 读取数据主要使用Read10X这个函数
pbmc.data <- Read10X(data.dir = "data/filtered_gene_bc_matrices/hg19")

# Initialize the Seurat object with the raw (non-normalized data).
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
pbmc

An object of class Seurat 
13714 features across 2700 samples within 1 assay 
Active assay: RNA (13714 features, 0 variable features)

总共得到2700个细胞和13714个基因

count matrix长什么样？

# 首先看看三个基因的count值
pbmc.data[c("CD3D", "TCL1A", "MS4A1"), 1:30]

## 3 x 30 sparse Matrix of class "dgCMatrix"
##                                                                    
## CD3D  4 . 10 . . 1 2 3 1 . . 2 7 1 . . 1 3 . 2  3 . . . . . 3 4 1 5
## TCL1A . .  . . . . . . 1 . . . . . . . . . . .  . 1 . . . . . . . .
## MS4A1 . 6  . . . . . . 1 1 1 . . . . . . . . . 36 1 2 . . 2 . . . .

其中的.表示0，即no molecules detected。当然，这个地方还有另外一种含义就是这个基因是真的没有表达。

由于单细胞测序数据中大多数的值都为0，因此，seurat使用一个稀疏矩阵来保存测序得到的count matrix，这样有利于数据存储空间的节省。

我们来看看使用稀疏矩阵和使用0来存储两种方式的大小对比

dense.size <- object.size(as.matrix(pbmc.data))
dense.size

sparse.size <- object.size(pbmc.data)
sparse.size

dense.size/sparse.size
23.7 bytes

dense为转换为0后的matrix存储大小，709591472 bytes

image-20210404003451736.png

sparse为.即稀疏矩阵的大小，29905192 bytes

两者比值为23.7 倍，即使用系数矩阵来存储单细胞水平的基因表达值非常节省空间。

3.单细胞数据分析预处理

预处理主要包括基于QC指标的细胞和基因过滤，数据标准化和归一化，高变基因选择。

3.1首先是QC来筛选高质量的细胞

一般筛选条件有三个：

• 1.每个细胞中检测到的唯一基因数
  • 低质量的细胞或者空的droplet液滴通常含有很少的基因
  • Cell doublets双胞体或多胞体含有很高的异常的gene counts
• 2.每个细胞中检测到的分子总数
• 3.线粒体基因含量比例
  • 低质量或者死亡细胞含有很高的线粒体基因
  • 使用PercentageFeatureSet()计算一个特征的比例
  • MT-开头的基因认为是线粒体基因

# The [[ operator can add columns to object metadata. This is a great place to stash QC stats
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

# 查看QC指标
# Show QC metrics for the first 5 cells
head(pbmc@meta.data, 5)

                 orig.ident nCount_RNA nFeature_RNA percent.mt
AAACATACAACCAC-1     pbmc3k       2419          779  3.0177759
AAACATTGAGCTAC-1     pbmc3k       4903         1352  3.7935958
AAACATTGATCAGC-1     pbmc3k       3147         1129  0.8897363
AAACCGTGCTTCCG-1     pbmc3k       2639          960  1.7430845
AAACCGTGTATGCG-1     pbmc3k        980          521  1.2244898

我们将使用以下标准进行基因过滤：

• We filter cells that have unique feature counts over 2,500 or less than 200
• We filter cells that have >5% mitochondrial counts

过滤前三个指标可视化

# Visualize QC metrics as a violin plot
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

image-20210404004816633.png

每个特征之间的相关性

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

image-20210404005036149.png

过滤

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

3.2数据标准化

过滤了不想要的cells之后，下一步就是标准化。标准化方法为：LogNormalize。标准化的数据存在pbmc[["RNA"]]@data中。

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

Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

3.3识别高变基因

高标基因：在一些细胞中高表达，在另外的细胞中低表达。

使用均值与方差之间的关系，来挑选高变基因，默认返回前2000个高变基因进入下游分析，如PCA。

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

# Identify the 10 most highly variable genes
top10 <- head(VariableFeatures(pbmc), 10)

# plot variable features with and without labels
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1 + plot2

image-20210404130613422.png

3.4数据归一化

在数据降维之前，我们使用ScaleData()进行数据归一化。归一化的功能：

• Shifts the expression of each gene, so that the mean expression across cells is 0
• Scales the expression of each gene, so that the variance across cells is 1
• The results of this are stored in pbmc[["RNA"]]@scale.data

总的来说就是使得每一个基因在所有cell中的表达均值为0，方差为1.

# 默认只对2千个高变基因进行数据归一化
pbmc <- ScaleData(pbmc)

# 也可以使用全部的基因进行归一化，但是耗时会比较久
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)

使用2000个高边基因和所有的基因进行归一化对后续的PCA和聚类分析并不会有太大影响，有影响的地方是后面绘制热图的时候 DoHeatmap()，因为绘制热图需要输入归一化后的基因表达值。

3.5 进行线性降维

接下来，我们对归一化后的数据进行PCA分析。默认的，只是用前面决定的高变基因进行PCA分析，也可以使用features参数设置用户自己选择的基因进行PCA分析。

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

# Examine and visualize PCA results a few different ways
print(pbmc[["pca"]], dims = 1:5, nfeatures = 5)

VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")

image-20210404133815898.png

DimPlot(pbmc, reduction = "pca")

image-20210404133917621.png

此外，DimHeatmap()允许我们探索数据中的最初的异质性，然后决定下游使用多少个PC进行分析。

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

image-20210404134234949.png

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

image-20210404134337643.png

3.6 确定数据的维度

Seurat对细胞进行聚类主要基于他们的PCA打分，每一个PC代表一个综合特征，它综合了数据中相关基因表达的一些信息。前几的主成分代表了一个数据集的稳定的综合的信息。那么我们应该选择多少个PC进行分析呢？

我们使用了一个样本随机扰动的方法：the JackStraw procedure

我们随机选择一个数据的subset（默认为数据得1%）进行PCA分析，返回PCA打分，然后重复这个过程，形成一个打分的零分布。

然后我们根据这个零分布，计算PC的显著性，根据P值，选择显著的PC。

# NOTE: This process can take a long time for big datasets, comment out for expediency. 
# More approximate techniques such as those implemented in ElbowPlot() can be used to 
# reduce computation time
pbmc <- JackStraw(pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(pbmc, dims = 1:20)

# 对PC打分进行可视化
JackStrawPlot(pbmc, dims = 1:15)

image-20210404140452215.png

还可以根据PC的贡献度进行排序绘制一个散点图

ElbowPlot(pbmc)

image-20210404141358293.png

根据上图，我们可以看到在第10个PC处有一个拐点。

3.7 聚类

我们使用k最近邻法来对细胞进行聚类（K-nearest neighbor (KNN) graph），具有相似基因表达模式的细胞之间绘制边，然后尝试将这张图划分为高度相互关联的“准派系”或“群体”。

PhenoGrap中，我们基于PCA空间的欧氏距离，构建一个KNN图，然后根据两个细胞的局部邻居的共享重叠部分(Jaccard相似性)来优化它们之间的边权值。这一步使用FindNeighbors()函数，使用前面决定的10个PCs作为输入。

为了对细胞聚类，我们接下来使用模块化优化技术如the Louvain algorithm (default) or SLM [SLM, Blondel et al., Journal of Statistical Mechanics]，迭代地将细胞分组。这一步主要使用FindClusters()，其中有个函数resolution会随着值变大，分类数变多。对于单细胞数据3k个细胞，resolution值为0.4-1.2就可以返回一个比较好的聚类结果。对于较大的数据集，最佳分辨率通常会增加。

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

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 2638
Number of edges: 96033

Running Louvain algorithm...
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Maximum modularity in 10 random starts: 0.8720
Number of communities: 9
Elapsed time: 0 seconds

这里聚成了9个类。

# Look at cluster IDs of the first 5 cells
head(Idents(pbmc), 5)

3.8 非线性降维(UMAP/tSNE)

Seurat提供几个非线性降维方法，如tSNE和UMAP，来可视化和探索这些数据集。这些算法的目标是学习数据的底层流形，以便将相似的细胞放在低维空间中。

# If you haven't installed UMAP, you can do so via reticulate::py_install(packages =
# 'umap-learn')
pbmc <- RunUMAP(pbmc, dims = 1:10)

# note that you can set `label = TRUE` or use the LabelClusters function to help label
# individual clusters
DimPlot(pbmc, reduction = "umap", label=T)

image-20210404144600590.png

保存结果下次使用

dir.create("output")
saveRDS(pbmc, file = "output/pbmc_tutorial.rds")

3.9 差异表达分析

这一步主要使用FindMarkers和FindAllMarkers。

# find all markers of cluster 1
cluster1.markers <- FindMarkers(pbmc, ident.1 = 2, min.pct = 0.25)
head(cluster1.markers, n = 5)

# find all markers distinguishing cluster 5 from clusters 0 and 3
cluster5.markers <- FindMarkers(pbmc, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25)
head(cluster5.markers, n = 5)

# find markers for every cluster compared to all remaining cells, report only the positive ones
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
pbmc.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_log2FC)

cluster1.markers <- FindMarkers(pbmc, ident.1 = 0, logfc.threshold = 0.25, test.use = "roc", only.pos = TRUE)

可视化差异表达基因：VlnPlot，FeaturePlot()，RidgePlot(),CellScatter()，DotPlot()

VlnPlot(pbmc, features = c("MS4A1", "CD79A"))

image-20210404151249117.png

# you can plot raw counts as well
VlnPlot(pbmc, features = c("NKG7", "PF4"), slot = "counts", log = TRUE)

image-20210404151629784.png

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

image-20210404151855546.png

这里作者展示的marker基因都是根据他的背景知识来的，他事先就知道这些个基因是哪些细胞类型里面特异表达的基因，可以用来鉴定细胞类型。

如果是我们，可能就比较困难，需要看大量的文献，或者借助已知marker的数据库。

绘制么个cluster的差异top20的基因的热图

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

image-20210404152831073.png

3.10 细胞类型注释

这个地方也是整个单细胞数据分析里面比较难的部分，有很多种方法可以进行细胞类型注释。这里作者采用了已知的canonical markers进行注释。

image-20210404153007569.png

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)
DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()

image-20210404153231048.png

保存最终的结果

saveRDS(pbmc, file = "output/pbmc3k_final.rds")

总结

综合以上内容，去掉可视化部分，标准分析步骤有

pbmc.counts <- Read10X(data.dir = "pbmc3k/filtered_gene_bc_matrices/hg19/")
pbmc <- CreateSeuratObject(counts = pbmc.counts)
pbmc <- NormalizeData(object = pbmc)
pbmc <- FindVariableFeatures(object = pbmc)
pbmc <- ScaleData(object = pbmc)
pbmc <- RunPCA(object = pbmc)
pbmc <- FindNeighbors(object = pbmc)
pbmc <- FindClusters(object = pbmc)
pbmc <- RunTSNE(object = pbmc)
pbmc <- RunUMAP(object = pbmc)
DimPlot(object = pbmc, reduction = "tsne")