单细胞分析 | Seurat基础流程 | 保姆级教程

来自牧星人小伙伴（小杨想去看极光）的投稿，投稿目的是督促自己学习，突然觉得自己又多了一丢丢价值哈哈哈哈哈哈哈！开心！我觉得他写得超棒的，大家一起加油！！！冲！

单细胞的内容我鸽的有点久呜呜呜呜呜呜呜呜，给大家道个歉！抱歉最近实在是太太太太太忙了！！！咱们今天就从这里开启吧嘿嘿嘿！臭不要脸地借花献佛咯！

---

• 如果小伙伴们有需求的话，可以加入我们的交流群：一定要知道 | 永久免费的环境友好型生信学习交流群又双叒叕来啦！| 伴随不定期群友好物分享！在这里，你可以稍有克制地畅所欲言！

• 超级建议大家在入群前或入群后可以看一下这个：干货满满 | 给生信小白的入门小建议 | 掏心掏肺版！绝对干货满满！让你不虚此看！

• 如果有需要个性化定制分析服务的小伙伴，可以看看这里：你要的个性化生信分析服务今天正式开启啦！定制你的专属解决方案！全程1v1答疑！！绝对包你满意！

---

Seurat的分析是从单细胞计数矩阵开始的。单细胞计数矩阵是由cellranger流程对下机数据进行处理而产生的，里面存储的是UMI计数矩阵，矩阵中的值表示的是在每个细胞（列）中检测到的每个特征（即基因，行）的分子数。

00 加载分析所需要的包

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

01 创建Seurat对象

Seurat的对象将作为一个容器，其中包含单细胞数据集的数据和分析的结果。

加载数据集

首先，需要加载我们的数据集：

# 设置目录
setwd("D:\\Document\\Programe\\RProjrct\\R-4.2.3")
# 加载测试数据集
pbmc.data <- Read10X(data.dir = "./filtered_gene_bc_matrices/hg19/")

数据集文件解析

测试数据集中包含barcodes.tsv、genes.tsv、matrix.mtx三个文件。如下图所示：

在解释三个文件之前，我们先统计一下这三个文件的行数，一会儿有用哟！

$ wc -l *
    2700 barcodes.tsv
   32738 genes.tsv
 2286887 matrix.mtx
 2322325 total

接下来，我们开始解析每个文件：

1.barcodes.tsv文件存储的是用于标记不同细胞的条形码序列，即barcode，只有一列信息，由于其标记的是细胞，所以该文件的行数对应的就是细胞数量，由行数统计结果可知细胞数量为2700个；

$ head barcodes.tsv
AAACATACAACCAC-1
AAACATTGAGCTAC-1
AAACATTGATCAGC-1
AAACCGTGCTTCCG-1
AAACCGTGTATGCG-1

2.genes.tsv文件存储的是基因的信息，分别为gene id（第一列）和gene symbol（第二列），列数不固定，主要作为矩阵的行名，可进行匹配识别基因即可。与barcodes.tsv文件同理，该文件的行数对应的是基因数量，基因数量为32738个；

$ head genes.tsv
ENSG00000243485  MIR1302-10
ENSG00000237613  FAM138A
ENSG00000186092  OR4F5
ENSG00000238009  RP11-34P13.7
ENSG00000239945  RP11-34P13.8

3.matrix.mtx文件存储的主要是有表达量的基因的counts数量（第三列信息，即表达量信息），前三行为注释信息，其中第三行跟我们上面统计的每个文件的行数信息基本对应，实际上这三个数字对应的便分别是基因数量、细胞数量、有表达量的基因的数量（任一基因在任一细胞中出现都算一次计数），去掉三行注释行，刚好是2286884行。

$ head 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

我看别人说，这种方式可以高效地存储数据。我们随便找一行举例，看看这个文件怎么解析啊。比如"32706 1 10"这一行，表示的是barcodes.tsv文件中的第1个细胞的在genes.tsv文件中的第32706个基因的counts数为10。

创建Seurat对象

接下来，我们就来创建一个Seurat对象：

# 创建Seurat对象
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
#查看对象信息

CreateSeuratObject的几个主要参数：

• counts：输入对应的矩阵信息
• project：Seurat对象的项目名
• min.cells：规定基因的表达范围，即单个基因至少在多少个细胞中表达方可保留
• min.features：规定细胞表达的基因数范围，即单个细胞中至少表达多少个基因方可保留

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

我们可以看到数据集变成了一个2700个细胞，13714个基因的矩阵，说明上面过滤掉了一部分基因的信息，细胞的表达的基因数量都高于过滤阈值，所以细胞数量并没有减少。

我们也可以简单看一下不加这些参数会是什么样！

# 创建Seurat对象，不加过滤条件
pbmc <- CreateSeuratObject(counts = pbmc.data)
pbmc
# An object of class Seurat 
# 32738 features across 2700 samples within 1 assay 
# Active assay: RNA (32738 features, 0 variable features)

确实是所有基因和细胞的信息都被保留了下来。

