Seurat官方教程 | pbmc3k_tutorial最新版
#### 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
数据下载下来后进行解压，主要包含以下三个文件：
### 2.加载包和读取数据
```R
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长什么样？
```R
# 首先看看三个基因的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来存储两种方式的大小对比
```R
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
sparse为.即稀疏矩阵的大小，29905192 bytes
两者比值为23.7 倍，即使用系数矩阵来存储单细胞水平的基因表达值非常节省空间。
### 3.单细胞数据分析预处理
预处理主要包括基于QC指标的细胞和基因过滤，数据标准化和归一化，高变基因选择。
#### 3.1首先是QC来筛选高质量的细胞
一般筛选条件有三个：
- 1.每个细胞中检测到的唯一基因数
- 低质量的细胞或者空的droplet液滴通常含有很少的基因
- Cell doublets双胞体或多胞体含有很高的异常的gene counts
- 2.每个细胞中检测到的分子总数
- 3.线粒体基因含量比例
- 低质量或者死亡细胞含有很高的线粒体基因
- 使用[PercentageFeatureSet()]( https://satijalab.org/seurat/reference/PercentageFeatureSet.html) 计算一个特征的比例
- MT-开头的基因认为是线粒体基因
```R
# 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
过滤前三个指标可视化
```R
# Visualize QC metrics as a violin plot
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
```
每个特征之间的相关性
```R
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
```
过滤
```R
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
```
#### 3.2数据标准化
过滤了不想要的cells之后，下一步就是标准化。标准化方法为：LogNormalize。标准化的数据存在pbmc[["RNA"]]@data中。
```R
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。
```R
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
```
#### 3.4数据归一化
在数据降维之前，我们使用[ScaleData()]( https://satijalab.org/seurat/reference/ScaleData.html) 进行数据归一化。归一化的功能：
- 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.
```R
# 默认只对2千个高变基因进行数据归一化
pbmc <- ScaleData(pbmc)
# 也可以使用全部的基因进行归一化，但是耗时会比较久
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
```
使用2000个高边基因和所有的基因进行归一化对后续的PCA和聚类分析并不会有太大影响，有影响的地方是后面绘制热图的时候 `DoHeatmap()`，因为绘制热图需要输入归一化后的基因表达值。
#### 3.5 进行线性降维
接下来，我们对归一化后的数据进行PCA分析。默认的，只是用前面决定的高变基因进行PCA分析，也可以使用features参数设置用户自己选择的基因进行PCA分析。
```R
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")
```
```R
DimPlot(pbmc, reduction = "pca")
```
此外，[DimHeatmap()]( https://satijalab.org/seurat/reference/DimHeatmap.html) 允许我们探索数据中的最初的异质性，然后决定下游使用多少个PC进行分析。
```R
DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE)
```
```R
DimHeatmap(pbmc, dims = 1:15, cells = 500, balanced = TRUE)
```
#### 3.6 确定数据的维度
Seurat对细胞进行聚类主要基于他们的PCA打分，每一个PC代表一个综合特征，它综合了数据中相关基因表达的一些信息。前几的主成分代表了一个数据集的稳定的综合的信息。那么我们应该选择多少个PC进行分析呢？
我们使用了一个样本随机扰动的方法：the JackStraw procedure
我们随机选择一个数据的subset（默认为数据得1%）进行PCA分析，返回PCA打分，然后重复这个过程，形成一个打分的零分布。
然后我们根据这个零分布，计算PC的显著性，根据P值，选择显著的PC。
```R
# 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)
```
还可以根据PC的贡献度进行排序绘制一个散点图
```R
ElbowPlot(pbmc)
```
根据上图，我们可以看到在第10个PC处有一个拐点。
#### 3.7 聚类
我们使用k最近邻法来对细胞进行聚类（K-nearest neighbor (KNN) graph），具有相似基因表达模式的细胞之间绘制边，然后尝试将这张图划分为高度相互关联的“准派系”或“群体”。
PhenoGrap中，我们基于PCA空间的欧氏距离，构建一个KNN图，然后根据两个细胞的局部邻居的共享重叠部分(Jaccard相似性)来优化它们之间的边权值。这一步使用[FindNeighbors()]( https://satijalab.org/seurat/reference/FindNeighbors.html) 函数，使用前面决定的10个PCs作为输入。
为了对细胞聚类，我们接下来使用模块化优化技术如the Louvain algorithm (default) or SLM [[SLM, Blondel *et al*., Journal of Statistical Mechanics\]]( http://dx.doi.org/10.1088/1742-5468/2008/10/P10008) ，迭代地将细胞分组。这一步主要使用[FindClusters()]( https://satijalab.org/seurat/reference/FindClusters.html) ，**其中有个函数resolution会随着值变大，分类数变多。对于单细胞数据3k个细胞，resolution值为0.4-1.2就可以返回一个比较好的聚类结果。对于较大的数据集，最佳分辨率通常会增加。**
```R
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个类。
```R
# Look at cluster IDs of the first 5 cells
head(Idents(pbmc), 5)
```
#### 3.8 非线性降维(UMAP/tSNE)
Seurat提供几个非线性降维方法，如tSNE和UMAP，来可视化和探索这些数据集。这些算法的目标是学习数据的底层流形，以便将相似的细胞放在低维空间中。
```R
# 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)
```
保存结果下次使用
```R
dir.create("output")
saveRDS(pbmc, file = "output/pbmc_tutorial.rds")
```
#### 3.9 差异表达分析
这一步主要使用FindMarkers和[FindAllMarkers]( https://satijalab.org/seurat/reference/FindAllMarkers.html) 。
```R
# 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()]( https://satijalab.org/seurat/reference/FeaturePlot.html) ，RidgePlot()`, `CellScatter()，[DotPlot()]( https://satijalab.org/seurat/reference/DotPlot.html)
```R
VlnPlot(pbmc, features = c("MS4A1", "CD79A"))
```
```R
# you can plot raw counts as well
VlnPlot(pbmc, features = c("NKG7", "PF4"), slot = "counts", log = TRUE)
```
```R
FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"))
```
这里作者展示的marker基因都是根据他的背景知识来的，他事先就知道这些个基因是哪些细胞类型里面特异表达的基因，可以用来鉴定细胞类型。
如果是我们，可能就比较困难，需要看大量的文献，或者借助已知marker的数据库。
绘制每个cluster的差异top20的基因的热图
```R
top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
DoHeatmap(pbmc, features = top10$gene) + NoLegend()
```
#### 3.10 细胞类型注释
这个地方也是整个单细胞数据分析里面比较难的部分，有很多种方法可以进行细胞类型注释。这里作者采用了已知的canonical markers进行注释。
```R
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()
```
保存最终的结果
```R
saveRDS(pbmc, file = "output/pbmc3k_final.rds")
```
### 总结
综合以上内容，去掉可视化部分，标准分析步骤有
```R
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")
```
