Single-cell Workshop 2021 - 01 - scRNASeq 分析基础

学习资料：https://holab-hku.github.io/Fundamental-scRNA/

准备

软件包

载入下面的包，如果没有就先安装下。

library(Seurat)

Attaching SeuratObject

library(ggplot2)
library(patchwork)
library(dplyr)

Attaching package: 'dplyr'

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

library(hdf5r)

数据

下载数据：

Data used in this material is a 10k PBMC data getting from 10x Genomics website

https://github.com/holab-hku/Fundamental-scRNA/blob/master/data/10k_PBMC.h5

单细胞技术介绍

技术流程

单细胞RNA测序（scRNA-seq）包括一系列技术，以产生许多单个细胞的全基因组表达数据。

细胞分离

根据技术主要可以分为两类：Droplet-based 和 Non droplet-based。

Droplet-based：

• 组织样品必须解离成悬浮液
• 细胞将被单独封装成一个油包水滴
• 高通量，低成本
• 相关技术：Drop-seq (Macosko et al. 2015), inDrop (Klein et al. 2015), Chromium 10X (Zheng et al. 2017)

Non droplet-based：

• Smart-seq2 (Ramsköld et al. 2012)：用微型毛细管移液器手工细胞取样
• CEL-seq (Hashimshony et al. 2012)：单个细胞被添加到试管中；第一个介绍了条形码和RNA聚合
• MARS-seq (Jaitin et al. 2014)是第一个使用FACS将单细胞分离到单孔中的方法，优化版MARS-seq2 (Keren-Shaul et al. 2019)的推出成本更低，可重复性更高，井间污染更少

下图显示了一些常见的单细胞分离技术（Hwang, Lee, and Bang 2018）

Single Cell Isolation (modified from Hwang, Lee, and Bang 2018)

条形码和唯一分子标识符（UMI）

双端测序输出两个 fastq 文件，分别对应测序的5′和3′方向。使用这种测序技术，配对的第一个read总是与引物的细胞（条形码+UMI）部分一致。

Biased paired-end reads (David Tse et.al)

根据获得的细胞条形码、UMI 和 cDNA 的 reads，我们可以估计转录物的丰度。这允许 mapping 算法区分哪些序列是条形码，哪些是转录序列。因此，在确定细胞条形码和 UMI 条形码序列的长度和位置时，识别用于测序的文库制备化学方法非常重要。

为了获得 UMIs 的计数，我们可以首先通过细胞条形码对 reads 进行分组，然后对齐 cDNA reads 并使用 UMIs 对每个细胞每个基因的独特分子进行计数。

Grouping barcodes to assign reads to cells (modified from David Tse et.al)

对细胞条形码和 UMIs 的分析包括在校准过程中，我们将在下一节中介绍更多内容。

广泛应用的scRNA-seq技术汇总

Methods	Transcript coverage	UMI possibility	Strand specific	References
Tang method	Nearly full-length	No	No	Tang et al. (2009)
Quartz-Seq	Full-length	No	No	Sasagawa et al. (2013)
SUPeR-seq	Full-length	No	No	X. Fan et al. (2015)
Smart-seq	Full-length	No	No	Ramsköld et al. (2012)
Smart-seq2	Full-length	No	No	Picelli et al. (2013)
MATQ-seq	Full-length	Yes	Yes	Sheng et al. (2017)
STRT-seq STRT/C1	5′-only	Yes	Yes	Islam et al. (2011)
CEL-seq	3′-only	Yes	Yes	Hashimshony et al. (2012)
CEL-seq2	3′-only	Yes	Yes	Hashimshony et al. (2016)
MARS-seq	3′-only	Yes	Yes	Jaitin et al. (2014)
CytoSeq	3′-only	Yes	Yes	H. C. Fan, Fu, and Fodor (2015)
Drop-seq	3′-only	Yes	Yes	Macosko et al. (2015)
InDrop	3′-only	Yes	Yes	Klein et al. (2015)
Chromium	3′-only	Yes	Yes	Zheng et al. (2017)
SPLiT-seq	3′-only	Yes	Yes	Rosenberg et al. (2018)
sci-RNA-seq	3′-only	Yes	Yes	Cao et al. (2017)
Seq-Well	3′-only	Yes	Yes	Gierahn et al. (2017)
DroNC-seq	3′-only	Yes	Yes	Habib et al. (2017)
Quartz-Seq2	3′-only	Yes	Yes	Sasagawa et al. (2018)