补充：那么矩阵加载完成后，矩阵中的数据到底是什么样的呢？

# 我们参照官网示例，看看部分基因在前30个细胞中的表达矩阵信息。
# 三个基因：CD3D、TCL1A、MS4A1在前30个细胞中的表达矩阵信息
pbmc.data[c("CD3D", "TCL1A", "MS4A1"), 1:30]
# 3 x 30 sparse Matrix of class "dgCMatrix"
#   [[ suppressing 30 column names ‘AAACATACAACCAC-1’, ‘AAACA# TTGAGCTAC-1’, ‘AAACATTGATCAGC-1’ ... ]]
#                                                             #        
# 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，表示未在该细胞中检测到该基因。

为什么要用.来代表0呢？

这是由于scRNA-seq矩阵中的大多数值都是0，0属于无效信息，存储他们的值和位置信息需要占用大量内存，而Seurat尽可能使用稀疏矩阵表示，只需要存储非0数据的值和位置信息，因此，这将大大节省计算数据的内存和速度。

这样讲太笼统了，我们直接上数据：

# 密集矩阵
dense.size <- object.size(as.matrix(pbmc.data))
dense.size
# 709591472 bytes

# 稀疏矩阵
sparse.size <- object.size(pbmc.data)
sparse.size
# 29905192 bytes

# 两者相差倍数
dense.size/sparse.size
# 23.7 bytes

可以看到使用稀疏矩阵存储数据可以极大的减少数据的内存占用（23.7倍），提升运行的速度。

02 标准预处理工作流程

该步骤主要是进行质控以及选择细胞进行进一步分析。

Seurat中一些常用的质控（QC）指标：

nFeature_RNA：每个细胞中检测到的基因的数量。

• 低质量的细胞或空液滴通常含有很少的基因；
• 细胞双胞体或多胞体可能表现出异常高的基因计数。

nCount_RNA：细胞内检测到的分子总数。即检测到的基因的总数量（一个基因可能会因为有多个转录本而被检测到多次），因此分子总数与基因数量密切相关。

percent.mt：对应线粒体基因组的序列的百分比。低质量/濒死的细胞通常表现出广泛的线粒体污染。我们通常使用PercentageFeatureSet()函数来计算线粒体基因相关的reads的百分比。代码如下：

pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

"[[]]"运算符可以向对象pbmc的元数据meta.data中添加列。meta.data是一个存储QC统计数据的好地方。关于seurat V4的数据结构以及和V5版本数据结构的区别后续会做讲解。

我们先来看一下计算这一步之前的meta.data长什么样子：

# 计算前的meta.data
# head(pbmc@meta.data)
#                  orig.ident nCount_RNA nFeature_RNA
# AAACATACAACCAC-1     pbmc3k       2419          779
# AAACATTGAGCTAC-1     pbmc3k       4903         1352
# AAACATTGATCAGC-1     pbmc3k       3147         1129
# AAACCGTGCTTCCG-1     pbmc3k       2639          960
# AAACCGTGTATGCG-1     pbmc3k        980          521
# AAACGCACTGGTAC-1     pbmc3k       2163          781

这里可能有人会问，这里面的信息除了第一列orig.ident好像是在CreateSeuratObject()的时候定义的project，后面两个参数nCount_RNA和nFeature_RNA是哪来的呢？这是由于在CreateSeuratObject()过程中会自动计算基因的数量nFeature_RNA和分子总数nCount_RNA，他们会被存储在对象的meta.data中。

我们再来看一下计算这一步之后的meta.data：

# 计算后的meta.data
> pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
> head(pbmc@meta.data)
#                  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
# AAACGCACTGGTAC-1     pbmc3k       2163          781  1.6643551

多了这样一列 percent.mt的数据，存储的便是每个细胞中线粒体基因的占比情况。

注：上面是作为示例数据，他的基因格式比较标准。

本教程中的线粒体基因如下：

rownames(pbmc)[grep("^MT-", rownames(pbmc))]
# [1] "MT-ND1"  "MT-ND2"  "MT-CO1"  "MT-CO2"  "MT-ATP8" "MT-ATP6" "MT-CO3"  "MT-ND3"  "MT-ND4L" "MT-ND4" "MT-ND5"  "MT-ND6"  "MT-CYB"

那么当我们使用常规的数据，只知道线粒体基因的id，他们没有共同特征时该怎么办呢？方法就是把pattern展开来写，只有13个基因，稍微费点事就写下来了：

pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "ATP6|ATP8|COX1|COX2|COX3|CYTB|ND1|ND2|ND3|ND4|ND4L|ND5|ND6")

那还有人会问了，如果不知道线粒体基因的信息怎么办呢？这个线粒体基因信息在gtf或者gff文件中都会有的，你去找一下基本都会找到的，找到后再像上面一样匹配就好啦！

上面我们介绍完质控的这些指标，那么就进入正题，开始质控吧！

设置阈值

首先呢，我们要了解如何设置质控指标的阈值。

