Compiled: May 23, 2023

设置Seurat对象

对于本教程，我们将分析来自10X Genomics的一个外周血单个核细胞（PBMC）数据集。该数据集包含2,700个单个细胞，采用Illumina NextSeq 500进行测序。原始数据可以在这里找到。 here.

我们首先读取数据。Read10X()函数读取10X的cellranger流程的输出结果，返回一个唯一分子标识（UMI）计数矩阵。该矩阵中的值表示在每个细胞（列）中检测到的每个特征（即基因；行）的分子数量。

接下来，我们使用计数矩阵创建一个Seurat对象。该对象作为一个容器，包含了单细胞数据集的数据（如计数矩阵）和分析结果（如PCA或聚类结果）。关于Seurat对象结构的技术讨论，请参阅我们的GitHub Wiki。例如，计数矩阵存储在pbmc[["RNA"]]@counts中。

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

# 加载PBMC数据集
pbmc.data <- Read10X(data.dir = "../data/pbmc3k_filtered_gene_bc_matrices/filtered_gene_bc_matrices/hg19/")
# 用原始（非归一化）数据初始化(创建)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)

<details style="box-sizing: border-box; display: block;"><summary style="box-sizing: border-box; display: list-item;">一个计数矩阵中的数据是怎样的？</summary></details>

标准的预处理流程

以下步骤包含Seurat中scRNA-seq数据的标准预处理工作流程。这些步骤涉及基于质量控制指标选择和过滤细胞、数据归一化和缩放，以及检测高度可变特征。

质控（Quality Control，QC）和选择进一步分析的细胞

Seurat允许您轻松地探索质控指标，并根据任何用户定义的标准筛选细胞。一些指控指标 常用的

• 每个细胞中检测到的唯一基因数量。
  • 低质量的细胞或空滴在通常情况下会检测到非常少的基因。
  • 细胞双胞体或多胞体可能会表现出异常高的基因计数。
• 类似地，细胞内检测到的分子总数（与唯一基因强相关）也会呈现出相似的趋势。
• 映射到线粒体基因组的读取比例。
  • 低质量/濒临死亡的细胞通常表现出严重的线粒体污染。
  • 我们使用PercentageFeatureSet()函数计算线粒体的质控指标，该函数计算来自一组特征的计数所占的百分比。
  • 我们使用以MT-开头的所有基因作为线粒体基因集。

# [[ 操作符可以向object metadata添加列。这是一个存储质控统计信息好地方。
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

<details style="box-sizing: border-box; display: block;"><summary style="box-sizing: border-box; display: list-item;">质控指标存储在Seurat中的哪里？</summary>

在下面的示例中，我们可视化质控指标，并使用它们来筛选细胞。

• 我们过滤掉具有唯一特征（基因）计数超过2,500或小于200的细胞。
• 我们过滤掉具有大于5%线粒体计数的细胞。

# 将质控指标可视化为小提琴图。
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

image.png

# FeatureScatter通常用于可视化特征之间的关系，但可以用于对象计算的任何内容，例如对象元数据中的列，PC分数等。

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)

对数据进行归一化处理。

在从数据集中去除不需要的细胞之后，下一步是对数据进行归一化处理。默认情况下，我们使用全局缩放归一化方法“LogNormalize”，该方法通过每个细胞的总表达量对特征表达测量值进行归一化，然后乘以一个比例因子（默认为10,000），并对结果进行对数转换。归一化后的值存储在pbmc[["RNA"]]@data中。

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

为了清晰起见，在这个先前的代码行（以及将来的命令）中，我们提供了函数调用中某些参数的默认值。然而，这不是必需的，可以通过以下方式实现相同的行为：

pbmc <- NormalizeData(pbmc)

高变异特征（基因）的识别（特征/基因选择）