预处理和质量控制

FASTQ 文件

原始 RNA 测序数据可能在 FASTQ 文件中。它是一种基于文本的格式，用于存储由单字母代码表示的读序列。 FASTQ 文件中的序列以@符号开头的readID开始，然后是序列数据行，一个简单的加号+分隔符和碱基质量分数。

它以以下格式表示:

@ReadID
READ SEQUENCE
+
SEQUENCING QUALITY SCORES

一般来说，fastq 文件是使用质量控制工具（如 FastQC）进行预处理的。这将输出一系列评估序列 reads 质量的指标。

Cell Ranger

Cell Ranger 是一组分析管道，用于处理 Chromium 单细胞数据以对齐 reads，生成特征条形码矩阵，执行聚类和其他二次分析等等。 它帮助我们生成 RNA reads 计数矩阵，我们将在学习中使用。

一些概念：

• GEM 孔（以前称为 GEM 组）：来自单个 10x Chromium™ 芯片通道的分隔单元（凝胶颗粒悬浮液）集合。可以从一个 GEM 孔中获得一个或多个测序文库。
• 文库（或测序文库）：从单个 GEM 孔中制备的带有 10x 条形码的测序文库。借助特征条形码或 V(D)J 分析，可以从同一个 GEM 孔中创建多个文库。文库类型可能包括基因表达、抗体捕获、CRISPR 引导捕获、TCR 富集等。
• 测序 Run（或 Flowcell）: A flowcell containing data from one sequencing instrument run.（这个从英文直译上很难理解，通俗的说就是一次上机测序得到的数据流）

通用流程

从 Chromium 10X 管道获得的单细胞数据可以使用 cellranger 通过以下工作流程进行处理。

Cellranger Workflow (taken from the cellranger website)

输出文件

下面是我们将从 cellranger 获得的输出文件夹。outs 文件夹包含最终的管道输出文件，其中包括我们需要用于下游分析的内容。

Overview of the folder generate from cellranger

Overview of the outs folder

以上是我们可以在outs文件夹中找到的内容。它包含了一些测序数据的总结信息，注释的读取序列，以及我们通常工作的基因表达矩阵。下面是我们想看的一些重要的输出文件。

• Matrices
• Web Summary .html
• Secondary Analysis CSV
• BAM
• Molecule Info (h5)
• Loupe File (.cloupe)

矩阵

测序时，Chromium 10X 不仅对转录组进行测序，还对任何可能的分子进行测序。这就导致了背景条形码的存在。细胞相关条形码是cellranger 认为标记来自真实细胞的转录组而不是背景的条形码。

对于不同版本的cellranger，使用不同的算法来检测细胞相关条形码。一般的想法是，细胞的条形码应该比背景条形码有更多的转录本计数。更多信息请访问：https://kb.10xgenomics.com/hc/en-us/articles/115003480523-How-are-barcodes-classified-as-cell-associated-。

cellranger管道输出两种类型的特征条形码矩阵：

• 未经过滤的特征条形码矩阵存储在raw_feature_bc_matrix（1-1）文件夹下。它包含了来自已知良好条形码序列固定列表的每个条形码，其中至少有一个读取，这意味着它包括背景和与细胞相关的条形码。
• 经过过滤的特征条形码矩阵存储在filtered_feature_bc_matrix（1-2）文件夹下。它只包括已检测到的与细胞相关的条形码。

raw_feature_bc_matrix和filtered_feature_bc_matrix文件夹中都包含三个文件。