质控的目的是挑选质量更好的细胞进行进一步的分析，那么我们如何有针对性地进行挑选呢？当然是需要在质控前先对质控的指标进行一波统计，作为设置阈值的参考。统计的方法很简单，小提琴图VlnPlot和散点图FeatureScatter （VlnPlot这个函数有很多的功能，这里就先略过，下次专门整理！），我们就按照官网上的来：

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

第一个图nFeature_RNA是对所有细胞检测到的基因数量的小提琴图，横坐标为样品名，纵坐标为每个细胞中包含的基因数量，图中的每个点代表一个细胞；

第二个图nCount_RNA为样品所有细胞内检测到的分子总数的小提琴图，横坐标为样品名，纵坐标为每个细胞内检测到的分子总数，图中的每个点代表一个细胞。

第三个图percent.mt为样品所有细胞的线粒体基因比例的小提琴图，横坐标为样品名，纵坐标为每个细胞的线粒体基因比例，图中的每个点代表一个细胞。

小提琴图展示了任意位置数值的密度，可以知道哪些位置的数值的密度较高。

# FeatureScatter是一种典型的可视化特征与特征之间关系的方法
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2

随着测序深度的增加，单细胞所检测到的基因数量也在增加，两者是有一定关联性的。

对含有线粒体基因的细胞进行筛选是因为线粒体基因比例过高的细胞通常为低质量的细胞，这些数据杂糅进来会干扰细胞分群。

第一个图nCount_RNA vs percent.mt，是线粒体基因比例与细胞中检测到的分子总数的关系的散点图，通常两者没有什么相关关系，我们通常对这两个指标的标准是：线粒体基因比例应相对较低以排除低质量/濒死的细胞，细胞中检测到的分子总数也不宜过高以排除双胞和多胞细胞的影响。

第二个图nCount_RNA vs nFeature_RNA，是基因数量与细胞中检测到的分子总数的关系的散点图，即测序深度与基因数量的关系。高质量的测序数据中两者基本处于正相关的关系，但要排除由于双胞和多胞造成的分子数量过大的部分数据，即右上方离群点。

开始过滤

接着呢，我们就开始过滤！

根据小提琴图和散点图的情况，个人认为可以选择如下的质控过滤指标，即：

• nFeature_RNA > 200（最开始的min.features = 200）& nFeature_RNA < 2000
• nCount_RNA < 7000
• percent.mt < 7
  但是呢，为了和官方教程相统一，我们还是采用官方的过滤阈值，实际上相差并不多!

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

过滤完成

最后，过滤完成，让我们检查一些指标看看吧！

# 细胞数量
pbmc
# An object of class Seurat 
# 13714 features across 2638 samples within 1 assay 
# Active assay: RNA (13714 features, 0 variable features)

可以看到这一通简单的过滤之后，细胞数量由2700减少到了2638，基因数量不变。

# 小提琴图可视化QC指标
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
# FeatureScatter可视化QC指标之间的关系
plot1 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2

03 数据标准化

通过质控过滤完不需要的细胞之后呢，就要对数据进行标准化处理了！

LogNormalize：全局缩放的标准化方法（其实还有另一种标准化方法SCTransform，随后咱单独整理一篇！）

1. 为什么要做normalization呢？

   因为细胞间的异质性除了可能来自于细胞自身基因的差异表达，还有可能是由测序技术误差造成的。数据标准化的目的便是排除技术误差，让测序深度和文库大小不同的细胞的基因表达量具有可比性，进而更好地呈现基因水平的差异。

2. 计算方法

   we employ a global-scaling normalization method “LogNormalize” that normalizes the feature expression measurements for each cell by the total expression, multiplies this by a scale factor (10,000 by default), and log-transforms the result.

   我们直接上公式，这个简单易懂，其中的log1p指的是log(x + 1)，只看到log就表示以e为底的log值。

   Value：该基因在该细胞中的表达量

   

3. 执行代码

   # scale.factor这里我们选择10000
   pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
   # 结果储存在pbmc[["RNA"]]@data中，可以通过pbmc[["RNA"]]@data查看
   head(pbmc[["RNA"]]@data[1:6, 1:25])
   # 6 x 25 sparse Matrix of class "dgCMatrix"
   #   [[ suppressing 25 column names ‘AAACATACAACCAC-1’, ‘AAACATTGAGCTAC-1’, ‘AAACATTGATCAGC-1’ ... ]]
                                                                            
   # AL627309.1    . . . . . . . . . . . .        . . . . . . . . .        . . . .
   # AP006222.2    . . . . . . . . . . . .        . . . . . . . . .        . . . .
   # RP11-206L10.2 . . . . . . . . . . . .        . . . . . . . . .        . . . .
   # RP11-206L10.9 . . . . . . . . . . . .        . . . . . . . . .        . . . .
   # LINC00115     . . . . . . . . . . . .        . . . . . . . . .        . . . .
   # NOC2L         . . . . . . . . . . . 1.646272 . . . . . . . . 1.398186 . . . .
   
   # 我们对比一下标准化之前的结果
   head(pbmc@assays$RNA@counts[1:6, 1:25])
   # 6 x 25 sparse Matrix of class "dgCMatrix"
   #   [[ suppressing 25 column names ‘AAACATACAACCAC-1’, ‘AAACATTGAGCTAC-1’, ‘AAACATTGATCAGC-1’ ... ]]
                                                              
   # AL627309.1    . . . . . . . . . . . . . . . . . . . . . . . . .
   # AP006222.2    . . . . . . . . . . . . . . . . . . . . . . . . .
   # RP11-206L10.2 . . . . . . . . . . . . . . . . . . . . . . . . .
   # RP11-206L10.9 . . . . . . . . . . . . . . . . . . . . . . . . .
   # LINC00115     . . . . . . . . . . . . . . . . . . . . . . . . .
   # NOC2L         . . . . . . . . . . . 1 . . . . . . . . 1 . . . .
   

