Seurat standard pipeline(10核心流程)

创建Seurat对象 Read10X CreateSeuratObjecet
质控 PercentageFeatureSubject subset
标准化Normalization NormalizeData
高变基因选择 FindVariableFeatures
数据缩放 ScaleData
线性降维 RunPCA
维数选择 FindNeighbors
细胞聚类 FindClusters
非线性降维（UMAP/tSNE）RunUMAP DimPlot
鉴定差异表达特征（cluster markers）
细胞注释 CellMarker---FeaturePlot---RenameIdents

下载示例数据：

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

一、创建 Seurat 分析对象

1.1 Read10X() 函数可以自动读入和解析数据

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

#读取PBMC数据
pbmc.data <- Read10X(data.dir = "./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")

1.2 pbmc 数据初始化 Seurat 对象

使用 CreateSeuratObject 创建对象。

sce<- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
sce
# An object of class Seurat 
# 13714 features across 2700 samples within 1 assay 
# Active assay: RNA (13714 features, 0 variable features)

head(sce$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                        .                .                .                .                .

#过滤检测 
min.features = 200：一个细胞最少要检测到200个基因，
min.cells = 3：一个基因最少得在4个细胞中表达

其他参数解释

• counts：未标准化的数据，如原始计数或TPMs
• project = "CreateSeuratObject"：设置Seurat对象的项目名称
• assay = "RNA"：与初始输入数据对应的分析名称
• names.field = 1：对于每个cell的初始标识类，从cell的名称中选择此字段。例如，如果cell在输入矩阵中被命名为BARCODE_CLUSTER_CELLTYPE，则设置名称。字段设置为3以将初始标识设置为 #CELLTYPE。
• names.delim = "_"：对于每个cell的初始标识类，从cell的列名中选择此分隔符。例如，如果cell命名为bar - cluster - celltype，则将此设置为“-”，以便将cell名称分离到其组成部分中，以选择相关字段。
• meta.data = NULL：要添加到Seurat对象的其他单元级元数据。应该是data.frame，其中行是单元格名称，列 #是附加的元数据字段。

这里还涉及到的知识点是R data中的S3类和S4类。
list一般情况下被认为是S3类，S4是指使用slots储存数据的格式。

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

#1个数据集，包含2700个细胞，13714个基因。

.在矩阵中的值表示0(未检测到分子)。

由于scRNA-seq矩阵中的大多数值为0,因此Seurat在任何可能的情况下都使用稀疏矩阵表示。这为Drop-seq/inDrop/10x数据节省了大量内存和速度。

所谓稀疏矩阵，也就是在矩阵中，若数值为0的元素数目远多于非0元素的数目，并且非0元素分布没有规律时，则称该矩阵为稀疏矩阵；与之相反，若非0元素数目占大多数时，则称该矩阵为稠密矩阵。

二、数据预处理

这部分是基于数据质控方法，标准化和检测到的变化基因对数据进行筛选。

1. 对细胞的质控 subset()

可以参考文章：Classification of low quality cells from single-cell RNA-seq data

常用的质控指标：
（1）每个细胞在检测到的特异基因数

• 低质量的细胞以及空泡油滴中一般检测到很少的基因
• 包含多个细胞的油滴会检测到异常多的基因

（2）每个细胞检测到的分子总数（与基因密切相关）
（3）每个细胞的线粒体基因比例。低质量/濒死细胞常表现出广泛的线粒体污染

使用PercentageFeatureSet()函数计算线粒体QC指标，使用所有以MT-开头的基因作为一组线粒体基因

（1）线粒体基因占比计算

PercentageFeatureSet（）参数意义:
pattern = "^MT-"为正则表达式写法，意为匹配"MT-"开头的基因，即线粒体基因。区分MT大小写，此处如果小写则代表另一个物种。

为什么将线粒体基因作为一个过滤条件？“因为线粒体是参与细胞凋亡启动和执行的主要细胞器之一。线粒体基因高表达意为着样品质量差，导致大量细胞凋亡或裂解。以及特定样本的生物学，例如肿瘤活检，可能由于代谢活动和/或坏死而增加线粒体基因表达。对于线粒体基因比例过高的细胞，会干扰细胞分群，在分析过程中为了避免这些细胞的影响，通常会根据实际情况设计一个阈值进行细胞的过滤

###向 sce 新增一列 percent.mt 数据
sce[['percent.mt']] <- PercentageFeatureSet(sce,pattern='^MT-')

#展示前5个细胞的QC指标
head(sce@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

QC指标储存在哪？
每个细胞基因数和总分子数在建立seurat对象时就已经自动计算好了。

可以使用sce@meta.data$percent.mt或者 sce['percent.mt']查看每个细胞的线粒体比例

（2）画图查看基因数目, UMI数目, 线粒体基因占比

# 使用小提琴图可视化QC指标
VlnPlot(sce, 
        features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), 
        ncol = 3)

nFeature_RNA： 代表每个细胞测到的基因数目。
nCount_RNA ：代表每个细胞测到所有基因的表达量之和。
percent.mt：代表测到的线粒体基因的比例

（3）查看基因数目, 线粒体基因占比与UMI数目的关系
FeatureScatter通常用于可视化 feature-feature 相关性，

#nCount_RNA 与 percent.mt的相关性
plot1 <- FeatureScatter(sce, 
                        feature1 = "nCount_RNA", 
                        feature2 = "percent.mt")

#nCount_RNA 与 nFeature_RNA的相关性
plot2 <- FeatureScatter(sce, 
                        feature1 = "nCount_RNA", 
                        feature2 = "nFeature_RNA")

plot1 + plot2 #合并两图

拿 nCounts_RNA 与 nFeature_RNA 的散点图（scatter）来说：

• 每个点代表一个细胞，斜率代表随着count的增加gene的增加程度。
• count和gene一般呈现线性关系，斜率越大也就是较少的count就可以检出较多的gene，说明这批细胞基因的丰度偏低（普遍贫穷）；反之，斜率较小，说明这批细胞基因丰度较高（少数富有）。
• 有的时候不是一条可拟合的线，或者是两条可拟合的线，也反映出细胞的异质性。
• 总之，他就是一个散点图，描述的是两个变量的关系**。

过滤线粒体基因表达比例过高的细胞，和一些极值细胞（可以根据小提琴图判断，查看两端离群值）。
（4）质控 subset()

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

nCount_RNA:每个细胞的UMI数量
nFeature_RNA:每个细胞的基因数
选取 2500 > nFeature_RNA >200 和percent.mt < 5的数据

三、数据标准化 NormalizeData()

因为在测序之前会对抓取到的RNA进行PCR扩增，所以需要考虑文库深度对测序的影响。例：细胞A中总reads数10，基因a的reads数为2；细胞B中总reads数1000，基因a的reads数为20。并不能说基因a在细胞B中表达量大于A，故需要进行标准化。所以需要对上一步得到的稀疏矩阵进行Normalize。

Normalize方式：每个细胞每个基因的特征计数除以该细胞（一列）的特征总计数，再乘以scale.factor(默认10,000)，然后对结果进行log1p对数转换（log1p=log(n+1)）。标准化的数值存储在sce[["RNA"]]@data中。

scale.factor = 10,000 的原因是

进行log转换的原因是：
以bulk-rna-seq为例，用FPKM或者TPM绘制热图的时候，因为数值变化的范围太过巨大，都需要进行一个log转换，让数据压缩在一个区间内。其次，也是最重要的改变数据分布：测序数值本身不符合正态分布，log转换能让数据趋近于正态分布，以便于后续进一步分析

默认使用数据标准化方法是LogNormalize, 每个细胞总的表达量都标准化到10000，然后log取对数；标准化后的数据保存在结果存放于sce[["RNA"]]@data。，原始数保存在sce[["RNA"]]@counts中

标准化前，每个细胞总的表达量

hist(colSums(sce$RNA@data),
     breaks = 100,
     main = "Total expression before normalisation",
     xlab = "Sum of expression")

标准化

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

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

若所有调用的参数都是默认值，则可省去

sce <- NormalizeData(sce)

标准化后，每个细胞总的表达量

hist(colSums(sce$RNA@data),
     breaks = 100,
     main = "Total expression after normalisation",
     xlab = "Sum of expression")  

标准化前和标准化后对比

四、 鉴定高变基因 FindVariableFeatures()

接下来，鉴定在细胞间表达高度变化的基因 (即，它们在某些细胞中高表达，而在其他细胞中低表达)。后续研究需要集中于这部分基因这些基因有助于突出单细胞数据集中的生物信号。

用FindVariableFeatures()函数实现，首先计算每一个基因的均值和方差，并且直接模拟其关系。默认返回2000个features 。这些将用于下游分析，如PCA。

高变基因方法选择vst
1.使用loss(局部加权回归)拟合平滑曲线模型
2.获取模型计算的值作为y = var.exp值
3.var.standarlized = get variance after feature standardization: (每个基因-mena)/sd后取var()，注意sd=sqrt(var.exp)
4.按照var.standalized降序排列，取前n(比如2000)个基因作为高变基因

sce <- FindVariableFeatures(sce, selection.method = "vst", nfeatures = 2000)
Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|


# 查看最高变的10个基因
top10 <- head(VariableFeatures(sce), 10)
#head(sce$RNA@var.features,10)
# "PPBP"   "LYZ"    "S100A9" "IGLL5"  "GNLY"   "FTL"    "PF4"    "FTH1"   "GNG11"  "S100A8"


# 画出不带标签或带标签基因点图
plot1 <- VariableFeaturePlot(sce)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1 + plot2

五、数据缩放（归一化）

线性变换(“缩放”)，是在PCA降维之前的一个标准预处理步骤。最终每个基因均值为0，方差为1。

ScaleData()函数功能:

• 转换每个基因的表达值，使每个细胞的平均表达为0
• 转换每个基因的表达值，使细胞间的方差为1

此步骤在下游分析中具有相同的权重，因此高表达的基因不会占主导地位
结果存储在sce[["RNA"]]@scale.data中

all.genes <- rownames(sce)
sce<- ScaleData(sce, features = all.genes)
Centering and scaling data matrix
  |================================================================| 100%

做了 Normalize 还做 scale 的原因是：
normalize改变的是数据的分布，scale改变的是数据的范围。
normalize是将偏态分布转换成趋近于正态分布，sclae是将数据的分布限定在一个范围内。
R语言的Z score计算是通过[scale()]函数求得，Seurat包中ScaleData()函数也基本参照了scale()函数的功能。
scale方法中的两个参数：center和scale

• center和scale默认为真，即T或者TRUE
• center为真表示数据中心化：数据集中的各项数据减去数据集的均值。如有数据集1, 2, 3, 6, 3，其均值为3 那么中心化之后的数据集为1-3,2-3,3-3,6-3,3-3，即：-2,-1,0,3,0 *scale为真表示数据标准化：指中心化之后的数据在除以数据集的标准差，即数据集中的各项数据减去数据集的均值再除以数据集的标准差。如有数据集1, 2, 3, 6, 3，其均值为3，其标准差为1.87 那么标准化之后的数据集为(1-3)/1.87,(2-3)/1.87,(3-3)/1.87,(6-3)/1.87,(3-3)/1.87，即：-1.069,-0.535,0,1.604,0

Z score的概念是指原始数据距离均值有多少个标准差。当以标准差为单位进行测量时，Z score 衡量的是一个数值偏离总体均值以上或以下多少个标准差。如果原始数值高于均值，则 Z score得分为正，如果低于均值，则Z score为负。Z score其实是标准正态分布（Standard Normal Distribution），即平均值μ=0，标准差 σ=1 的正态分布。SND标准正态分布的直方图如下所示：

六、线性降维

1、PCA

接下来，对缩放的数据执行PCA。默认情况下，只使用前面确定的变量特性作为输入，但是如果想选择不同的子集，可以使用features参数来定义(选择高变基因)。

> sce <- RunPCA(sce,features = VariableFeatures(object = sce))
PC_ 1
Positive:  CST3, TYROBP, LST1, AIF1, FTL, FTH1, LYZ, FCN1, S100A9, TYMP
           FCER1G, CFD, LGALS1, S100A8, CTSS, LGALS2, SERPINA1, IFITM3, SPI1, CFP
           PSAP, IFI30, SAT1, COTL1, S100A11, NPC2, GRN, LGALS3, GSTP1, PYCARD
Negative:  MALAT1, LTB, IL32, IL7R, CD2, B2M, ACAP1, CD27, STK17A, CTSW
           CD247, GIMAP5, AQP3, CCL5, SELL, TRAF3IP3, GZMA, MAL, CST7, ITM2A
           MYC, GIMAP7, HOPX, BEX2, LDLRAP1, GZMK, ETS1, ZAP70, TNFAIP8, RIC3
PC_ 2
Positive:  CD79A, MS4A1, TCL1A, HLA-DQA1, HLA-DQB1, HLA-DRA, LINC00926, CD79B, HLA-DRB1, CD74
           HLA-DMA, HLA-DPB1, HLA-DQA2, CD37, HLA-DRB5, HLA-DMB, HLA-DPA1, FCRLA, HVCN1, LTB
           BLNK, P2RX5, IGLL5, IRF8, SWAP70, ARHGAP24, FCGR2B, SMIM14, PPP1R14A, C16orf74
Negative:  NKG7, PRF1, CST7, GZMB, GZMA, FGFBP2, CTSW, GNLY, B2M, SPON2
           CCL4, GZMH, FCGR3A, CCL5, CD247, XCL2, CLIC3, AKR1C3, SRGN, HOPX
           TTC38, APMAP, CTSC, S100A4, IGFBP7, ANXA1, ID2, IL32, XCL1, RHOC
PC_ 3
Positive:  HLA-DQA1, CD79A, CD79B, HLA-DQB1, HLA-DPB1, HLA-DPA1, CD74, MS4A1, HLA-DRB1, HLA-DRA
           HLA-DRB5, HLA-DQA2, TCL1A, LINC00926, HLA-DMB, HLA-DMA, CD37, HVCN1, FCRLA, IRF8
           PLAC8, BLNK, MALAT1, SMIM14, PLD4, LAT2, IGLL5, P2RX5, SWAP70, FCGR2B
Negative:  PPBP, PF4, SDPR, SPARC, GNG11, NRGN, GP9, RGS18, TUBB1, CLU
           HIST1H2AC, AP001189.4, ITGA2B, CD9, TMEM40, PTCRA, CA2, ACRBP, MMD, TREML1
           NGFRAP1, F13A1, SEPT5, RUFY1, TSC22D1, MPP1, CMTM5, RP11-367G6.3, MYL9, GP1BA
PC_ 4
Positive:  HLA-DQA1, CD79B, CD79A, MS4A1, HLA-DQB1, CD74, HLA-DPB1, HIST1H2AC, PF4, TCL1A
           SDPR, HLA-DPA1, HLA-DRB1, HLA-DQA2, HLA-DRA, PPBP, LINC00926, GNG11, HLA-DRB5, SPARC
           GP9, AP001189.4, CA2, PTCRA, CD9, NRGN, RGS18, GZMB, CLU, TUBB1
Negative:  VIM, IL7R, S100A6, IL32, S100A8, S100A4, GIMAP7, S100A10, S100A9, MAL
           AQP3, CD2, CD14, FYB, LGALS2, GIMAP4, ANXA1, CD27, FCN1, RBP7
           LYZ, S100A11, GIMAP5, MS4A6A, S100A12, FOLR3, TRABD2A, AIF1, IL8, IFI6
PC_ 5
Positive:  GZMB, NKG7, S100A8, FGFBP2, GNLY, CCL4, CST7, PRF1, GZMA, SPON2
           GZMH, S100A9, LGALS2, CCL3, CTSW, XCL2, CD14, CLIC3, S100A12, CCL5
           RBP7, MS4A6A, GSTP1, FOLR3, IGFBP7, TYROBP, TTC38, AKR1C3, XCL1, HOPX
Negative:  LTB, IL7R, CKB, VIM, MS4A7, AQP3, CYTIP, RP11-290F20.3, SIGLEC10, HMOX1
           PTGES3, LILRB2, MAL, CD27, HN1, CD2, GDI2, ANXA5, CORO1B, TUBA1B
           FAM110A, ATP1A1, TRADD, PPA1, CCDC109B, ABRACL, CTD-2006K23.1, WARS, VMO1, FYB

查看对每个主成分影响比较大的基因集

> print(sce [["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(sce , dims = 1:2, reduction = "pca")

VizDimReduction, DimPlot, 和 DimHeatmap可以从基因或细胞角度可视化pca结果

可视化对每个主成分影响比较大的基因集

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

DimPlot(sce , reduction = "pca")

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

DimHeatmap(sce, dims = 1, cells = 500, balanced = TRUE) #1个PC 500个细胞

展示多个主成分

DimHeatmap(sce,dims = 1:15,cells = 500,balanced = TRUE) #15个PC

2、数据维度

主成分分析的原理非常简单，概括来说就是选择包含信息量大的维度（features），去除信息量少的“干扰”维度。所以这里会有个问题——如何知道保留几个维度是最佳的呢？我们希望通过保留尽可能少的维度来留存尽可能多的信息。Seurat有两种方法来确定维度。

JackStrawPlot()函数提供可视化方法，用于比较每一个主成分的p-value的分布，虚线是均匀分布；显著的主成分富集有小p-Value基因，实线位于虚线左上方。下图表明保留10个pca主成分用于后续分析是比较合理的。

JackStraw && Elboew plot

sce <- JackStraw(sce,num.replicate = 100)
sce <- ScoreJackStraw(sce,dims = 1:20)

> JackStrawPlot(sce,dims = 1:15)
Warning message:
Removed 23516 rows containing missing values (`geom_point()`).

主要看在哪个PC之后，显著性大幅下降，也就是前多少个维度包含了大部分的样本信息。可以看出，在10-12个PC之后，显著性大幅下降，也就是前10-12个维度包含了大部分的样本信息。

ElbowPlot(sce)

主要看PC在哪个附近有拐点（“elbow”），拐点表明大部分真实信号是在前多少个pc中捕获的。

可以看出，PC9-10附近有一个拐点（“elbow”），这表明大部分真实信号是在前10个pc中捕获的。

综合以上方法，选择10个主成成分作为参数用于后续分析。

启发式评估方法生成一个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）重复下游分析。其实也将观察到的，结果通常没有显著差异。因此，在选择此参数时，可以尽量选大一点的维度，维度太小的话对结果会产生不好的影响。

七、 细胞聚类

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。

> sce<- FindNeighbors(sce, dims = 1:10) #选取前10个主成分来分类细胞
Computing nearest neighbor graph
Computing SNN


> sce<- FindClusters(sce, resolution = 0.5)
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 2638
Number of edges: 95965

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

#查看细胞属于那个类群
head(Idents(sce), 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

八、非线性降维(UMAP/tSNE)

Seurat 提供了一些非线性降维度分析的方法，如 tSNE 和 UMAP，以可视化和探索这些数据集；这一步使用的数据建议与聚类分析使用的pc维度一致。

# install UMAP： reticulate::py_install(packages ='umap-learn')

sce <- RunUMAP(sce,dims = 1:20)
DimPlot(sce,reduction = 'umap')

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

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

saveRDS(sce, file = "../output/pbmc_tutorial.rds")

九、找差异表达基因 ( 聚类标志cluster biomarkers )

利用 FindMarkers 命令，可以找到各个细胞类型中与其他类别的差异表达基因，作为该细胞类型的生物学标记基因。

• ident.1：设置待分析的细胞类别
• min.pct：设定在两个细胞群中任何一个被检测到的百分比，通过此设定不检测很少表达基因来缩短程序运行时间。默认0.1。
• thresh.test：设定在两个细胞群中基因差异表达量（平均）。可以设置为0 ，程序运行时间会更长。
• max.cells.per.ident：每个类群细胞抽样设置；也可以缩短程序运行时间。通过降低每个类的采样值，提高计算速度

# find all markers of cluster 1
cluster1.markers <- FindMarkers(sce, 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(sce, 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


# 找出每个cluster的标记与所有剩余的细胞相比较，只报告阳性细胞
# find markers for every cluster compared to all remaining cells, report only the positive ones
sce.markers <- FindAllMarkers(sce, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
sce.markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_log2FC)


       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(sce, ident.1 = 0, logfc.threshold = 0.25, test.use = "roc", only.pos = TRUE)
head(cluster1.markers, n = 5)
  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=00s
      myAUC  avg_diff power avg_log2FC pct.1 pct.2
RPS12 0.824 0.5094804 0.648  0.7350248 1.000 0.991
RPS6  0.821 0.4712446 0.642  0.6798622 1.000 0.995
RPS27 0.819 0.4996079 0.638  0.7207819 0.999 0.992
RPL32 0.815 0.4238952 0.630  0.6115515 0.999 0.995
RPS14 0.808 0.4296946 0.616  0.6199183 1.000 0.994

FindAllMarkers（）参数意义：

• only.pos = TRUE：只寻找上调的基因
• min.pct = 0.1：某基因在细胞中表达的细胞数占相应cluster细胞数最低10%
• logfc.threshold = 0.25 ：fold change倍数为0.25
  该函数是决速步，执行比较耗时。

有多种方法可视化标记基因的方法

• VlnPlot: 基于细胞类群的基因表达概率分布
• FeaturePlot：在tSNE 或 PCA图中画出基因表达情况
• RidgePlot，CellScatter，DotPlot

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

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

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

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

#每个聚类前10个差异基因表达热图(如果小于10，则绘制所有标记)
top10 <- sce.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC)
DoHeatmap(sce, features = top10$gene) + NoLegend()

十、鉴定细胞类型

数据集的 markers 可以通过查阅相关文献及网站 (CellMarker) 人工注释，或者利用singleR自动注释。singleR 比较鸡肋，最好的办法目前还是查阅文件和cellmarker(http://biocc.hrbmu.edu.cn/CellMarker/)进行区分，比较考验个人及课题组经验及搜索能力。

比如研究多形性胶质母细胞瘤(GBM)，下图中的文献给出了利用单细胞转录组测序确定胶质母细胞瘤中肿瘤、基质及免疫细胞亚群：通过scRNA测序确定了GBM中CD45-和CD45+细胞(即非免疫和免疫细胞群)。表达EGFR和iCre的簇定为肿瘤细胞。采取Johnson-Verhaak命名法标记EGFR+肿瘤簇。通过marker确定细胞类型，如Oligodendrocytes(Mog，Plp1)，Microglia(P2ry12，Tmem119)，Astrocytes(Aqp4)，Perivascular macrophages(Cd163，Mrc1)，cDC1(Xcr1)等。

那我们就可以根据文中的基因试一下：

FeaturePlot(sce,
            features = c('MOG','PIP1','AQP4','CD163','P2RY12','TMEM119','XCR1'),
            reduction = 'umap')

列表

Cluster ID Markers Cell Type
8 P2ry12 Microglia

十一、细胞注释

Assigning cell type identity to clusters
根据已知细胞类型与基因标记的对应关系，可以为细胞类群找到对应的细胞类型。

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

# 保存分析的结果
saveRDS(sce, file = "../output/sce_final.rds")

官网此部分介绍：
https://satijalab.org/seurat/articles/pbmc3k_tutorial.html

参考：
https://www.jianshu.com/p/728125be7c53
https://www.jianshu.com/p/5b26d7bc37b7
https://zhuanlan.zhihu.com/p/359020041
https://www.jianshu.com/p/728125be7c53
https://www.jianshu.com/nb/52917361
https://www.jianshu.com/p/fef17a1babc2
Single cell