• matrix.mtx.gz文件以稀疏矩阵的形式存储 reads 计数，其中每一行表示scRNA-seq数据中的一个基因，每一列表示一个细胞。行索引的信息存储在features.tsv.gz文件中，而列索引的信息存储在barcodes.tsv.gz文件中。
• features.tsv.gz文件对应于行索引。在scRNA-seq数据中，“features”实际上指的是基因。features.tsv.gz文件包含三列：feature ID：参考GTF文件的注释字段中的gene_id，表示特征的ID。 name：参考GTF文件的注释字段中的gene_name，如果参考GTF中没有gene_name字段，则基因名称等同于基因ID。 type of feature：描述特征类型，可以是Gene Expression、Antibody Capture、CRISPR或CUSTOM之一。对于scRNA-seq数据，它将是Gene Expression（基因表达）。
• barcodes.tsv.gz对应于列索引，它包含了每个列的条形码。有关条形码序列格式的更多详细信息，请参考条形码BAM部分。

Web Summary .html

一个概要的HTML文件包含了摘要指标和自动化二次分析结果。如果在流程运行期间检测到问题，将在此页面上显示警报。

该HTML文件包括两个部分，SUMMARY（摘要）和ANALYSIS（分析）。您还可以点击每个仪表板顶部的“?”以获取有关每个指标的更多信息。

1. SUMMARY（摘要）指标描述了测序质量和检测到的细胞的各种特征。在页面顶部附近醒目地显示了检测到的细胞数量、每个细胞的平均reads数以及每个细胞检测到的基因的中位数。

Cells dashboard (modified from the cellranger website)

2. 在“Cells”仪表板下的“Barcode Rank Plot”显示了条形码计数的分布以及被推断与细胞相关的条形码。y轴表示每个条形码映射到的UMI计数，x轴表示低于该值的条形码数量。陡峭的下降表示细胞关联的条形码与与空分区相关的条形码之间有良好的分离。条形码可以通过其UMI计数或其RNA配置文件确定为与细胞相关，因此图表的某些区域可能同时包含细胞关联和背景关联的条形码。图表的颜色表示与细胞关联的条形码的局部密度。

Sequencing and Mapping dashboards (modified from the cellranger website)

其他用于评估的指标：

• Estimated Number of Cells: 500-10,000
• Mean Reads per Cell: 20,000 reads/cell minimum recommended reads per cell around 10,000
• Valid barcodes: greater than 75%
• Q30 bases in RNA read: ideally greater than 65%
• Reads mapped confidently to transcriptome: ideally greater than 30%
• Reads mapped antisense to gene: ideally smaller than 10%

3. 分析（ANALYSIS）指标包括以下自动化的二次分析：

• 降维分析：将细胞投影到二维空间（t-SNE），以展示它们之间的关系。
• 自动聚类分析：将具有相似表达特征的细胞分组在一起，形成聚类。
• 差异表达基因列表：列出在所选聚类之间表达差异显著的基因。
• 显示测序深度降低对观察到的文库复杂性的影响的图表。
• 显示测序深度降低对每个细胞检测到的平均基因数的影响的图表。

BAM

BAM文件以二进制压缩格式保存了关于 mapping reads 的信息。它由可选的头部部分和对齐部分组成。如果存在头部部分，它将通过第一列中的@与对齐部分区分，并位于对齐部分之前。

当解压缩成SAM文件时，信息以制表符分隔的表格形式存储，其中包含一些标准列和由 Cell Ranger 软件生成的特定列。Cell Ranger特定的列包含与 BAM 条形码、BAM 对齐和特征条形码相关的信息。标准列对应以下内容：

QNAME : read name (generally will include UMI barcode if applicable)
FLAG : number tag indicating the “type” of alignment, link to explanation of all possible “types”
RNAME : reference sequence name (i.e. chromosome read is mapped to).
POS : leftmost mapping position
MAPQ : Mapping quality
CIGAR : string indicating the matching/mismatching parts of the read (may include soft-clipping).
RNEXT : reference name of the mate/next read
PNEXT : POS for mate/next read
TLEN : Template length (length of reference region the read is mapped to)
SEQ : read sequence
QUAL : read quality