那么，问题来了，为什么原始数据都是1，在标准化之后数值就不一样了呢？原因很简单，因为不同细胞的总counts（即all_gene_counts）发生了改变，导致最终的值发生了改变。

04 鉴定高可变基因

我们接下来计算在数据集中表现出细胞间高差异性的基因子集（即它们在一些细胞中高度表达，在另一些细胞中低度表达）。研究发现，在下游分析中这些基因有助于突出单细胞数据集中的生物信号。

FindVariableFeatures：鉴定高可变基因的函数，默认返回每个数据集的前2000个基因用于像PCA这样的下游分析。这一部分没啥好说的，后面分析就主要用这2000个基因了，这样可以减小内存和计算资源。

# FindVariableFeatures鉴定高可变基因
pbmc <- FindVariableFeatures(pbmc, 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(pbmc), 10)

# 绘制高可变基因的火山图
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1 + plot2

05 数据缩放（数据中心化）

在计算完高可变基因后，需要对数据集进行一个线性转换（缩放）。

这个步骤是降维（如PCA）数据处理之前的一个标准的预处理过程，由ScaleData()这个函数执行。该函数的作用解析如下：

• 调整每个基因的表达量，使得细胞间的平均表达量为0
• 缩放每个基因的表达量，使得细胞间的方差为1
• 目的：这一步骤赋予每个基因在下游分析中同样的权重，不至于使高表达基因占据主导地位
• 结果存储在pbmc[[“RNA”]]$scale.data中

运行代码如下：

# 将Seurat对象的行名（所有基因）抓取出来
all.genes <- rownames(pbmc)
# 以基因为单位对矩阵数据进行缩放
pbmc <- ScaleData(pbmc, features = all.genes)
# Centering and scaling data matrix
#   |=======================================================================================| 100%

总结：概括起来讲，scale的作用是使数据具有一个统一的尺度，这并不会影响每个点的相对位置，只是使他们的表达量尺度统一起来。是PCA分析的必要步骤。

官网关于这一步讨论了两个问题：

1. 第一个问题是这一步的运行时间太长，是否可以加速运行？
   答案是肯定的。首先我们要明确的一点是虽然缩放是Seurat分析过程中不可缺少的一步，但是也只是需要将PCA分析需要的基因作为输入即可，因此，在这一点上ScaleData()默认使用的是先前鉴定的2000个高变基因。所以将features参数去掉，将all.genes变成高变基因即可减少运行时间。

   pbmc <- ScaleData(pbmc)
   # 也可以这样运行
   scale.genes <-  VariableFeatures(pbmc)
   pbmc <- ScaleData(pbmc, features = scale.genes)
   

2. 第二个问题是如何删除不需要的变异来源？

   pbmc <- ScaleData(pbmc, vars.to.regress = "percent.mt")
   # 如果希望使用此功能，建议使用新的标准化方法SCTransform()
   # 关于这个问题，我的了解较少，在此处不做过多展开。
   

06 PCA线性降维

PCA分析

接下来，对缩放的数据进行PCA分析。

通常情况下，使用上述计算的高变基因进行PCA分析，基因集的选择可以通过features参数进行设置。

pbmc <- RunPCA(pbmc, features = VariableFeatures(pbmc))
# 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(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

接下来进行PCA结果的可视化（其实不太重要）。

Seurat提供了几种方法来对PCA结果的细胞和基因进行可视化，包括VizDimReduce()、 DimPlot()和DimHeatmap()。

第一种：VizDimLoadings()：VizDimLoadings是一个R包中的函数，用于绘制维度加载图。维度加载图是一种用于可视化主成分分析（PCA）或其他降维方法中的变量与主成分之间的关系的图形。它显示了每个变量在每个主成分上的贡献程度，可以帮助我们理解变量在不同主成分上的权重和相关性。这种图形通常以散点图或箭头图的形式呈现。

# 绘制前两个主成分，降维方式采用pca方法。
VizDimLoadings(pbmc, dims = 1:2, reduction = "pca") # 没啥用其实

