原文见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

hist(colSums(pbmc.data),
     breaks = 100, main = "Expression sum per cell",
     xlab = "Sum expression")

hist

##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)

VlnPlot

#查看基因数目, 线粒体基因占比与UMI数目的关系

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

FeatureScatter

#质控

• 筛选检测到基因数目超过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")

Total expression before normalisation

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")  

Total expression after normalisation

##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

VariableFeaturePlot

##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")

VizDimLoadings

两个主成分的展示

DimPlot(pbmc, reduction = "pca",split.by = 'ident')

DimPlot

DimHeatmap绘制基于单个主成分的热图，细胞和基因的排序都是基于他们的主成分分数。对于数据异质性的探索是很有帮助的，可以帮助用户选择用于下游分析的主成分维度。

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

DimHeatmap

展示多个主成分

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

DimHeatmap

##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)

JackStrawPlot

• ElbowPlot

ElbowPlot(pbmc)

ElbowPlot

启发式评估方法生成一个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")

DimPlot

#添加细胞类群xiba的标签
DimPlot(pbmc, reduction = "umap",label = TRUE)
LabelClusters(DimPlot(pbmc, reduction = "umap"),id = 'ident')

DimPlot

此时可以保存数据，方便下次直接导入数据修改图形。

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"))

VlnPlot

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

VlnPlot

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

FeaturePlot

DoHeatmap为指定的细胞和基因花表达热图。每个类群默认展示top 20标记基因。

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

DoHeatmap

#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()

DimPlot

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

参考：
Seurat - Guided Clustering Tutorial, Compiled: April 17, 2020
Getting started with Seurat