可以使用 SAMtools 查看 BAM 文件：

samtools view output.bam

为了进行基于 RNA 速率的轨迹分析，需要 bam 文件。

分子信息（h5）

它是一个 HDF5 文件，包含有效条形码和有效 UMI 的所有分子的每分子信息，并以高置信度分配给基因或特征条形码。该 HDF5 文件包含与观察到的分子相对应的数据，以及有关所使用的库和特征集的数据。这个文件的结构是：

HDF5 File Hierarchy(taken from the cellranger website)

二级分析 CSV 文件

包含自动辅助分析结果的几个 CSV 文件。它包含有关降维、t-SNE、UMAP、聚类和差分表达的信息。它也通过Web Summary.html文件在ANALYSIS度量中可视化。

Loupe 文件 (.cloupe)

Loupe Browser 是一款桌面应用程序，提供与10x Genomics解决方案的数据进行交互式可视化分析的功能。它可以帮助寻找感兴趣的细胞、发现重要基因、识别细胞类型、探索细胞亚结构、研究细胞亚型、集成基因表达和V(D)J分析，并共享结果。

通过 Loupe Browser 可以查看 Loupe 文件，该文件包含以下信息：

• 样本中细胞的基因表达信息。
• 细胞的各种基于基因表达的聚类信息，包括t-SNE和UMAP投影以及差异基因表达情况。
• 来自转录组参考的基因注释信息。

STARsolo

STARsolo（Kaminow，Yunusov和Dobin 2021）是一个专为液滴式单细胞RNA测序数据（例如10X Genomics Chromium系统）设计的分析工具，直接内嵌在STAR代码中。STARsolo的输入是一个FASTQ文件，它可以以与Cell Ranger几乎相同的格式输出基因计数，但速度约快10倍。

STARsolo 程序输出大量反映 reads 比对过程细节的文件。在这里，我们只讨论其中一些关键文件。

STARsolo Output

• BAM文件包含有关比对 reads 的信息，与通过Cell Ranger生成的BAM文件非常相似。
• 比对摘要文件Features.stats和Summary.csv包含有关基本比对细节的信息。这可以作为对比对过程进行简单初步质量控制检查的便捷工具。
• 特征矩阵文件matrix.mtx包含每个个体细胞中映射的基因计数信息。列名对应于每个个体细胞的条形码，行名对应于所有注释的基因。由于该数据的规模较大，它以稀疏矩阵的形式存储。
• 辅助文件 barcodes.tsv 和 features.tsv 提供了下游分析所需的额外元数据。这些文件连同矩阵文件在功能上扮演着与Cell Ranger输出中的 Matrices 部分相同的角色，并且在金标准的 scRNA-seq 数据分析软件 seurat 中进行分析时是必需的。

Doublets

双细胞是指虽然设计为由一个细胞生成，但却是由两个细胞生成的人工文库。通常这是由于细胞分选或捕获过程中的错误引起的。

可以使用几种实验策略来去除双细胞：

• 基于不同供体个体之间的自然遗传变异（Kang等人，2018）。
• 使用与不同寡核苷酸结合的抗体标记一组细胞（例如，来自一个样本的所有细胞）。在混合后，观察到具有不同寡核苷酸的文库被视为双细胞并移除（Bach等人，2017）。
• 仅基于表达谱的计算方法（如模拟）推断双细胞并进行去除。

使用 Seurat 进行分析

本节的内容是根据Seurat - Guided Clustering Tutorial进行了一些修改和适应。我们使用的数据是从10x Genomics网站获取的10k PBMC数据。

我们将学习如何读取10X测序数据并将其转换为Seurat对象，进行质控和选择用于进一步分析的细胞，对数据进行归一化处理，识别高可变特征（特征选择），对数据进行缩放，进行线性降维和可视化分析。

Seurat 对象

Seurat对象充当一个容器，其中包含单个单细胞数据集的数据（如计数矩阵）和分析（如PCA或聚类结果）。在使用 Seurat 分析scRNA-seq数据之前，我们可以先从这里对Seurat对象有一些基本的了解。