第二种：DimPlot()：DimPlot函数可以将单细胞数据降维后的结果可视化，以便于我们观察不同细胞之间的相似性和差异性。在散点图中，每个点代表一个单细胞，不同颜色或形状的点表示不同的细胞类型或状态。通过观察散点图，我们可以发现细胞之间的聚类模式、相似性和差异性，进而对细胞类型或状态进行分类和注释。此外，DimPlot函数还可以根据其他降维方法（如t-SNE和UMAP）进行可视化，以便于我们比较不同降维方法的效果。

DimPlot(pbmc, reduction = "pca") # 也没啥用，但是分析基本都会展示一下

第三种：DimHeatmap()：DimHeatmap函数可以将单细胞数据在降维空间中的表达模式可视化，以便于我们观察不同基因或基因集之间的表达模式和差异性。在热图中，每一列代表一个细胞，每一行代表一个基因或基因集，颜色的深浅表示该基因或基因集在该细胞中的表达水平。通过观察热图，我们可以发现不同基因或基因集之间的相关性、表达模式和差异性，进而对细胞类型或状态进行分类和注释。在DimHeatmap函数中，我们可以指定降维的主成分数量和用于绘制热图的细胞数量，以便于我们控制可视化结果的精度和速度。此外，如果存在批次效应或其他技术因素导致的表达偏差，可以使用平衡化方法来消除这些效应，从而提高可视化结果的准确性和可靠性。

DimHeatmap(pbmc, dims = 1:15, cells = 500, balanced = TRUE) # 也没啥用

07 确定后续分析需要的“维度”

为了克服scRNA-seq数据的任何单一基因的广泛技术噪音，Seurat基于PCA分数对细胞进行聚类，每个PC基本上代表一个“metafeature”（metafeature可以理解为一组基因，这群基因具有相似的表达pattern），因此，数据集的前几个主成分被稳定（robust）地压缩（上面这一句是官网翻译过来的，Robust：稳健的，音译为鲁棒性，鲁棒性是指在面对异常情况或噪声时，系统或方法仍能保持稳定和可靠的能力。在数据分析中，鲁棒性通常指的是某个统计模型或算法对于数据中的离群点、缺失值、错误值等异常情况的容忍度。具有较高鲁棒性的模型或算法能够有效地识别和处理异常值，从而提高分析结果的准确性和可靠性）。因此，如何准确地确定降维的“维度”对后续的分析非常重要！

通常情况下，我们使用最多的就是肘形图（ElbowPlot），原因很简单，就是结果展现形式很直观（在V4官网上还有一种JackStraw图，肘形图实际上也来源于这种图，只不过JackStraw图并不直观，个人感觉在实际操作中也很少人用，这里就不做阐述啦）。

ElbowPlot(pbmc)

关于这个图的解析，官网给的解析是这样的：基于每个元素(PC)解释的方差百分比对主要元素进行排序。在这个例子中，我们可以观察到在 PC9-10周围有一个“肘”，这表明大部分真实信号是在前10个 PC 中捕获的。

我认为可以这样解释，横坐标的PC是你接下来要选择分析的维度数，即多少组基因表达pattern。当选择一组的时候，即所有基因都在一个表达pattern里面（我们从生物学角度看也肯定是不可能的），我们很容易想到，所有细胞都在一个组内，肯定是存在很大差异的，所以方差值（方差=标准差（Standard Deviation）的平方）就会很大。当PC增多，即分组增多后，每个组内的差异也就慢慢变小了，方差值自然也就变小了。再往后，随着选的PC的增多，基本上区分不出来组内的差异了，这个时候方差基本就稳定了下来，此时在选取更大的PC值也就没有更大的意义了，并且还会增加计算压力。因此，个人认为关于PC选取的最合适原则是在拐点右侧3-5个维度即可（本数据拐点在7-8左右，我们选10刚刚好），宽松一点的话，可以设置在拐点右侧10以内。不过最重要的一点还是尽可能地多尝试，选取最适合自己的PC值，这一点也是官网教程所建议的。

08 细胞聚类

Seurat 应用了一种基于图表的聚类方法。重要的是，驱动聚类分析的距离度量(基于先前识别的 PC)保持不变。将细胞嵌入到图形结构中-例如 K 最近邻(KNN)图，在具有相似特征表达模式的细胞之间绘制边界，然后将该图划分为高度相互关联的不同簇。

首先基于PCA空间中的欧几里得度量构建一个KNN图，然后基于局部邻域中的共享重叠(Jaccard 相似性)来精确任意两个细胞之间的边缘权重。这个步骤使用FindNeighbors()函数执行，并将先前定义的前10个PC作为输入。

pbmc <- FindNeighbors(pbmc, dims = 1:10)
# Computing nearest neighbor graph
# Computing SNN

为了将细胞聚类，我们接下来应用模块化优化技术，如Louvain算法(默认值)或SLM，迭代地将细胞聚类在一起。该步骤通过FindClusters()函数实现，其包含一个解析参数resolution，可以理解为就是分辨率，分辨率数值越高，对细胞分群就越细，分到的群就可能会越多；反之则越少。官网这个3K的细胞集，基本上分辨率会选0.4-1.2.然后还跟我们说，对于较大的数据集，最佳分辨率往往会提高。

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

