第七计 无中生有
本指本来没有却硬说有。现形容凭空捏造。

设置Seurat对象

对于本教程，我们将分析可从10X Genomics免费获得的外周血单个核细胞（PBMC）数据集。在Illumina NextSeq 500上已对2700个单细胞进行了测序。可以在此处找到原始数据。

我们从读取数据开始。该[Read10X()](https://satijalab.org/seurat/reference/Read10X.html)函数从10X读取cellranger管道的输出，返回唯一的分子识别（UMI）计数矩阵。该矩阵中的值表示在每个单元格（列）中检测到的每个特征（即基因；行）的分子数。

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

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

# Load the PBMC dataset
pbmc.data <- Read10X(data.dir = "../data/pbmc3k/filtered_gene_bc_matrices/hg19/")
# Initialize the Seurat object with the raw (non-normalized data).
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)

标准的预处理流程

以下步骤涵盖了Seurat中scRNA-seq数据的标准预处理工作流程。这些代表基于质量控制指标，数据归一化和缩放以及高度可变特征的检测来选择和过滤单元。

质控和选择细胞进行进一步分析

Seurat允许您轻松浏览质量控制指标并根据任何用户定义的标准过滤单元。社区常用的一些质量控制指标包括

• 每个细胞中检测到的独特基因的数量。
  • 低质量的细胞或空液滴通常很少有基因
  • 细胞双峰或多重峰可能显示异常高的基因计数
• 同样，细胞内检测到的分子总数（与独特基因高度相关）
• 映射到线粒体基因组的读段的百分比
  • 低质量/濒临死亡的细胞通常表现出广泛的线粒体污染
  • 我们使用[PercentageFeatureSet()](https://satijalab.org/seurat/reference/PercentageFeatureSet.html)函数计算线粒体QC指标，该函数计算源自一组功能的计数的百分比
  • 我们使用所有MT-以线粒体基因开头的基因

# The [[ operator can add columns to object metadata. This is a great place to stash QC stats
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

在下面的示例中，我们可视化QC指标，并使用它们来过滤单元格。

• 我们过滤具有超过2500或少于200的唯一特征计数的像元
• 我们过滤线粒体计数> 5％的细胞

# Visualize QC metrics as a violin plot
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

image

# FeatureScatter is typically used to visualize feature-feature relationships, but can be used
# for anything calculated by the object, i.e. columns in object metadata, PC scores etc.

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

image

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)

识别高度可变的特征（特征选择）

接下来，我们计算特征的子集，这些特征在数据集中显示出较高的单元间差异（即，它们在某些单元格中高表达，而在其他单元格中低表达）。我们和其他人发现，在下游分析中关注这些基因有助于在单细胞数据集中突出显示生物信号。