设置Seurat对象

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

我们从读入数据开始。Read10X_h5从10X CellRanger hdf5文件读取计数矩阵，返回唯一的分子识别（UMI）计数矩阵。该矩阵中的值表示每个特征（即基因；行）在每个细胞（列）中检测到的表达量。它可用于读取scATAC-seq和scRNA-seq矩阵。

接下来我们使用计数矩阵来创建一个Seurat对象。

# Load the PBMC dataset
pbmc.data <- Read10X_h5("/Users/wsx/Library/CloudStorage/OneDrive-shanghaitech.edu.cn/Public/data/10k_PBMC.h5")
# Initialize the Seurat object with the raw (non-normalized data).
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc10k", min.cells = 3, min.features = 200)
pbmc

An object of class Seurat 
22432 features across 10813 samples within 1 assay 
Active assay: RNA (22432 features, 0 variable features)

如果我们想直接使用 cellranger 管道从 10X 读取数据，我们可以使用 Read10X()。

标准预处理流程

以下步骤包含了 Seurat 中 scRNA-seq 数据的标准预处理工作流程。它们基于我们将从 Cell Ranger 或 STARsolo 输出中获得的 RNA reads 计数矩阵。标准的预处理工作流程代表了基于 QC 指标的细胞选择和过滤，数据归一化和缩放，以及高度可变特征的检测。

QC 和细胞过滤

Seurat 允许你根据任何用户定义的标准轻松探索 QC 指标和过滤细胞。社区常用的一些质量控制指标（Ilicic et al. 2016）包括：

• 在每个细胞中检测到的独特基因的数量
  • 低质量的细胞或空液滴通常只有很少的基因
  • 细胞双胞胎或多胞胎可能表现出异常高的基因计数
  • 同样，细胞内检测到的分子总数（与独特的基因密切相关）
• reads 占到线粒体基因组的百分比
  • 低质量/垂死细胞通常表现出广泛的线粒体污染
  • 我们使用PercentageFeatureSet()函数计算线粒体 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-")

在 CreateSeuratObject() 过程中自动计算唯一基因和总分子的数量。它们存储在对象元数据中。

# Show QC metrics for the first 5 cells in the control group
head(pbmc@meta.data, 5)

                   orig.ident nCount_RNA nFeature_RNA percent.mt
AAACCCAGTATATGGA-1    pbmc10k        860          350 44.1860465
AAACCCAGTATCGTAC-1    pbmc10k       1548          729  0.4521964
AAACCCAGTCGGTGAA-1    pbmc10k       6387         1827 10.4117739
AAACCCAGTTAGAAAC-1    pbmc10k      16664         3744  5.2808449
AAACCCAGTTATCTTC-1    pbmc10k       3352         1464 13.8424821

我们可以可视化 nFeature_RNA, nCount_RNA 和 percent.mt 作为 QC 指标。

与 Cell Ranger 输出一样，以下结果中的 feature 表示基因。 nFeature_RNA 是每个细胞中检测到的基因数量。 nCount_RNA 是细胞内检测到的分子总数。 下图中的每个点代表一个细胞。

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

我们接着使用点图来表示 3 者之间的关系。

# 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") + 
  theme(legend.position="none")
plot2 <- FeatureScatter(pbmc, feature1 = "nCount_RNA", feature2 = "nFeature_RNA") + 
  theme(legend.position="none")
plot1 + plot2

在这里，我们过滤掉具有独特特征计数（基因）超过5000或少于200的细胞。我们也过滤掉线粒体数量大于15%的细胞。

pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 5000 & percent.mt < 15)

然后再一次用提琴图查看过滤后的质量图：

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

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

标准化数据

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

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

在这里，我们从基因表达矩阵中采样 10,000 个reads计数，分别可视化标准化前后的基因表达分布（不包括零）。

# set seed and put two plots in one figure
set.seed(123)
par(mfrow=c(1,2))
# original expression distribution
raw_geneExp = as.vector(pbmc[['RNA']]@counts) %>% sample(10000)