# 分组情况可以使用Idents()函数进行查看
head(Idents(pbmc))
# AAACATACAACCAC-1 AAACATTGAGCTAC-1 AAACATTGATCAGC-1 AAACCGTGCTTCCG-1 AAACCGTGTATGCG-1 
#                2                3                2                1                6 
# Levels: 0 1 2 3 4 5 6 7 8

09 执行非线性降维（UMAP/tSNE）

Seurat 提供了一些非线性降维技术，如 tSNE 和 UMAP，来可视化和探索这些数据集。这些算法的目标是学习数据集中的底层结构，以便在低维空间中将相似的细胞放在一起。虽然像 tSNE 和 UMAP 这样的二维可视化技术是探索数据集的有价值的工具，但是所有的可视化技术都有局限性，不能完全代表底层数据的复杂性。特别是，这些方法旨在保持数据集中的局部距离(即确保具有非常相似的基因表达谱的细胞共定位) ，但通常不能保持更多的全局关系。（这些是官网翻译过来的原话，关于这一部分没啥好说的，主要作用就是可视化）。

# 官方推荐用UMAP，咱先来UMAP，先降维后可视化
pbmc <- RunUMAP(pbmc, dims = 1:10)
# 11:00:02 UMAP embedding parameters a = 0.9922 b = 1.112
# 11:00:02 Read 2638 rows and found 10 numeric columns
# 11:00:02 Using Annoy for neighbor search, n_neighbors = 30
# 11:00:02 Building Annoy index with metric = cosine, n_trees = 50
# 0%   10   20   30   40   50   60   70   80   90   100%
# [----|----|----|----|----|----|----|----|----|----|
# **************************************************|
# 11:00:02 Writing NN index file to temp file C:\Users\YangQi\AppData\Local\Temp\RtmpEzVJIW\file3af0f97a09
# 11:00:02 Searching Annoy index using 1 thread, search_k = 3000
# 11:00:03 Annoy recall = 100%
# 11:00:05 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
# 11:00:06 Initializing from normalized Laplacian + noise (using irlba)
# 11:00:06 Commencing optimization for 500 epochs, with 105124 positive edges
# Using method 'umap'
# 0%   10   20   30   40   50   60   70   80   90   100%
# [----|----|----|----|----|----|----|----|----|----|
# **************************************************|
# 11:00:14 Optimization finished

DimPlot(pbmc, reduction = "umap")

# 咱再试试tSNE
pbmc <- RunTSNE(pbmc, dims = 1:10)
DimPlot(pbmc, reduction = "umap")

关于UMAP和TSNE图的展示简单做一个描述：可以说UMAP展现的距离更真实，TSNE图展示的结果似乎更美观（TSNE图呈现出来的结果中，各个细胞团都很紧促，最终呈现出来细胞都抱团分布，出来的图个人感觉更美观一丢丢）。

# 保存结果为rds文件（这里虽然可以保存，但是并非最好的保存位置）
saveRDS(pbmc, file = "../output/pbmc_tutorial.rds")

10 寻找差异表达基因

Seurat可以通过差异表达(DE)来获取用于定义簇的标记。默认情况下，与所有其他细胞簇相比，它标识单个簇(在 id.1中指定)的正标记和负标记。FindAllMarkers()为所有簇自动完成这个过程，但是您也可以对任意簇之间进行比较，或者是某一簇对其他所有簇进行比较。（官网的话术比较专业，这一步没什么深奥的原理，我们自己写几个就知道了）。

计算所有簇的markers

首先第一个，我们就算所有簇的marker基因吧，计算所有簇需要出动FindAllMarkers()，这里官网教程仅设置了一个参数，就是only.pos，这个参数的意思是只保留正标记（positive markers），换句话说也就是只保留相对于其他簇而言，目标簇显著高表达的基因标记。

pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE)
# Calculating cluster 0
# |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=03s  
# Calculating cluster 1
# |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=07s  
# Calculating cluster 2
# |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=03s  
# Calculating cluster 3
# |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=03s  
# Calculating cluster 4
# |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=04s  
# Calculating cluster 5
# |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=07s  
# Calculating cluster 6
# |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=09s  
# Calculating cluster 7
# |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=07s  
# Calculating cluster 8
# |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=06s 

看这个计算过程就知道，每个簇都算一遍，所以簇数多的话，这一步会比较费时间。接下来，就把他们按簇展示一下！

