【原】scRNA-seq聚类分析(一)

回顾

单细胞RNA-seq分析介绍
单细胞RNA-seq的设计和方法
从原始数据到计数矩阵
差异分析前的准备工作
scRNA-seq——读入数据详解
scRNA-seq——质量控制
为什么需要Normalization和PCA分析

学习目标

• 对每个样本执行归一化、方差估计和最大差异基因的识别

• 使用最大差异的基因跨条件整合细胞，以鉴定彼此最相似的细胞

clustering analysis: aligning cells across conditions

现在我们有了高质量的细胞，在将细胞聚类并确定不同的潜在细胞类型之前，我们需要执行一些步骤。我们的数据集包含来自两个不同条件（Control and Stimulated）的两个样本，因此整合这些样本有助于更好地进行比较。在此之前，我们需要归一化我们的基因表达值，并根据我们数据集中最大的变异来源跨条件排列我们的细胞。在本节中，我们将在聚类之前讨论并执行这些初始步骤。

目标  

• 为了准确地规范和缩放基因表达值，以解决测序深度和过度分散的计数值的差异

• 找出最有可能指示不同细胞类型的变异基因

• 跨条件排序相似细胞

挑战

• 检查并消除不必要的差异，避免下游人为原因导致的细胞聚集

• 排列相似细胞类型的细胞，避免在下游分析中因细胞类型偏向性的聚类

建议

• 在执行聚类之前，要很好地了解您对要呈现的细胞类型的期望。了解你所期望的是低复杂性的细胞类型还是高线粒体含量的细胞类型，以及这些细胞是否正在分化。

• Regress out  number of UMIs (default using sctransform), mitochondrial content, and cell cycle, if needed and appropriate for experiment, so not to drive clustering downstream

Clustering workflow

从差异中往往能得到更多的信息，但并不是所有的差异都是有用的信息。聚类分析的目的是在我们的数据集中保留应该定义我们的细胞类型的主要变异来源，同时限制由于不感兴趣的变异来源(测序深度、细胞周期差异、线粒体表达、批量效应等)而产生的变异。

然后，为了确定存在的细胞类型，我们将使用变化最大的基因进行聚类分析，以定义数据集中的主要变异来源。

分析流程如下：

• Satija Lab: Seurat v3 Guided Integration Tutorial

• Paul Hoffman: Cell-Cycle Scoring and Regression

执行以下步骤进行聚类：

1. 对每个样品进行Normalization、variance stabilization和regression of unwanted variation(例如，线粒体转录本丰度、细胞周期时相等)。

2. 使用共享的高度可变基因整合样本(可选，但如果细胞类型按样本/条件分离，则建议对齐来自不同样本/条件的细胞)

3. Clustering cells based on top PCs (metagenes)

4. 探索质量控制指标：确定clusters是否与UMI、基因、细胞周期、线粒体含量、样本等不平衡

5. 使用已知的细胞类型特异性基因标记搜索预期的细胞类型

Set-up

为了执行此分析，我们将主要使用Seurat软件包中提供的功能。因此，我们需要加载Seurat库，如果尚未加载tidyverse库，也需加载。创建脚本SCT_integration_analysis.R并加载库：

1# Single-cell RNA-seq analysis - clustering analysis
2
3# Load libraries
4library(Seurat)
5library(tidyverse)
6library(RCurl)
7library(cowplot)


为了执行分析，Seurat要求数据作为seurat对象存在。我们已经在质量控制课程（filtered_seurat）中创建了该对象，因此我们可以直接使用它。

Normalization, variance stabilization, and regression of unwanted variation for each sample

分析的第一步是将原始计数归一化，以解决每个样品每个细胞的测序深度差异。Seurat最近介绍了一种新的scRNA-seq数据 normalization and variance stabilization方法，称为sctransform。

sctransform方法使用regularized negative binomial model 对UMI计数建模以去除由于测序深度(每个细胞的总nUMIs)引起的变化，同时根据具有相似丰度的基因之间的信息集合来调整方差(类似于一些批量RNA-seq方法)。

模型的输出(残差)是测试的每个转录本的归一化表达水平。

Sctransform会自动回归序列深度（nUMI）；但是，通常在特定于数据集的数据中还有其他无意义的变化来源。例如，对于某些数据集，细胞周期阶段可能是显着变化的来源，而对于其他数据集则不是。在归因于细胞周期阶段的变化之前，您需要检查细胞周期阶段是否是数据变化的主要来源。

Cell cycle scoring

建议在执行sctransform方法之前检查细胞周期阶段。由于细胞之间的计数需要具有可比性，并且每个细胞具有不同的UMI总数，因此我们通过除以每个细胞的总计数并取自然对数进行粗略归一化。这种方法的准确性不如我们最终用于识别细胞簇的sctransform方法，但是它足以探索我们数据中的变异源。

1# Normalize the counts
2seurat_phase <- NormalizeData(filtered_seurat)


将数据针对测序深度进行了归一化后，我们就可以根据细胞的G2/M和S期标记的表达情况给每个细胞一个评分。

我们提供了人类细胞周期标志物列表供您下载。右键单击此链接，将“另存为…”直接保存到data目录中。但是，如果您不使用人类数据，我们将提供其他材料，详细介绍如何获取其他感兴趣生物的细胞周期标记。