Warning in .sparse2v(x): sparse->dense coercion: allocating vector of size 1.6
GiB

raw_geneExp = raw_geneExp[raw_geneExp != 0]
hist(raw_geneExp)
# expression distribution after normalization
logNorm_geneExp = as.vector(pbmc[['RNA']]@data) %>% sample(10000)

Warning in .sparse2v(x): sparse->dense coercion: allocating vector of size 1.6
GiB

logNorm_geneExp = logNorm_geneExp[logNorm_geneExp != 0]
hist(logNorm_geneExp)

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

接下来，我们计算数据集中表现出高细胞间差异的特征子集（即，它们在一些细胞中高表达，而在其他细胞中低表达）。研究表明(Brennecke et al. 2013)在下游分析中关注这些基因有助于突出单细胞数据集中的生物信号。

Seurat中的过程详细描述在文献 Stuart et al. 2019 中，FindVariableFeatures()函数中实现直接建模单细胞数据中固有的均值-方差关系来改进以前的版本。默认情况下，Seurat每个数据集返回2000个特征，这些特征将用于下游分析，如PCA。

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

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

# plot variable features with and without labels
plot1 <- VariableFeaturePlot(pbmc) + 
  theme(legend.position="top")
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE) + 
  theme(legend.position="none")

When using repel, set xnudge and ynudge to 0 for optimal results

plot1 + plot2

Warning: Transformation introduced infinite values in continuous x-axis
Transformation introduced infinite values in continuous x-axis

缩放数据

接下来，我们应用线性变换(“缩放”)，这是在PCA等降维技术之前的标准预处理步骤。函数ScaleData()：

• 移动每个基因的表达，使细胞间的平均表达为 0。
• 缩放每个基因的表达，使细胞间的方差为 1。这一步在下游分析中具有同等的权重，因此高表达基因不会占主导地位。
• 结果存储在 pbmc[["RNA"]]@scale.data。

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

在使用上述命令时，我们使用所有基因来缩放数据。缩放是 Seurat 工作流程中的重要步骤，但仅限于将用作 PCA 输入的基因。因此，ScaleData()中的默认仅对先前确定的可变特征（默认值为2000）执行缩放。这将使这一步更快。

在这种情况下，我们的PCA和聚类结果将不受影响。然而，Seurat 热图（用DoHeatmap()生成，如下图所示）需要对热图中的基因进行缩放，以确保高表达的基因不会在热图中占主导地位。我们在本教程中缩放了所有基因。

ScaleData()函数还允许我们从单个单细胞数据集中删除不需要的变异源。例如，我们可以通过回归移除细胞周期阶段或线粒体污染相关的异质性。这个特性可以通过指定vars.to.regress来实现。即：

# skip here
pbmc <- ScaleData(pbmc, vars.to.regress = "percent.mt")

但是，特别是对于想要使用此功能的高级用户，Seurat 建议使用他们新的标准化工作流 SCTransform()。该方法在 Seurat 论文(Hafemeister and Satija 2019)中进行了描述，并有一个单独的技术文档。与ScaleData() 一样，函数SCTransform()也包含一个vars.to.regress参数。

执行线性降维

接下来，我们对缩放后的数据执行 PCA。默认情况下，只使用先前确定的可变特征作为输入，但是如果你希望选择不同的子集，则可以使用 features 参数定义。

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

Seurat提供了几种有用的方法来可视化定义 PCA 的细胞和特征，包括VizDimLoadings()、DimPlot() 和 DimHeatmap()。

# Examine and visualize PCA results a few different ways
print(pbmc[["pca"]], dims = 1:5, nfeatures = 5)