pbmc.markers %>% group_by(cluster)
# A tibble: 4,667 × 7
# Groups:   cluster [9]
#         p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene 
#         <dbl>      <dbl> <dbl> <dbl>     <dbl> <fct>   <chr>
#   1 1.81e-144      0.735 1     0.991 2.48e-140 0       RPS12
#   2 7.14e-142      0.680 1     0.995 9.79e-138 0       RPS6 
#   3 5.26e-140      0.721 0.999 0.992 7.21e-136 0       RPS27
#   4 4.23e-136      0.612 0.999 0.995 5.80e-132 0       RPL32
#   5 1.80e-130      0.620 1     0.994 2.47e-126 0       RPS14
#   6 5.51e-123      0.744 0.997 0.975 7.55e-119 0       RPS25
#   7 1.41e-117      0.609 1     0.994 1.94e-113 0       RPS3 
#   8 3.66e-117      0.744 0.994 0.971 5.02e-113 0       RPL9 
#   9 7.37e-117      0.565 1     0.996 1.01e-112 0       RPL13
#  10 2.75e-116      0.740 0.996 0.964 3.77e-112 0       RPL31
# ℹ 4,657 more rows

# 看起来还不少，4657个基因呢

head(pbmc.markers %>% group_by(cluster))
# A tibble: 6 × 7
# Groups:   cluster [1]
#       p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene 
#       <dbl>      <dbl> <dbl> <dbl>     <dbl> <fct>   <chr>
# 1 1.81e-144      0.735 1     0.991 2.48e-140 0       RPS12
# 2 7.14e-142      0.680 1     0.995 9.79e-138 0       RPS6 
# 3 5.26e-140      0.721 0.999 0.992 7.21e-136 0       RPS27
# 4 4.23e-136      0.612 0.999 0.995 5.80e-132 0       RPL32
# 5 1.80e-130      0.620 1     0.994 2.47e-126 0       RPS14
# 6 5.51e-123      0.744 0.997 0.975 7.55e-119 0       RPS25

tail(pbmc.markers %>% group_by(cluster))
# A tibble: 6 × 7
# Groups:   cluster [1]
#     p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene  
#     <dbl>      <dbl> <dbl> <dbl>     <dbl> <fct>   <chr> 
# 1 0.00791      0.996 0.214 0.054         1 8       PRKCD 
# 2 0.00841      1.10  0.286 0.09          1 8       ARF3  
# 3 0.00872      1.36  0.357 0.134         1 8       TUBB4B
# 4 0.00955      1.01  0.286 0.09          1 8       CERS2 
# 5 0.00985      0.285 0.143 0.027         1 8       CTDP1 
# 6 0.00991      1.12  0.143 0.029         1 8       DGCR2 

计算特定簇间的markers

# 只计算一下第2簇的markers
cluster2.markers <- FindMarkers(pbmc, ident.1 = 2)
# |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=06s 

head(cluster2.markers, n = 5)
#             p_val avg_log2FC pct.1 pct.2    p_val_adj
# IL32 2.593535e-91  1.2154360 0.949 0.466 3.556774e-87
# LTB  7.994465e-87  1.2828597 0.981 0.644 1.096361e-82
# CD3D 3.922451e-70  0.9359210 0.922 0.433 5.379250e-66
# IL7R 1.130870e-66  1.1776027 0.748 0.327 1.550876e-62
# LDHB 4.082189e-65  0.8837324 0.953 0.614 5.598314e-61

# 计算一下第5簇和0，3两簇的差异基因
cluster5.markers <- FindMarkers(pbmc, ident.1 = c(5,6), ident.2 = c(0, 3))
# |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=09s

# 也可以算某几个簇与另外几个簇的差异基因
cluster5.markers <- FindMarkers(pbmc, ident.1 = c(5,6), ident.2 = c(0, 3))
# |++++++++++++++++++++++++++++++++++++++++++++++++++| 100% elapsed=09s

FindMarkers中的一些参数

以下面这个代码为例，计算第0簇的markers。
logfc.threshold：表示基因在两组细胞间的平均表达量的差异倍数，该参数的默认值设置在0.1，增加该值可以加速运算速度，但可能会错过较弱的信号。通常应用中我使用的是0.25。
test.use：对数据进行差异表达测试，可以理解为参数检验，一般笔者采用的方法是wilcox，筛选阈值就是大家常用的p值，也就是0.05。

cluster0.markers <- FindMarkers(pbmc, ident.1 = 0, logfc.threshold = 0.25, test.use = "roc", only.pos = TRUE)

保存markers

# 先把这个markers存入变量，写入变量为了后面好用
All_markers <- pbmc.markers %>% group_by(cluster)

# 保存
write.table(All_markers[All_markers$p_val_adj <= 0.05,], file = "./pbmc_all_markers.xls", sep = "\t", col.names = TRUE, row.names = F, quote = F)

实际上，到了这里单细胞Seurat流程的分析基本上走完了，在这里保存RDS是比较合适的，因为里面也存好了计算的差异基因的情况，到时候重新加载RDS文件时，不需要运算即可获取。同时与注释获得完整的细胞类型还有一定的时间要求。所以，我们选择在这里保存第一版分析数据的RDS。

# 保存结果为rds文件
saveRDS(pbmc, file = "../output/pbmc_tutorial.rds")

11 基因表达可视化

这一部分是专门添加进来的内容，特意和第10步的内容相区别。关于基因表达可视化，俺有五种方法，任君挑选。至于图的含义，我就不做阐述了，已经挺直观啦！And，后续我们也会出基因表达可视化绘图的内容，可以期待一下！