接下来，我们计算数据集中呈现细胞间变异性较高的一部分特征/基因（即，在某些细胞中高表达，在其他细胞中低表达）。我们和[其他研究人员] (https://www.nature.com/articles/nmeth.2645)发现，在后续的分析中聚焦于这些基因有助于突出显示单细胞数据集中的生物学信号%E5%8F%91%E7%8E%B0%EF%BC%8C%E5%9C%A8%E5%90%8E%E7%BB%AD%E7%9A%84%E5%88%86%E6%9E%90%E4%B8%AD%E8%81%9A%E7%84%A6%E4%BA%8E%E8%BF%99%E4%BA%9B%E5%9F%BA%E5%9B%A0%E6%9C%89%E5%8A%A9%E4%BA%8E%E7%AA%81%E5%87%BA%E6%98%BE%E7%A4%BA%E5%8D%95%E7%BB%86%E8%83%9E%E6%95%B0%E6%8D%AE%E9%9B%86%E4%B8%AD%E7%9A%84%E7%94%9F%E7%89%A9%E5%AD%A6%E4%BF%A1%E5%8F%B7).

我们在Seurat中的处理过程在这里(https://doi.org/10.1016/j.cell.2019.05.031),中详细描述，并通过直接建模单细胞数据中固有的均值-方差关系来改进之前的版本。该处理过程实现在FindVariableFeatures()函数中。默认情况下，我们每个数据集返回2,000个特征。这些特征将用于后续的分析，比如PCA。

pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)

# 识别出最高变异的10个基因
top10 <- head(VariableFeatures(pbmc), 10)

# 绘制带有和不带有标签的变异特征（基因）图
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1 + plot2

image.png

对数据进标准化处理

接下来，我们对数据应用线性变换（“scaling”），这是在进行诸如PCA之类的降维技术之前的标准预处理步骤。ScaleData()函数：

• 将每个基因的表达值进行转换，使得所有细胞中基因表达的平均值为0。
• 对每个基因的表达值进行标准化，使得所有细胞中基因表达的方差为1。
  • 这一步在后续的分析中给予相等的权重，以确保高表达基因不会主导分析结果。
• 这些结果存储在 pbmc[["RNA"]]@scale.data 中。

all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)

<details style="box-sizing: border-box; display: block;"><summary style="box-sizing: border-box; display: list-item;">这个步骤太耗时，我能让它更快一下吗</summary></details> <details style="box-sizing: border-box; display: block;"><summary style="box-sizing: border-box; display: list-item;">如何去除不需要的变异源，就像在Seurat v2中一样？</summary>%E4%B8%AD%E8%BF%9B%E8%A1%8C%E4%BA%86%E6%8F%8F%E8%BF%B0%EF%BC%8C%E5%B9%B6%E5%9C%A8%E4%BD%BF%E7%94%A8Seurat)</details>

进行线性降维

接下来，我们对已经缩放的数据进行PCA（主成分分析）。默认情况下，只使用之前确定的高变基因作为输入，但如果您希望选择不同的子集，可以使用features参数进行定义。

pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))

Seurat提供了几种有用的方法来可视化PCA定义的细胞和特征，包括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")

image.png

DimPlot(pbmc, reduction = "pca")

image.png

特别是DimHeatmap()函数允许轻松探索数据集中的主要异质性来源，并且在决定进一步下游分析中应包含哪些主成分时非常有用。细胞和特征根据它们的PCA得分进行排序。通过将cells参数设置为一个数字，可以在谱系（话外音：spectrum在这里翻译成什么比较好呢？）的两端绘制’极端’细胞，这对于大型数据集的绘图速度大大加快。虽然明显是一种有监督的分析，但我们发现这是一种有价值的工具，可以探索相关特征集。

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

image.png

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

image.png

确定数据集的’维度性’

为了克服scRNA-seq数据中任何单个特征（基因）的大量技术噪音，Seurat根据它们的PCA得分对细胞进行聚类，每个PC基本上代表了跨相关特征集合的’元特征’的组合信息。因此，顶部主成分代表了数据集的稳健压缩。然而，我们应该选择包含多少个主成分呢？10个？20个？100个？