PC_ 1 
Positive:  LTB, IL32, TRAC, CD3D, TRBC2 
Negative:  FCN1, FGL2, CST3, IFI30, TYMP 
PC_ 2 
Positive:  MS4A1, CD79A, BANK1, IGHM, NIBAN3 
Negative:  IL32, GZMM, CD3D, CD7, CD247 
PC_ 3 
Positive:  GZMB, CLIC3, NKG7, GNLY, KLRD1 
Negative:  CCR7, LEF1, TRABD2A, TCF7, LTB 
PC_ 4 
Positive:  CD79B, MS4A1, GNLY, CD79A, LINC00926 
Negative:  LILRA4, CLEC4C, SERPINF1, TPM2, SCT 
PC_ 5 
Positive:  CDKN1C, HES4, CTSL, TCF7L2, BATF3 
Negative:  S100A12, ITGAM, VCAN, CES1, MGST1 

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

DimPlot(pbmc, reduction = "pca")

特别是 `DimHeatmap() 允许轻松地探索数据集中异构的主要来源，并且在试图决定将哪些pc包括在进一步的下游分析中时非常有用。

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

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

为了克服scRNA-seq数据中任何单个特征中广泛的技术噪声，Seurat基于它们的PCA分数对细胞进行聚类，每个 PC 本质上代表一个“元特征”，该元特征结合了相关特征集中的信息。因此，顶部主成分代表了数据集的鲁棒压缩。在这里，我们选择前20 个 PC。

细胞聚类

Seurat v3采用了基于图的聚类方法，基于(Macosko等人，2015)中的初始策略进行改进。重要的是，驱动聚类分析的距离度量（基于先前确定的主成分）保持不变。然而，Seurat在将细胞距离矩阵划分为聚类的方法有了显著改进。这种方法受到最近一些研究的启发，这些研究将基于图的聚类方法应用于scRNA-seq数据（Xu和Su，2015）和CyTOF数据（Levine等人，2015）。简而言之，这些方法将细胞嵌入到一个图结构中，例如K最近邻（KNN）图，其中类似特征表达模式的细胞之间有边连接，然后尝试将该图划分为高度相互连接的“拟团”或“社区”。

就像PhenoGraph一样，首先根据PCA空间中的欧氏距离构建了一个K最近邻图，然后根据局部邻域中的重叠情况（Jaccard相似性）调整任意两个细胞之间的边权重。这个步骤是使用FindNeighbors()函数完成的，它的输入是之前定义的数据集的维度（前20个主成分）。

为了对细胞进行聚类，接下来我们应用模块度优化技术，如Louvain算法（默认）或SLM（Blondel等人，2008），以迭代地将细胞分组在一起，目标是优化标准的模块度函数。FindClusters()函数实现了这个过程，并包含一个分辨率参数，用于设置下游聚类的“粒度”，增大的值会导致更多的聚类。我们发现，在大约3K个细胞的单细胞数据集中，将此参数设置在 0.4-1.2 之间通常会得到良好的结果。对于更大的数据集，通常需要增加最佳分辨率。可以使用Idents()函数找到这些聚类。

pbmc <- FindNeighbors(pbmc, dims = 1:20, verbose = FALSE)
pbmc <- FindClusters(pbmc, resolution = 0.5, verbose = FALSE)

# Look at cluster IDs of the first 5 cells
head(Idents(pbmc), 5)

AAACCCAGTATCGTAC-1 AAACCCAGTCGGTGAA-1 AAACCCAGTTAGAAAC-1 AAACCCAGTTATCTTC-1 
                 2                  4                  2                  8 
AAACCCAGTTGCCGAC-1 
                 0 
Levels: 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

运行非线性降维 UMAP/tSNE

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

pbmc <- RunUMAP(pbmc, dims = 1:20, verbose = FALSE)

Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session

然后得到单细胞聚类结果的UMAP图。

DimPlot(pbmc, reduction = "umap")

类似地，tSNE。

pbmc <- RunTSNE(pbmc, dims = 1:20, verbose = FALSE)
DimPlot(pbmc, reduction = "tsne")

贴上聚类标签：

DimPlot(pbmc, reduction = "umap", label = TRUE)

plot <- DimPlot(object = pbmc)
LabelClusters(plot = plot, id = 'ident')

是时候保存我们的分析结果了：

saveRDS(pbmc, 
file = "/Users/wsx/Library/CloudStorage/OneDrive-shanghaitech.edu.cn/Public/data/pbmc_processed.rds")