1# Load cell cycle markers
2load("data/cycle.rda")
3
4# Score cells for cell cycle
5seurat_phase <- CellCycleScoring(seurat_phase, 
6                                 g2m.features = g2m_genes, 
7                                 s.features = s_genes)
8
9# View cell cycle scores and phases assigned to cells                                 
10View([e-mail])    


在对细胞周期进行评分后，我们可以使用主成分分析（PCA）来确定细胞周期是否是我们数据集中的一个主要变异来源。要执行PCA，我们需要首先选择最大变异的特征，然后对数据进行缩放。由于高表达的基因表现出最高的变异量，我们不希望我们的“高可变基因”只反映高表达，我们需要将数据缩放到随表达水平变化的规模。Seurat ScaleData()函数将按以下方式缩放数据：

• 调整每个基因的表达，使细胞间的平均表达量为0

• 缩放每个基因的表达以使细胞间的方差为1

1# Identify the most variable genes
2seurat_phase <- FindVariableFeatures(seurat_phase, 
3                     selection.method = "vst",
4                     nfeatures = 2000, 
5                     verbose = FALSE)
6
7# Scale the counts
8seurat_phase <- ScaleData(seurat_phase)


注意：对于selection.method和nfeatures参数，指定的值是默认设置。因此，您不一定需要在代码中包含这些内容。为了透明起见，我们将其包含在这里，并通知您正在使用的内容。

现在，我们可以执行PCA分析并绘制排名靠前的PCs：

1# Perform PCA
2seurat_phase <- RunPCA(seurat_phase)
3
4# Plot the PCA colored by cell cycle phase
5DimPlot(seurat_phase,
6        reduction = "pca",
7        group.by= "Phase",
8        split.by = "Phase")


由于细胞周期阶段的不同，我们没有看到很大的差异。根据这个图，我们不能得出细胞周期的变异。

SCTransform

现在我们可以使用SCTransform方法作为更精确的归一化方法，估计原始过滤数据的方差，并识别最可变的基因。默认情况下，sctransform会计算细胞测序深度或nUMI。

我们已经检查了细胞周期，并认为它不代表我们数据中的主要变异来源，但线粒体的表达是另一个可以很大程度上影响聚类的因素。通常，排除由于线粒体表达而引起的变异是有用的。但是，如果线粒体基因表达的差异代表可能有助于区分细胞簇的生物学现象，那么我们建议不要使线粒体表达regressing 。

我们可以使用'for循环’在每个样本上运行NormalizeData()、CellCycleScoring()和SCTransform()，并通过SCTransform()函数的vars.to.regress参数来指定 regress 线粒体表达。

在运行此for循环之前，我们知道输出可能会在内存方面生成较大的R对象/变量。如果我们有一个大型数据集，则可能需要使用以下代码调整R内允许的对象大小限制(默认值为500*1024^2=500Mb)：

1options(future.globals.maxSize = 4000 * 1024^2)


现在，对所有样本执行细胞周期评分和sctransform。这可能需要一些时间(约10分钟)：

1# Split seurat object by condition to perform cell cycle scoring and SCT on all samples
2split_seurat <- SplitObject(filtered_seurat, split.by = "sample")
3
4split_seurat <- split_seurat[c("ctrl", "stim")]
5
6for (i in 1:length(split_seurat)) {
7    split_seurat[[i]] <- NormalizeData(split_seurat[[i]], verbose = TRUE)
8    split_seurat[[i]] <- CellCycleScoring(split_seurat[[i]], g2m.features=g2m_genes, s.features=s_genes)
9    split_seurat[[i]] <- SCTransform(split_seurat[[i]], vars.to.regress = c("mitoRatio"))
10    }


注：默认情况下，SCTransform在归一化、调整方差、消除无用的变异源后，会根据残差方差对基因进行排序，并输出3000个变异最多的基因。如果数据集具有较大的细胞编号，则使用variable.feature ures.n参数，将此参数调整得更高可能是有益的。

注意，输出的最后一行指定“将默认分析设置为SCT”。我们可以查看存储在Seurat对象中的不同分析。

1# Check which assays are stored in objects
2split_seurat$ctrl@assays


现在我们可以看到，除了原始RNA计数外，我们的assays插槽中还有一个SCT组件。最具可变性的特征将是SCT分析中存储的唯一基因。在进行scRNA-seq分析时，我们将选择最合适的分析方法用于分析的不同步骤。如果我们在一个Seurat对象中同时对这两个条件执行归一化，并可视化细胞之间的相似性，我们将看到特定于条件的聚类：

在特定条件下的细胞聚类表明我们需要跨条件整合细胞。

注意：Seurat有一个关于如何在没有集成的情况下运行工作流程的vignette。该工作流与此工作流非常相似，但是样本在开始时不一定要拆分，也不一定要执行整合。

未完待续……

References

[1] Seurat v3 Guided Integration Tutorial: https:///seurat/v3.0/immune_alignment.html
[2] Cell-Cycle Scoring and Regression: http:///seurat/cell_cycle_vignette.html
[3] 此链接，: https://www./s/hus4mrkueh1tfpr/cycle.rda?dl=1
[4] 其他材料，: https://hbctraining./scRNA-seq/lessons/cell_cycle_scoring.html
[5] vignette: https:///seurat/v3.1/sctransform_vignette.html

---

注：以上内容来自哈佛大学生物信息中心(HBC)_的教学团队的生物信息学培训课程。原文链接：https://hbctraining./scRNA-seq/schedule/ 点击 “阅读原文” 可直达