Seurat包------标准流程
单细胞学习入门教程，里面的东西可以反复琢磨，标准分析基本上就这么多内容。教程中提到了大量的资源介绍，参数含义解读，甚至高级部分的简洁版算法介绍，结果解读，可视化方法，以及数据变换后使用什么函数来提取对应的变换后的数据。
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.png
2、安装Seurat包
安装Seurat包
install.packages('Seurat')
library('Seurat')
packageVersion("Seurat")
相关包下载
install.packages(c('dplyr','patchwork'))
3、加载包和读取数据
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)
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
as.matrix(pbmc.data)[c("CD3D", "TCL1A", "MS4A1"), 1:30]
AAACATACAACCAC AAACATTGAGCTAC AAACATTGATCAGC AAACCGTGCTTCCG
CD3D 4 0 10 0
TCL1A 0 0 0 0
MS4A1 0 6 0 0
AAACCGTGTATGCG AAACGCACTGGTAC AAACGCTGACCAGT AAACGCTGGTTCTT
CD3D 0 1 2 3
TCL1A 0 0 0 0
MS4A1 0 0 0 0
AAACGCTGTAGCCA AAACGCTGTTTCTG AAACTTGAAAAACG AAACTTGATCCAGA
CD3D 1 0 0 2
TCL1A 1 0 0 0
MS4A1 1 1 1 0
AAAGAGACGAGATA AAAGAGACGCGAGA AAAGAGACGGACTT AAAGAGACGGCATT
CD3D 7 1 0 0
TCL1A 0 0 0 0
MS4A1 0 0 0 0
4、单细胞数据分析预处理
预处理主要包括基于QC指标的细胞和基因过滤，数据标准化和归一化，高变基因选择。
4.1首先是QC来筛选高质量的细胞
一般筛选条件有三个：
1.每个细胞中检测到的唯一基因数
低质量的细胞或者空的droplet液滴通常含有很少的基因
Cell doublets双胞体或多胞体含有很高的异常的gene counts
2.线粒体基因含量比例
低质量或者死亡细胞含有很高的线粒体基因
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
过滤前三个指标可视化
# Visualize QC metrics as a violin plot
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
image.png
每个特征之间的相关性
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
image.png
过滤
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
4.2数据标准化
过滤了不想要的cells之后，下一步就是标准化。标准化方法为：LogNormalize。标准化的数据存在pbmc[["RNA"]]@data中。
归一化和标准化的区别
首先去除样本/细胞效应：因为不同样本或者细胞的测序数据量不一样，所以同样的一个基因在不同细胞，哪怕看到的表达量是一样的，但是背后细胞整体测序数据量的差异其实反而说明了这个基因在不同细胞表达量其实是有差异的。
然后去除基因效应，这个主要是在绘制热图的时候会需要使用，因为个别基因表达量超级高，在热图里面一枝独秀，实际上我们并不会关心不同基因的表达量高低，我们仅仅是想看指定基因在不同细胞的高低而已，这样的话，就把该基因的表达量在不同细胞的数值，进行z-score这样的转换。这样处理后的表达矩阵，就可以进行后续的降维聚类分群。
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
Performing log-normalization
0% 10 20 30 40 50 60 70 80 90 100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
4.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.png
4.4识别高变基因
在数据降维之前，我们使用[ScaleData()进行数据归一化。归一化的功能：使得每一个基因在所有cell中的表达均值为0，方差为1。
# 默认只对2千个高变基因进行数据归一化
pbmc <- ScaleData(pbmc)
# 也可以使用全部的基因进行归一化，但是耗时会比较久
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
使用2000个高变基因和所有的基因进行归一化对后续的PCA和聚类分析并不会有太大影响，有影响的地方是后面绘制热图的时候 DoHeatmap()，因为绘制热图需要输入归一化后的基因表达值。
4.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.png
DimPlot(pbmc, reduction = "pca")
image.png
此外，[DimHeatmap()]允许我们探索数据中的最初的异质性，然后决定下游使用多少个PC进行分析。
DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE)
image.png
DimHeatmap(pbmc, dims = 1:15, cells = 500, balanced = TRUE)
image.png
4.6 确定数据的维度
这里我们需要选择出主成分的数目，用于后续细胞分类。这里定义的“维度”并不代表细胞类型的数目，而是对细胞分类时需要用到的一个参数。
Seurat对细胞进行聚类主要基于他们的PCA打分，每一个PC代表一个综合特征，它综合了数据中相关基因表达的一些信息。前几的主成分代表了一个数据集的稳定的综合的信息。那么我们应该选择多少个PC进行分析呢？
我们使用了一个样本随机扰动的方法：the JackStraw procedure
我们随机选择一个数据的subset（默认为数据得1%）进行PCA分析，返回PCA打分，然后重复这个过程，形成一个打分的零分布。
# 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.png
还可以根据PC的贡献度进行排序绘制一个散点图,碎石图
ElbowPlot(pbmc)
image.png
根据上图，我们可以看到在第10个PC处有一个拐点。
4.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)
4.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.png
保存结果下次使用
dir.create("output")
saveRDS(pbmc, file = "output/pbmc_tutorial.rds")
4.9 差异表达分析
这一步主要使用FindMarkers和[FindAllMarkers]
# find all markers of cluster 1
cluster1.markers <- FindMarkers(pbmc, ident.1 = 1, 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 = 1, logfc.threshold = 0.25, test.use = "roc", only.pos = TRUE)
可视化差异表达基因：VlnPlot、RidgePlot()、CellScatter()，
VlnPlot(pbmc, features = c("MS4A1", "CD79A"))
image.png
image-20210404151249117.png
# you can plot raw counts as well
VlnPlot(pbmc, features = c("NKG7", "PF4"), slot = "counts", log = TRUE)
image.png
FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP", "CD8A"))
image.png
image-20210404151855546.png
这里作者展示的marker基因都是根据他的背景知识来的，他事先就知道这些个基因是哪些细胞类型里面特异表达的基因，可以用来鉴定细胞类型。
如果是我们，可能就比较困难，需要看大量的文献，或者借助已知marker的数据库。
绘制9个cluster的差异top20的基因的热图
top20 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 20, wt = avg_log2FC)
DoHeatmap(pbmc, features = top20$gene) + NoLegend()
image.png
3.10 细胞类型注释
这个地方也是整个单细胞数据分析里面比较难的部分，有很多种方法可以进行细胞类型注释。这里作者采用了已知的canonical markers进行注释。 NovoCM辅助细胞定义+SingleR细胞定义。
image.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.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")