山脊图

features <- c("LYZ", "CCL5", "IL32", "PTPRCAP")
RidgePlot(pbmc, features = features, ncol = 2)
# Picking joint bandwidth of 0.319
# Picking joint bandwidth of 0.178
# Picking joint bandwidth of 0.161
# Picking joint bandwidth of 0.151

小提琴图

VlnPlot(pbmc, features = features, ncol = 2) 

大家有没有发现山脊图和小提琴图有什么关系嘞？把其中一个图旋转90度试试看哈哈哈哈哈哈哈哈哈哈哈！

FeaturePlot

FeaturePlot(pbmc, features = features) 
# 这个图很惭愧，我画的最多，却不知道他中文名应该叫什么

点阵图

DotPlot(pbmc, features = features) + RotatedAxis() 
# 在文章中标注markers的时候很好用，一个图能放下很多基因。

热图

DoHeatmap(subset(pbmc, downsample = 100), features = features, size = 3) 
# 也是一个展示所有注释结果的markers的好形式。

这里记住一点，基因表达量的图是用来显示基因在各簇中的表达情况的，这个为该基因是否可以作为本簇的marker提供了直接的证据来源，如和已有报道相互印证（如XXX细胞特异性表达该基因），则可通过该基因对某一簇细胞进行细胞类型注释。细胞类型注释是一个人工主观的过程（自动注释的应用范围很少，只聚焦在极少数物种中，大多数的细胞类型注释还是依赖人工进行的），这一过程已游离于上述前10步分析之外，所以这就是为什么我在第10步结束的时候保存.rds的原因，不过这里也确实没有细讲如何进行细胞类型注释，咱们后续专门出一期和大家唠唠细胞注释的内容！

12 细胞簇命名

假定我们走完上述11步，已经完成了细胞类型注释工作，即知道了每个簇对应的细胞类型。那么现在我们要如何将这些簇更改成对应的细胞名称呢？看我手法：

# 把每个簇对应的新名字写成一个集合
# 这里注意，如果有两个簇注释到的名字一样，
# 要在对应数字位置都写一遍，
# 即一定保证这个集合的元素数和簇的个数是一致的
new.cluster.ids <- c("Naive CD4 T", "CD14+ Mono", "Memory CD4 T", "B", "CD8 T", "FCGR3A+ Mono", "NK", "DC", "Platelet")

# 我们看一下现在的new.cluster.ids里面是什么
new.cluster.ids
# [1] "Naive CD4 T"  "CD14+ Mono"   "Memory CD4 T" "B"            "CD8 T"        "FCGR3A+ Mono"
# [7] "NK"           "DC"           "Platelet"    

# 输出的结果显示的是9个单独的元素

# 这一步是把pbmc的levels传输给new.cluster.ids作为他的name属性
names(new.cluster.ids) <- levels(pbmc)

# 我们看一下现在的new.cluster.ids里面是什么
new.cluster.ids
#              0              1              2              3              4 
#  "Naive CD4 T"   "CD14+ Mono" "Memory CD4 T"            "B"        "CD8 T" 
#              5              6              7              8 
# "FCGR3A+ Mono"           "NK"           "DC"     "Platelet" 

attributes(new.cluster.ids)
# $names
# [1] "0" "1" "2" "3" "4" "5" "6" "7" "8"

# 明显new.cluster.ids中添加了names属性，每一个names对应了一种细胞类型信息

# 对pbmc的簇进行重命名
pbmc <- RenameIdents(pbmc, new.cluster.ids)

# 可以看到levels中的信息更新成了细胞名称

levels(pbmc)
# [1] "Naive CD4 T"  "CD14+ Mono"   "Memory CD4 T" "B"            "CD8 T"       
# [6] "FCGR3A+ Mono" "NK"           "DC"           "Platelet"    
# 最后，再展示一下注释后的细胞图谱
DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()

到这里，Seurat的基础分析算全部完成了，这时候我们可以保存一下最完整的信息。

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

本教程的内容是按照官网上的V4教程编写的，也是因为V4版本使用的最多，近期seurat更新了V5的教程，听说主要是数据结构放生了改变，不过基本上换汤不换药，下次可以在讲V4的单细胞数据结构的时候穿插一下V5版本的改动，以及seuratV5流程中的一些改变，能让大家实现V4和V5的精准切换。同时，本篇教程中如有什么问题，也欢迎小伙伴们指出，我们可以进行讨论和修改！

不一定要知道的趣味小知识

不一定要知道 | Seurat 包为什么叫 Seurat！| 单细胞 x 艺术家 梦幻联动！

文末碎碎念

那今天的分享就到这里啦！我们下期再见哟！

最后顺便给自己推荐一下嘿嘿嘿！

如果我的分享对你有用的话，欢迎关注点赞在看转发分享阿巴阿巴阿巴阿巴巴巴！这可是我的第一原动力！

蟹蟹你们的喜欢和支持！！！