我们在Seurat中的过程在此进行了详细描述，并通过直接建模单细胞数据中固有的均值-方差关系对以前的版本进行了改进，并在[FindVariableFeatures()](https://satijalab.org/seurat/reference/FindVariableFeatures.html)函数中实现了该过程。默认情况下，我们为每个数据集返回2,000个要素。这些将用于下游分析，例如PCA。

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

# Identify the 10 most highly variable genes
top10 <- head(VariableFeatures(pbmc), 10)

# plot variable features with and without labels
plot1 <- VariableFeaturePlot(pbmc)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1 + plot2

image

---

缩放数据

接下来，我们应用线性变换（“缩放”），这是像PCA这样的降维技术之前的标准预处理步骤。该[ScaleData()](https://satijalab.org/seurat/reference/ScaleData.html)函数：

• 移动每个基因的表达，以使整个细胞的平均表达为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></details>

---

执行线性尺寸缩减

接下来，我们对缩放后的数据执行PCA。默认情况下，仅将先前确定的变量特征用作输入，但是features如果您希望选择其他子集，则可以使用arguments进行定义。

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

瑟拉提供可视化细胞和定义PCA，包括功能的几种有用的方法VizDimReduction()，[DimPlot()](https://satijalab.org/seurat/reference/DimPlot.html)和[DimHeatmap()](https://satijalab.org/seurat/reference/DimHeatmap.html)

# Examine and visualize PCA results a few different ways
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

DimPlot(pbmc, reduction = "pca")

image

特别是[DimHeatmap()](https://satijalab.org/seurat/reference/DimHeatmap.html)可以轻松探索数据集中异质性的主要来源，并且在尝试确定要包括哪些PC以便进行进一步的下游分析时非常有用。单元和要素均根据其PCA分数排序。设置cells为数字会在频谱的两端绘制“极端”单元，从而极大地加快了大型数据集的绘制速度。尽管显然是监督分析，但我们发现这是探索相关特征集的有价值的工具。

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

image

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

image

确定数据集的“维数”

为了克服scRNA-seq数据的任何单个特征中的广泛技术噪声，Seurat根据其PCA分数对细胞进行聚类，每个PC本质上代表一种“元特征”，该特征将跨相关特征集的信息进行组合。因此，最主要的主成分代表了数据集的稳健压缩。但是，我们应该选择包括多少个组件？10个？20吗 100？

在Macosko等人中，我们实施了受JackStraw程序启发的重采样测试。我们随机置换数据的一部分（默认为1％），然后重新运行PCA，以构建特征分数的“零分布”，然后重复此过程。我们将“重要”的PC识别为具有丰富的低p值功能的PC。

# NOTE: This process can take a long time for big datasets, comment out for expediency. More
# approximate techniques such as those implemented in ElbowPlot() can be used to reduce
# computation time
pbmc <- JackStraw(pbmc, num.replicate = 100)
pbmc <- ScoreJackStraw(pbmc, dims = 1:20)

该[JackStrawPlot()](https://satijalab.org/seurat/reference/JackStrawPlot.html)功能提供了一个可视化工具，用于比较每个PC的p值分布与均匀分布（虚线）。“重要”的PC将显示出具有丰富的低p值功能（虚线上方的实线）。在这种情况下，似乎在前10到12台PC之后，重要性显着下降。

JackStrawPlot(pbmc, dims = 1:15)

image

另一种启发式方法会生成“肘形图”：基于每个分量（[ElbowPlot()](https://satijalab.org/seurat/reference/ElbowPlot.html)函数）所解释的方差百分比对主成分进行排序。在此示例中，我们可以在PC9-10周围观察到一个“弯头”，这表明大多数真实信号是在前10台PC中捕获的。

ElbowPlot(pbmc)

image

识别数据集的真实维度–对用户而言可能是挑战/不确定。因此，我们建议考虑这三种方法。第一个受到更多监督，探索PC以确定异质性的相关来源，例如可以与GSEA结合使用。第二种实现基于随机空模型的统计测试，但是对于大型数据集而言非常耗时，并且可能无法返回清晰的PC截止值。第三种是常用的启发式方法，可以立即进行计算。在此示例中，所有三种方法都产生了相似的结果，但是我们可能有理由在PC 7-12之间选择任何东西作为临界值。

我们在这里选择10，但鼓励用户考虑以下内容：

• 树突状细胞和NK爱好者可能认识到与PC 12和13紧密相关的基因定义了罕见的免疫亚群（即MZB1是浆细胞样DC的标记）。但是，这些组非常罕见，在没有先验知识的情况下，很难将它们与背景噪声区分开来。
• 我们鼓励用户使用不同数量的PC（10台，15台，甚至50台！）重复进行下游分析。正如您将看到的，结果通常并没有太大的不同。
• 我们建议用户选择此参数时偏高一点。例如，仅用5台PC进行下游分析确实会对结果产生不利影响。

---

聚集细胞

Seurat v3在（Macosko等人）的初始策略的基础上，应用了基于图的聚类方法。重要的是，驱动聚类分析的距离度量（基于先前确定的PC）保持不变。但是，我们将蜂窝距离矩阵划分为群集的方法已得到了极大的改进。我们的方法受到最近手稿的启发，这些手稿将基于图的聚类方法应用于scRNA-seq数据[SNN-Cliq，Xu和Su，Bioinformatics，2015]和CyTOF数据[PhenoGraph，Levine等人，Cell，2015]。简而言之，这些方法将单元格嵌入图结构（例如，K近邻图（KNN）），并在具有相似特征表达模式的单元格之间绘制边，然后尝试将该图划分为高度互连的“准斜体”或“社区”。

像在PhenoGraph中一样，我们首先基于PCA空间中的欧式距离构造一个KNN图，然后基于两个相邻像元在其局部邻域中的共享重叠来细化任意两个像元之间的边缘权重（Jaccard相似性）。使用该[FindNeighbors()](https://satijalab.org/seurat/reference/FindNeighbors.html)功能执行此步骤，并将先前定义的数据集维度（前10个PC）作为输入。

为了对单元进行聚类，我们接下来将应用模块化优化技术（例如Louvain算法（默认）或SLM [SLM，Blondel等，统计力学杂志]）将单元迭代地分组在一起，以优化标准模块化功能。该[FindClusters()](https://satijalab.org/seurat/reference/FindClusters.html)函数实现此过程，并包含一个分辨率参数，该参数设置下游集群的“粒度”，其值越大，导致集群的数量越多。我们发现，将此参数设置在0.4-1.2之间通常可为大约3K个单元的单单元数据集返回良好的结果。对于较大的数据集，最佳分辨率通常会增加。可以使用[Idents()](https://satijalab.org/seurat/reference/Idents.html)函数找到簇。

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

# Look at cluster IDs of the first 5 cells
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的输入，我们建议使用相同的PC作为聚类分析的输入。

# If you haven't installed UMAP, you can do so via reticulate::py_install(packages =
# 'umap-learn')
pbmc <- RunUMAP(pbmc, dims = 1:10)

# note that you can set `label = TRUE` or use the LabelClusters function to help label
# individual clusters
DimPlot(pbmc, reduction = "umap")

image

您可以在此时保存该对象，以便可以轻松地将其重新加载，而不必重新运行上面执行的计算量大的步骤，也可以轻松地与协作者共享该对象。

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

---

查找差异表达的特征（簇生物标志物）

Seurat可以帮助您找到通过差异表达定义聚类的标记。默认情况下，ident.1与所有其他单元格相比，它识别单个群集（在中指定）的正向和负向标记。[FindAllMarkers()](https://satijalab.org/seurat/reference/FindAllMarkers.html)自动执行所有集群的此过程，但是您也可以测试集群组之间的相互关系，或针对所有单元进行测试。

该min.pct参数要求在两组像元的任何一组中以最小百分比检测特征，而thresh.test参数要求在两组像元之间平均表达（平均）一定程度的特征。您可以将它们都设置为0，但是时间会大大增加-因为这将测试可能不太具有高度歧视性的大量功能。作为加快这些计算速度的另一种选择，max.cells.per.ident可以设置。这将降低每个身份类的采样率，使其没有更多的单元格。尽管通常会断电，但速度的提高可能会非常明显，而差异化程度最高的功能可能仍会上升到最高位置。

# find all markers of cluster 1
cluster1.markers <- FindMarkers(pbmc, ident.1 = 2, min.pct = 0.25)
head(cluster1.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

# 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_log2FC pct.1 pct.2     p_val_adj
## FCGR3A        2.150929e-209   4.267579 0.975 0.039 2.949784e-205
## IFITM3        6.103366e-199   3.877105 0.975 0.048 8.370156e-195
## CFD           8.891428e-198   3.411039 0.938 0.037 1.219370e-193
## CD68          2.374425e-194   3.014535 0.926 0.035 3.256286e-190
## RP11-290F20.3 9.308287e-191   2.722684 0.840 0.016 1.276538e-186

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

## # A tibble: 18 x 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.74e-109       1.07 0.897 0.593 2.39e-105 0       LDHB    
##  2 1.17e- 83       1.33 0.435 0.108 1.60e- 79 0       CCR7    
##  3 0\.              5.57 0.996 0.215 0\.        1       S100A9  
##  4 0\.              5.48 0.975 0.121 0\.        1       S100A8  
##  5 7.99e- 87       1.28 0.981 0.644 1.10e- 82 2       LTB     
##  6 2.61e- 59       1.24 0.424 0.111 3.58e- 55 2       AQP3    
##  7 0\.              4.31 0.936 0.041 0\.        3       CD79A   
##  8 9.48e-271       3.59 0.622 0.022 1.30e-266 3       TCL1A   
##  9 1.17e-178       2.97 0.957 0.241 1.60e-174 4       CCL5    
## 10 4.93e-169       3.01 0.595 0.056 6.76e-165 4       GZMK    
## 11 3.51e-184       3.31 0.975 0.134 4.82e-180 5       FCGR3A  
## 12 2.03e-125       3.09 1     0.315 2.78e-121 5       LST1    
## 13 1.05e-265       4.89 0.986 0.071 1.44e-261 6       GZMB    
## 14 6.82e-175       4.92 0.958 0.135 9.36e-171 6       GNLY    
## 15 1.48e-220       3.87 0.812 0.011 2.03e-216 7       FCER1A  
## 16 1.67e- 21       2.87 1     0.513 2.28e- 17 7       HLA-DPB1
## 17 7.73e-200       7.24 1     0.01  1.06e-195 8       PF4     
## 18 3.68e-110       8.58 1     0.024 5.05e-106 8       PPBP

Seurat可以使用test.use参数设置一些针对差异表达的测试（有关详细信息，请参见我们的DE小插图）。例如，ROC测试返回任何单个标记的“分类能力”（范围从0-随机，到1-完美）。

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

我们提供了几种可视化标记表达的工具。[VlnPlot()](https://satijalab.org/seurat/reference/VlnPlot.html)（显示整个集群中的表达概率分布），以及[FeaturePlot()](https://satijalab.org/seurat/reference/FeaturePlot.html)（在tSNE或PCA图上可视化特征表达）是我们最常用的可视化。我们还建议探索[RidgePlot()](https://satijalab.org/seurat/reference/RidgePlot.html)，[CellScatter()](https://satijalab.org/seurat/reference/CellScatter.html)和[DotPlot()](https://satijalab.org/seurat/reference/DotPlot.html)作为查看数据集的其他方法。

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

image

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

image

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

image

[DoHeatmap()](https://satijalab.org/seurat/reference/DoHeatmap.html)生成给定单元和特征的表达式热图。在这种情况下，我们将为每个聚类绘制前20个标记（如果小于20，则绘制所有标记）。

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

image

---

将单元类型标识分配给集群

幸运的是，对于该数据集，我们可以使用规范标记轻松地将无偏聚类与已知单元格类型进行匹配：

集群ID	标记物	电池类型
0	IL7R，CCR7	天真CD4 + T
1个	CD14，LYZ	CD14 +单声道
2个	IL7R，S100A4	记忆体CD4 +
3	MS4A1	乙
4	CD8A	CD8 + T
5	FCGR3A，MS4A7	FCGR3A +单声道
6	GNLY，NKG7	NK
7	FCER1A，CST3	直流电
8	PPBP	血小板

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

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