在Macosko et al中，我们实施了一种受JackStraw程序启发的重采样测试。我们随机排列数据的一个子集（默认为1%），重新运行PCA，构建特征得分的“零分布”，并重复此过程。我们将具有低p值特征的主成分视为“显著”主成分。

# 注意：对于大型数据集，这个过程可能需要很长时间，请注释掉以加快速度。可以使用更加近似的技术，例如ElbowPlot()中实现的技术，以减少计算时间。
pbmc <- JackStraw(pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(pbmc, dims = 1:20)

JackStrawPlot()函数提供了一种可视化工具，用于将每个主成分的p值分布与均匀分布（虚线）进行比较。’显著’的主成分将显示出具有低p值的特征的显著富集（实线曲线在虚线之上）。在这种情况下，似乎在前10-12个主成分之后，显著性迅速下降。

JackStrawPlot(pbmc, dims = 1:15)

image.png

另一种启发式方法生成一个“弯曲点图”：根据每个主成分解释的方差百分比对主成分进行排名（使用ElbowPlot()函数）。在这个例子中，我们可以观察到在PC9-10附近有一个“弯曲点”，这表明大部分真实信号被捕捉在前10个主成分中。

ElbowPlot(pbmc)

image.png

确定数据集的真实维度对用户来说可能是具有挑战性和不确定性的。因此，我们建议考虑以下三种方法。第一种方法更加监督，探索主成分以确定相关的异质性来源，可以与GSEA等方法结合使用。第二种方法基于随机空模型进行统计测试，但对于大型数据集来说计算时间较长，并且可能无法得出明确的主成分截断点。第三种方法是常用的启发式方法，可以立即计算得出。在这个例子中，这三种方法都得到了类似的结果，但我们可以选择在PC 7-12之间任何一个作为截断值。

我们在这里选择了10个主成分，但鼓励用户考虑以下几点：

• 擅长树突状细胞和自然杀伤细胞的专家可能会注意到，与第12和第13个主成分强相关的基因定义了罕见的免疫亚群（例如，MZB1是浆细胞样树突状细胞的标记物）。然而，这些亚群在这样规模的数据集中非常罕见，如果没有先前的知识，很难将它们与背景噪声区分开来。
• 我们鼓励用户使用不同数量的主成分（例如10个、15个，甚至50个）来重复进行后续分析。正如您所观察到的，结果通常并没有显著差异。
• 我们建议用户在选择该参数时保守一些，选择较高的值。例如，仅使用5个主成分进行后续分析会显着且不利地影响结果。

对细胞进行聚类

Seurat v3采用了基于图的聚类方法，构建在Macosko等人(http://www.cell.com/abstract/S0092-8674(15)00549-8))的研究的初始策略基础上。重要的是，驱动聚类分析的距离度量（基于先前确定的主成分）保持不变。然而，我们在将细胞距离矩阵分割成聚类的方法上有了显著改进。我们的方法受到了最近应用图形聚类方法于scRNA-seq数据SNN-Cliq, Xu和Su, Bioinformatics, 2015和CyTOF数据PhenoGraph, Levine等人, Cell, 2015的启发。简而言之，这些方法将细胞嵌入到图结构中 - 例如K最近邻（KNN）图，其中边连接具有相似特征表达模式的细胞，然后试图将该图分割成高度相互连接的“准团簇”或“社区”。

与PhenoGraph类似，我们首先根据PCA空间中的欧氏距离构建K最近邻（KNN）图，并根据两个细胞在局部邻域中的共享重叠（Jaccard相似度）来调整它们之间的边权重。这一步使用FindNeighbors()函数执行，并将之前定义的数据集的维度（前10个主成分）作为输入。

为了对细胞进行聚类，我们接下来应用模块度优化技术，如Louvain算法（默认）或SLM SLM, Blondel et al., Journal of Statistical Mechanics，以迭代地将细胞分组，以优化标准模块度函数。FindClusters()函数实现了这个过程，并包含一个解析度参数，用于设置下游聚类的“粒度”，增加值会导致更多的聚类。我们发现，在约3K个细胞的单细胞数据集中，将该参数设置在0.4-1.2之间通常会得到良好的结果。对于较大的数据集，最佳解析度通常会增加。可以使用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: 95965
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8723
## Number of communities: 9
## Elapsed time: 0 seconds

# 查看前5个细胞的聚类ID
head(Idents(pbmc), 5)

## AAACATACAACCAC-1 AAACATTGAGCTAC-1 AAACATTGATCAGC-1 AAACCGTGCTTCCG-1 
##                2                3                2                1 
## AAACCGTGTATGCG-1 
##                6 
## Levels: 0 1 2 3 4 5 6 7 8

运行非线性降维方法（UMAP或tSNE）

Seurat提供了几种非线性降维技术，例如tSNE和UMAP，用于可视化和探索这些数据集。这些算法的目标是学习数据的潜在流形，以便将相似的细胞放置在低维空间中的一起。在上述基于图的聚类中确定的细胞应该在这些降维图中共同定位。作为UMAP和tSNE的输入，我们建议使用与聚类分析相同的主成分。

# 如果您尚未安装UMAP，您可以通过reticulate::py_install(packages = 'umap-learn')来安装它。
pbmc <- RunUMAP(pbmc, dims = 1:10)

# 请注意，您可以设置label = TRUE或使用LabelClusters函数来帮助标记单个聚类。
DimPlot(pbmc, reduction = "umap")

image.png

您可以在此时保存对象，以便可以轻松地加载它，而无需重新运行上述计算密集型步骤（ intensive steps ），或者方便地与合作者共享。.

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

寻找差异表达的特征（簇标记物）

Seurat可以通过差异表达来帮助您找到定义簇的标记物。默认情况下，它会将一个簇（在ident.1中指定）与其他所有细胞进行比较，找出正向和负向的标记物。FindAllMarkers()函数可以自动化这个过程，适用于所有簇，但您也可以将一组簇与另一组簇进行比较，或将其与所有细胞进行比较。

min.pct参数要求在两组细胞中至少以一定百分比检测到某个特征，而thresh.test参数要求某个特征在两组细胞之间（平均而言）有一定程度的差异表达。您可以将这两个参数都设置为0，但这会显著增加计算时间，因为这将测试大量不太可能具有高度区分性的特征。为了加快计算速度，您还可以设置max.cells.per.ident参数。这将使每个标识类的细胞数不超过设置的数量。虽然通常会损失一定的功效，但速度的提升可能是显著的，并且最高差异表达的特征可能仍然能够排在前面。

# 寻找簇2的所有标记物
cluster2.markers <- FindMarkers(pbmc, ident.1 = 2, min.pct = 0.25)
head(cluster2.markers, n = 5)

# 寻找将簇5与簇0和簇3区分开的所有标记物。
cluster5.markers <- FindMarkers(pbmc, ident.1 = 5, ident.2 = c(0, 3), min.pct = 0.25)
head(cluster5.markers, n = 5)

# 对每个聚类与所有其他细胞进行比较，找出只有正向标记物的结果。
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
pbmc.markers %>%
    group_by(cluster) %>%
    slice_max(n = 2, order_by = avg_log2FC)

Seurat提供了几种用于差异表达的测试方法，可以使用test.use参数进行设置（请参阅我们的DE vignette获取详细信息）。例如，ROC测试可以返回任何单个标记物的“分类能力”（范围从0 - 随机到1 - 完美）。

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

我们提供了几种用于可视化标记物表达的工具。VlnPlot()（显示不同聚类间的表达概率分布）和FeaturePlot()（在tSNE或PCA图上可视化特征表达）是我们最常用的可视化方法。我们还建议尝试使用RidgePlot()、CellScatter()和DotPlot()这些附加方法来查看数据集。

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

image.png

# 您也可以绘制原始counts。
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

DoHeatmap()会为给定的细胞和特征生成一个表达热图。在这种情况下，我们将为每个簇绘制前20个标记物（如果少于20个，则为所有标记物）。

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

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