原文见Seurat - Guided Clustering Tutorial, Compiled: April 17, 2020
#1 Seurat安装
install.packages("Seurat")
#2 数据下载
Peripheral Blood Mononuclear Cells (PBMC) 是10X Genomics dataset page提供的一个数据，包含2700个单细胞，出自Illumina NextSeq 500平台。
PBMCs是来自健康供体具有相对少量RNA（around 1pg RNA/cell）的原代细胞。在Illumina NextSeq 500平台，检测到2700个单细胞，每个细胞获得69000 reads。
tar -xvf pbmc3k_filtered_gene_bc_matrices.tar
文件夹下包含3个文件
barcodes.tsv
genes.tsv
matrix.mtx
matrix.mtx：matrix.mtx 是 MatrixMarket格式文件；更多内容见：http://math.nist.gov/MatrixMarket/formats.html
文件中储存非零值;
注释使用%标记；
第一行包含文件中总行数，总列数，总的记录数
每行中提供记录的所处的行号和列号，已经记录的内容
head filtered_gene_bc_matrices/hg19/matrix.mtx
%%MatrixMarket matrix coordinate real general
%
32738 2700 2286884
32709 1 4
32707 1 1
32706 1 10
32704 1 1
32703 1 5
32702 1 6
32700 1 10
#3 数据导入
##3.1 Read10X()函数可以自动读入和解析数据。
library(dplyr)
library(Seurat)
library(patchwork)
#读取PBMC数据
pbmc.data <- Read10X(data.dir = "../data/pbmc3k/filtered_gene_bc_matrices/hg19/")
#查看数据
dim(pbmc.data)
# 32738 2700
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
summary(colSums(pbmc.data))
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 548 1758 2197 2367 2763 15844
#查看每个细胞有多少基因被检测到
at_least_one <- apply(pbmc.data, 2, function(x) sum(x>0))
hist(at_least_one, breaks = 100,
main = "Distribution of detected genes",
xlab = "Genes with at least one tag")
hist(colSums(pbmc.data),
breaks = 100, main = "Expression sum per cell",
xlab = "Sum expression")
##3.2 使用pbmc数据初始化Seurat对象
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)
head(pbmc$RNA@data[,1:5])
6 x 5 sparse Matrix of class "dgCMatrix"
AAACATACAACCAC-1 AAACATTGAGCTAC-1 AAACATTGATCAGC-1 AAACCGTGCTTCCG-1 AAACCGTGTATGCG-1
AL627309.1 . . . . .
AP006222.2 . . . . .
RP11-206L10.2 . . . . .
RP11-206L10.9 . . . . .
LINC00115 . . . . .
NOC2L . . . . .
#4 数据预处理
这部分是基于数据质控方法，标准化和检测到的变化基因对数据进行筛选。
##4.1 对细胞的质控
可以参考文章：Classification of low quality cells from single-cell RNA-seq data
- 单个细胞中检测到单个基因的数目
- 低质量的细胞以及空泡油滴中一般检测到很少的基因
- 包含多个细胞的油滴会检测到异常多的基因
- 类似，在一个单细胞中的基因总数
- 检测到线粒体的基因数目百分比
- 低质量/垂死细胞通常会有线粒体污染
- 使用PercentageFeatureSet函数函数可以计算线粒体QC，计算线粒体基因所占百分比
- 基因名以MT- 开始的基因定义为线粒体基因
#线粒体基因占比计算
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
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
#画图查看基因数目, UMI数目, 线粒体基因占比
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
#查看基因数目, 线粒体基因占比与UMI数目的关系
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2
#质控
- 筛选检测到基因数目超过2500或低于200的细胞
- 单个细胞中线粒体基因数目占比超过>5%
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
##4.2 数据标准化
默认使用数据标准化方法是LogNormalize, 每个细胞总的表达量都标准化到10000，然后log取对数；结果存放于pbmc[["RNA"]]@data。
#标准化前，每个细胞总的表达量
hist(colSums(pbmc$RNA@data),
breaks = 100,
main = "Total expression before normalisation",
xlab = "Sum of expression")
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
#标准化后，每个细胞总的表达量
hist(colSums(pbmc$RNA@data),
breaks = 100,
main = "Total expression after normalisation",
xlab = "Sum of expression")
##4.3 变化基因鉴定
鉴定在细胞间表达高度变化的基因，后续研究需要集中于这部分基因。Seurat内置的FindVariableFeatures()函数，首先计算每一个基因的均值和方差，并且直接模拟其关系。默认返回2000个基因。
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
# 10个表达变化最为剧烈的基因
top10 <- head(VariableFeatures(pbmc), 10) #head(pbmc$RNA@var.features,10)
# "PPBP" "LYZ" "S100A9" "IGLL5" "GNLY" "FTL" "PF4" "FTH1" "GNG11" "S100A8"
# 画出表达变化的基因，从而观察其分布
plot1 <- VariableFeaturePlot(pbmc)
# 画出表达变化的基因，标记前10个基因
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1
plot2
##4.4 数据缩放
线性转换缩放数据，ScaleData()函数可以实现此功能。
最终每个基因均值为0，方差为1。
结果存放于pbmc[["RNA"]]@scale.data。
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
设置参数features是因为ScaleData默认处理前面鉴定的差异基因。这一步怎么做都不会影响到后续pca和聚类，但是会影响做热图。
移除影响方差的因素
pbmc <- ScaleData(pbmc, vars.to.regress = "percent.mt")
#5 线性降维分析
##5.1 PCA
对缩放后的数据进行PCA分析，默认使用前面鉴定表达变化大的基因。使用features参数可以重新定义数据集。
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))
VizDimReduction, DimPlot, 和 DimHeatmap可以从基因或细胞角度可视化pca结果
#查看对每个主成分影响比较大的基因集
print(pbmc[["pca"]], dims = 1:5, nfeatures = 5)
## PC_ 1
## Positive: CST3, TYROBP, LST1, AIF1, FTL
## Negative: MALAT1, LTB, IL32, IL7R, CD2
## PC_ 2
## Positive: CD79A, MS4A1, TCL1A, HLA-DQA1, HLA-DQB1
## Negative: NKG7, PRF1, CST7, GZMB, GZMA
## PC_ 3
## Positive: HLA-DQA1, CD79A, CD79B, HLA-DQB1, HLA-DPB1
## Negative: PPBP, PF4, SDPR, SPARC, GNG11
## PC_ 4
## Positive: HLA-DQA1, CD79B, CD79A, MS4A1, HLA-DQB1
## Negative: VIM, IL7R, S100A6, IL32, S100A8
## PC_ 5
## Positive: GZMB, NKG7, S100A8, FGFBP2, GNLY
## Negative: LTB, IL7R, CKB, VIM, MS4A7
#可视化对每个主成分影响比较大的基因集
VizDimLoadings(pbmc, dims = 1:2, reduction = "pca")
两个主成分的展示
DimPlot(pbmc, reduction = "pca",split.by = 'ident')
DimHeatmap绘制基于单个主成分的热图，细胞和基因的排序都是基于他们的主成分分数。对于数据异质性的探索是很有帮助的，可以帮助用户选择用于下游分析的主成分维度。
DimHeatmap(pbmc, dims = 1, cells = 500, balanced = TRUE)
展示多个主成分
DimHeatmap(pbmc, dims = 1:15, cells = 500, balanced = TRUE)
##5.1 数据维度
为了避免单个基因影响，Seurat聚类细胞时使用pca结果。首先需要确定的是使用多少个主成分用于后续分析。常用有两种方法，一种是基于零分布的统计检验方法，这种方法耗时且可能不会返回明确结果。另一种是主成分分析常用的启发式评估。
- JackStraw()
在JackStraw()函数中, 使用基于零分布的置换检验方法。随机抽取一部分基因(默认1%)然后进行pca分析得到pca分数，将这部分基因的pca分数与先前计算的pca分数进行比较得到显著性p-Value，。根据主成分（pc）所包含基因的p-value进行判断选择主成分。最终的结果是每个基因与每个主成分的关联的p-Value。保留下来的主成分是那些富集小的p-Value基因的主成分。
处理大数据时会花费大量时间；ElbowPlot()内置了一些其它的方法可以减少运行时间。
pbmc <- JackStraw(pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(pbmc, dims = 1:20)
JackStrawPlot()函数提供可视化方法，用于比较每一个主成分的p-value的分布，虚线是均匀分布；显著的主成分富集有小p-Value基因，实线位于虚线左上方。下图表明保留10个pca主成分用于后续分析是比较合理的。
JackStrawPlot(pbmc, dims = 1:15)
- ElbowPlot
ElbowPlot(pbmc)
启发式评估方法生成一个Elbow plot图。在图中展示了每个主成分对数据方差的解释情况（百分比表示），并进行排序。根据自己需要选择主成分，图中发现第9个主成分是一个拐点，后续的主成分(PC)变化都不大了。
注意：鉴别数据的真实维度不是件容易的事情；除了上面两种方法，Serat官当文档还建议将主成分(数据异质性的相关来源有关)与GSEA分析相结合。Dendritic cell 和 NK aficionados可能识别的基因与主成分 12 和 13相关，定义了罕见的免疫亚群 (i.e. MZB1 is a marker for plasmacytoid DCs)。如果不是事先知道的情况下，很难发现这些问题。
Serat官当文档因此鼓励用户使用不同数量的PC（10、15，甚至50）重复下游分析。其实也将观察到的，结果通常没有显著差异。因此，在选择此参数时，可以尽量选大一点的维度，维度太小的话对结果会产生不好的影响。
#6 细胞聚类
Seurat v3应用基于图形的聚类方法，例如KNN方法。具有相似基因表达模式的细胞之间绘制边缘，然后将他们划分为一个内联群体。
在PhenoGraph中，首先基于pca维度中（先前计算的pca数据）计算欧式距离（the euclidean distance），然后根据两个细胞在局部的重合情况（Jaccard 相似系数）优化两个细胞之间的边缘权值。此步骤内置于FindNeighbors函数，输入时先前确定的pc数据。
为了聚类细胞，接下来应用模块化优化技术迭代将细胞聚集在一起。（the Louvain algorithm (default) or SLM [SLM, Blondel et al., Journal of Statistical Mechanics]），FindClusters函数实现这一功能，其中需要注意resolution参数，该参数设置下游聚类分析的“granularity”，更大的resolution会导致跟多的细胞类群。3000左右的细胞量，设置resolution为0.4-1.2是比较合适的。细胞数据集越大，需要更大的resolution参数, 会获得更多的细胞聚类。
查看细胞属于那个类群可以使用函数Idents。
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
#查看细胞属于那个类群
head(Idents(pbmc), 5)
AAACATACAACCAC-1 AAACATTGAGCTAC-1 AAACATTGATCAGC-1 AAACCGTGCTTCCG-1 AAACCGTGTATGCG-1
0 3 2 5 6
Levels: 0 1 2 3 4 5 6 7 8
#7 非线性降维分析
Seurat提供了一些非线性降维度分析的方法，如tSNE和UMAP，以可视化和探索这些数据集；这一步使用的数据建议与聚类分析使用的pc维度一致。
# install UMAP： reticulate::py_install(packages ='umap-learn')
pbmc <- RunUMAP(pbmc, dims = 1:10)
#画图展示
DimPlot(pbmc, reduction = "umap")
#添加细胞类群xiba的标签
DimPlot(pbmc, reduction = "umap",label = TRUE)
LabelClusters(DimPlot(pbmc, reduction = "umap"),id = 'ident')
此时可以保存数据，方便下次直接导入数据修改图形。
saveRDS(pbmc, file = "../output/pbmc_tutorial.rds")
#8 寻找差异表达基因 (cluster biomarkers)
Seurat可以通过差异表达分析寻找不同细胞类群的标记基因。FindMarkers函数可以进行此操作，但是默认寻找单个类群(参数ident.1)与其他所有类群阳性和阴性标记基因。FindAllMarkers函数会自动寻找每个类群和其他每个类群之间的标记基因。
min.pct参数:设定在两个细胞群中任何一个被检测到的百分比，通过此设定不检测很少表达基因来缩短程序运行时间。默认0.1
thresh.test参数：设定在两个细胞群中基因差异表达量。可以设置为0 ，程序运行时间会更长。
max.cells.per.ident参数：每个类群细胞抽样设置；也可以缩短程序运行时间。
# find all markers of cluster 1
cluster1.markers <- FindMarkers(pbmc, ident.1 = 1, min.pct = 0.25)
head(cluster1.markers, n = 5)
p_val avg_logFC pct.1 pct.2 p_val_adj
IL32 1.894810e-92 0.8373872 0.948 0.464 2.598542e-88
LTB 7.953303e-89 0.8921170 0.981 0.642 1.090716e-84
CD3D 1.655937e-70 0.6436286 0.919 0.431 2.270951e-66
IL7R 3.688893e-68 0.8147082 0.747 0.325 5.058947e-64
LDHB 2.292819e-67 0.6253110 0.950 0.613 3.144372e-63
# 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)
p_val avg_logFC pct.1 pct.2 p_val_adj
FCGR3A 7.583625e-209 2.963144 0.975 0.037 1.040018e-204
IFITM3 2.500844e-199 2.698187 0.975 0.046 3.429657e-195
CFD 1.763722e-195 2.362381 0.938 0.037 2.418768e-191
CD68 4.612171e-192 2.087366 0.926 0.036 6.325132e-188
RP11-290F20.3 1.846215e-188 1.886288 0.840 0.016 2.531900e-184
# 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_logFC)
p_val avg_logFC pct.1 pct.2 p_val_adj cluster gene
<dbl> <dbl> <dbl> <dbl> <dbl> <fct> <chr>
1 1.96e-107 0.730 0.901 0.594 2.69e-103 0 LDHB
2 1.61e- 82 0.922 0.436 0.11 2.20e- 78 0 CCR7
3 7.95e- 89 0.892 0.981 0.642 1.09e- 84 1 LTB
4 1.85e- 60 0.859 0.422 0.11 2.54e- 56 1 AQP3
5 0. 3.86 0.996 0.215 0. 2 S100A9
6 0. 3.80 0.975 0.121 0. 2 S100A8
7 0. 2.99 0.936 0.041 0. 3 CD79A
8 9.48e-271 2.49 0.622 0.022 1.30e-266 3 TCL1A
9 2.96e-189 2.12 0.985 0.24 4.06e-185 4 CCL5
10 2.57e-158 2.05 0.587 0.059 3.52e-154 4 GZMK
11 3.51e-184 2.30 0.975 0.134 4.82e-180 5 FCGR3A
12 2.03e-125 2.14 1 0.315 2.78e-121 5 LST1
13 7.95e-269 3.35 0.961 0.068 1.09e-264 6 GZMB
14 3.13e-191 3.69 0.961 0.131 4.30e-187 6 GNLY
15 1.48e-220 2.68 0.812 0.011 2.03e-216 7 FCER1A
16 1.67e- 21 1.99 1 0.513 2.28e- 17 7 HLA-DPB1
17 7.73e-200 5.02 1 0.01 1.06e-195 8 PF4
18 3.68e-110 5.94 1 0.024 5.05e-106 8 PPBP
Seurat可以通过参数test.use设定检验差异表达的方法(详情见 DE vignett)。
cluster1.markers <- FindMarkers(pbmc, ident.1 = 0, logfc.threshold = 0.25, test.use = "roc", only.pos = TRUE)
head(cluster1.markers, n = 5)
Seurat有多种方法可视化标记基因的方法
- VlnPlot: 基于细胞类群的基因表达概率分布
- FeaturePlot：在tSNE 或 PCA图中画出基因表达情况
- RidgePlot，CellScatter，DotPlot
VlnPlot(pbmc, features = c("MS4A1", "CD79A"))
# you can plot raw counts as well
VlnPlot(pbmc, features = c("NKG7", "PF4"), slot = "counts", log = TRUE)
FeaturePlot(pbmc, features = c("MS4A1", "GNLY", "CD3E", "CD14", "FCER1A", "FCGR3A", "LYZ", "PPBP",
"CD8A"))
DoHeatmap为指定的细胞和基因花表达热图。每个类群默认展示top 20标记基因。
top10 <- pbmc.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
DoHeatmap(pbmc, features = top10$gene) + NoLegend()
#9 Assigning cell type identity to clusters
根据已知细胞类型与基因标记的对应关系，可以为细胞类群找到对应的细胞类型。
| Cluster ID | Markers | Cell Type |
|---|---|---|
| 0 | IL7R, CCR7 | Naive CD4+ T |
| 1 | IL7R, S100A4 | Memory CD4+ |
| 2 | CD14, LYZ | CD14+ Mono |
| 3 | MS4A1 | B |
| 4 | CD8A | CD8+ T |
| 5 | FCGR3A, MS4A7 | FCGR3A+ Mono |
| 6 | GNLY, NKG7 | NK |
| 7 | FCER1A, CST3 | DC |
| 8 | PPBP | Platelet |
new.cluster.ids <- c("Naive CD4 T", "Memory CD4 T", "CD14+ Mono", "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()
saveRDS(pbmc, file = "../output/pbmc3k_final.rds")
参考：
Seurat - Guided Clustering Tutorial, Compiled: April 17, 2020
Getting started with Seurat
