scRNA-seq/06_SC_SCT_and_integration.md at master · hbctraining/scRNA-seq · GitHub
https://github.com/hbctraining/scRNA-seq/blob/master/lessons/06_SC_SCT_and_integration.md
学习目标:
对每个样本执行归一化,方差估计和最大变异基因的鉴定
使用最变异的基因跨条件整合细胞以鉴定彼此最相似的细胞
单细胞RNA-seq聚类分析:跨条件排列细胞
现在我们有了高质量的细胞,在将细胞聚类并确定不同的潜在细胞类型之前,我们需要执行一些步骤。我们的数据集包含来自两个不同条件(对照和刺激)的两个样本,因此将这些样本整合在一起以更好地进行比较可能会有所帮助。我们将需要归一化我们的基因表达值,并根据我们数据集中最大的变异来源跨条件排列细胞。在本课程中,我们将在聚类之前讨论并执行这些初始步骤。
目标:
为了准确地标准化和缩放基因表达值,以解决测序深度和过度分散的计数值的差异。
为了找出最变异基因可能成为指示呈现不同的细胞类型。
跨条件排列相似的细胞。
挑战:
检查并消除不必要的变化
对齐相似细胞类型的细胞
建议:
在执行群集之前,最好对将要出现的细胞类型有一个预估。了解您是否预测低复杂度或更高线粒体含量的细胞类型,以及细胞是否在分化
如果需要并且适合进行实验,则淘汰UMI的数量(默认使用sctransform),线粒体含量和细胞周期,以免导致下游聚集
集群工作流程
为了使某些信息具有信息性,它需要表现出差异性,但并非所有的差异性都具有信息性。聚类分析的目的是在应定义细胞类型的数据集中保留主要的变异来源,同时由于无用的变异来源(测序深度,细胞周期差异,线粒体表达,批次效应等)而限制变异。)。
然后,为了确定当前的细胞类型,我们将使用变化最大的基因进行聚类分析,以定义数据集中变异的主要来源。
此分析的工作流程来自以下来源:
Satija Lab:Seurat v3指导的集成教程
保罗·霍夫曼(Paul Hoffman):细胞周期评分和回归
要识别集群,将执行以下步骤:
每个样品的归一化,方差稳定和不需要的变异(例如线粒体转录本丰度,细胞周期阶段等)的回归
使用共享的高度可变的基因对样品进行整合(可选,但建议按细胞/样品/条件分离细胞类型以比对不同样品/条件的细胞)
基于高变基因对细胞进行聚类
探索质量控制指标:确定簇是否与UMI,基因,细胞周期,线粒体含量,样品等不平衡。
使用已知的细胞类型特异性基因标记物搜索预期的细胞类型
在本课程中,我们将介绍集群工作流程的前两个步骤。
设置
为了执行此分析,我们将主要使用Seurat软件包中提供的功能。因此,如果尚未加载tidyverse库,我们还需要加载Seurat库。创建脚本SCT_integration_analysis.R
并加载库:
# 加载库
library(Seurat)
library(tidyverse)
library(RCurl)
library(cowplot)
为了执行分析,Seurat要求数据作为seurat
对象存在。我们已经在质量控制课程(filtered_seurat
)中创建了该对象,因此我们可以使用它。
每个样本的归一化,方差稳定和不想要的变化的回归
分析的第一步是将原始计数标准化,以解决每个样品每个细胞的测序深度差异。Seurat最近推出了一种称为sctransform的用于对scRNA-seq数据进行规范化和方差稳定的新方法。
sctransform方法使用正则化的负二项式模型对UMI计数进行建模,以消除由于测序深度(每个细胞的总nUMI)所引起的变异,同时基于具有相似丰度的基因之间的合并信息来调整方差(类似于某些批量RNA测序方法) )。
图片来源: Hafemeister C和Satija R.使用正则化负二项式回归对单细胞RNA-seq数据进行归一化和方差稳定化,bioRxiv 2019(https://doi.org/10.1101/576827)
模型的输出(残差)是测试的每个转录本的标准化表达水平。
Sctransform会自动回归序列深度(nUMI);但是,数据中还有其他有趣的变化来源,这些变化通常是特定于数据集的。例如,对于某些数据集,细胞周期阶段可能是显著变化的来源,而对于其他数据集则不是。在归因于细胞周期阶段的变化之前,需要检查细胞周期阶段是否是数据变化的主要来源。
细胞周期评分
建议在执行sctransform方法之前检查细胞周期阶段。由于每个细胞之间的计数必须是可比较的,并且每个细胞具有不同数量的总UMI,因此我们通过除以每个细胞的总计数并取自然对数来进行粗略的归一化。这种方法不如我们最终将用于识别细胞簇的sctransform方法准确,但是足以探索我们数据的变异源。
# 归一化计数
seurat_phase <- NormalizeData(filtered_seurat)
将数据的测序深度归一化后,我们可以根据其细胞G2 / M和S期标志物的表达为其分配得分。
我们提供了人类细胞周期标志物列表供您下载。右键单击此链接,将“另存为...”直接保存到您的data
目录中。但是,如果您不使用人类数据,我们还会提供其他材料,详细介绍如何获取其他感兴趣生物的细胞周期标记。
#加载细胞周期基因
load("data/cycle.rda")
# Score cells for cell cycle
seurat_phase <- CellCycleScoring(seurat_phase, g2m.features = g2m_genes, s.features = s_genes)
# 查看分配给细胞的细胞周期得分和相位
View(seurat_phase@meta.data)
在对细胞进行细胞周期评分之后,我们要使用PCA确定细胞周期是否是我们数据集中变异的主要来源。要执行PCA,我们需要首先选择变化最大的功能,然后缩放数据。由于高表达的基因表现出最高的变异量,并且我们不希望我们的“高度可变的基因”仅反映高表达,因此我们需要对数据进行缩放以随表达水平对变异进行缩放。SeuratScaleData()
函数将按以下方式缩放数据:
调整每个基因的表达以使整个细胞的平均表达为0
缩放每个基因的表达以使细胞间的差异为1
# 识别变化最大的基因
seurat_phase <- FindVariableFeatures(seurat_phase,
selection.method = "vst",
nfeatures = 2000,
verbose = FALSE)
#缩放计数
seurat_phase <- ScaleData(seurat_phase)
注意:对于
selection.method
和nfeatures
参数,指定的值是默认设置。因此,您不一定需要在代码中包含这些内容。我们将其包括在此处,并告知您所使用的内容。
现在,我们可以执行PCA分析并绘制顶级PC:
# 执行PCA
seurat_phase <- RunPCA(seurat_phase)
# 绘制以细胞周期
DimPlot(seurat_phase,
reduction = "pca",
group.by= "Phase",
split.by = "Phase")
由于细胞周期阶段,我们没有看到较大的差异。基于此图,我们不会回归出由于细胞周期而引起的变化。
SCT标准化
现在,我们可以将sctransform方法用作标准化,估计原始过滤数据的方差并确定变化最大的基因的更准确方法。默认情况下,sctransform会说明细胞测序深度或nUMI。
我们已经检查了细胞周期,并决定它不代表我们数据的主要变异来源,但是线粒体表达是另一个可以极大影响聚类的因素。通常,将由于线粒体表达而引起的变异消除非常有用。但是,如果线粒体基因表达的差异代表可能有助于区分细胞簇的生物学现象,那么我们建议不要使线粒体表达消除。
我们可以使用一个“for循环”来运行NormalizeData()
,CellCycleScoring()
以及SCTransform()
在每个样品上,并回归出线粒体表达通过在指定vars.to.regress
的参数的SCTransform()
函数。
在运行此代码之前for loop
,我们知道输出可以在内存方面生成较大的R对象/变量。如果我们有一个很大的数据集,则可能需要使用以下代码在R(默认值为500 * 1024 ^ 2 = 500 Mb)内调整允许的对象大小的限制:
选项(future.globals.maxSize = 4000 * 1024 ^ 2)
现在,要对所有样本执行细胞周期评分和sctransform。这可能需要一些时间(约10分钟):
# 按条件拆分seurat对象,以对所有样本执行细胞周期评分和SCT
split_seurat <- SplitObject(filtered_seurat, split.by = "sample")
split_seurat <- split_seurat[c("ctrl", "stim")]
for (i in 1:length(split_seurat)) {
split_seurat[[i]] <- NormalizeData(split_seurat[[i]], verbose = TRUE)
split_seurat[[i]] <- CellCycleScoring(split_seurat[[i]], g2m.features=g2m_genes, s.features=s_genes)
split_seurat[[i]] <- SCTransform(split_seurat[[i]], vars.to.regress = c("mitoRatio"))
}
注意:默认情况下,在标准化,调整方差并回归出无用的变异源之后,SCTransform将按残差对基因进行排序,并输出3000个最变异的基因。如果数据集具有较大的细胞编号,则使用
variable.features.n
参数将此参数调整得更高可能是有益的。
注意,输出的最后一行指定“将默认检测设置为SCT”。我们可以查看存储在seurat对象中的不同测定。
# 检查哪些分析存储在对象中
split_seurat$ctrl@assays
现在我们可以看到,除了原始RNA计数之外,我们的assays
插槽中现在还有一个SCT组件。最具可变性的特征将是SCT分析中存储的唯一基因。随着我们对scRNA-seq分析的进行,我们将选择最合适的分析方法用于分析中的不同步骤。
通常,在决定是否需要执行集成之前,我们总是先查看细胞。如果我们在Seurat对象中同时对两个条件进行了归一化并可视化了细胞之间的相似性,我们将看到特定于条件的聚类:
细胞的特定于条件的聚类表明我们需要跨条件整合细胞。
注意: Seurat有一个小插图,说明如何在不集成的情况下运行工作流。该工作流程与本工作流程非常相似,但是样本不一定会在一开始就被拆分,也不会执行集成。
使用共享的高度可变的基因整合样品
如果细胞按样品,条件,数据集或模式进行聚类,则此步骤可以极大地改善聚类和下游分析。如果不确定要期望的簇或期望条件之间的某些不同细胞类型(例如肿瘤和对照样品),则可以先单独运行条件,然后一起运行以查看两种条件下是否存在针对特定细胞类型的条件特定簇。通常,当根据多个条件对细胞进行聚类时,会存在特定于条件的聚类,而集成可以帮助确保同一细胞类型一起聚类。
为了进行整合,我们将使用通过SCTransform确定的每个条件使用共享的高度可变的基因,然后,我们将对这些条件进行“整合”或“协调”,以覆盖相似或在组之间具有“共同的生物学特征集”的细胞。这些组可以代表:
不同的处理条件(例如对照和刺激)
不同的数据集(例如,使用相同样品的不同文库制备方法生成的数据集的scRNA-seq)
不同的方式(例如,scRNA-seq和scATAC-seq)
整合是一种强大的方法,它利用这些变化最大的共享源来识别条件或数据集之间的共享亚群[ Stuart和Bulter等。(2018) ]。整合的目的是确保一个条件/数据集的细胞类型与其他条件/数据集的相同细胞类型对齐(例如,对照巨噬细胞与受刺激的巨噬细胞对齐)。
特别地,这种整合方法期望在组中的单个细胞的至少一个子集之间具有“对应”或共享的生物学状态。下图概述了集成分析中的步骤:
图片来源: Stuart T和Butler A等。单细胞数据的全面整合,bioRxiv 2018(https://doi.org/10.1101/460147)
应用的不同步骤如下:
执行规范相关分析(CCA):
CCA识别条件/组之间共享的变异源。它是PCA的一种形式,它可以识别数据中最大的变异源,但前提是要在条件/组之间共享或保存(使用每个样本中的3000个变异最大的基因)。
此步骤使用最大的共享变异源粗略对齐细胞。
注意:使用共享的高度可变的基因是因为它们最有可能代表那些区分存在的不同细胞类型的基因。
在数据集中识别锚点或相互最近的邻居(MNN)(有时会识别出不正确的锚点):可以将MNN视为“最佳伙伴”。对于处于一种情况的每个细胞:
根据基因表达值确定细胞在其他情况下最接近的邻居-它是“最佳伙伴”。
进行reciprical分析,并且如果两个群在两个方向上是“死党”,那么这些细胞将被标记为锚为“锚定”的两个数据集在一起。
“ MNN对中细胞之间表达值的差异提供了批处理效果的估计,通过对许多这样的对进行平均,可以使批处理效果更加精确。获得校正向量并将其应用于表达值以执行批处理校正。” [ Stuart和Bulter等。(2018) ]。
过滤锚点以删除不正确的锚点:
通过其本地邻域中的重叠来评估锚点对之间的相似性(不正确的锚点得分较低)-相邻细胞是否具有彼此相邻的“最佳伙伴”?整合条件/数据集:
使用锚点和相应的分数来转换细胞表达值,从而可以整合条件/数据集(不同的样本,条件,数据集,模态)
注意:每个细胞的转换使用每个锚点的两个细胞在数据集锚点之间的加权平均值。权重由细胞相似度评分(细胞与k个最接近的锚点之间的距离)和锚点评分确定,因此同一邻域中的细胞应具有相似的校正值。
如果一个数据集中存在细胞类型,而另一个数据集中则没有,则这些细胞仍将显示为单独的特定于样本的群集。
现在,使用我们的SCTransform对象作为输入,让我们跨条件执行集成。
首先,我们需要指定我们要使用由SCTransform识别的3000个可变性最高的基因进行整合。默认情况下,此功能仅选择前2000个基因。
# 选择要用于集成的锚点
integ_features <- SelectIntegrationFeatures(object.list = split_seurat,
nfeatures = 3000)
然后,我们需要准备SCTransform对象以进行集成。
#准备用于集成的SCT列表对象
# 准备用于集成的SCT列表对象
split_seurat <- PrepSCTIntegration(object.list = split_seurat,
现在,我们将执行CCA,找到最好的伙伴或锚点,并过滤不正确的锚点。对于我们的数据集,这最多需要15分钟才能运行。另外,请注意,控制台中的进度条将保持为0%,但要知道它实际上正在运行。
# 找最佳锚点-需要运行一段时间
integ_anchors <- FindIntegrationAnchors(object.list = split_seurat, normalization.method = "SCT",anchor.features = integ_features)
最后,我们可以跨条件进行整合。
# 跨条件整合
seurat_integrated <- IntegrateData(anchorset = integ_anchors,normalization.method = "SCT")
通常需要保存R对象。
#保存集成的seurat对象
saveRDS(seurat_integrated, "results/integrated_seurat.rds")
UMAP可视化
集成后,为了可视化集成数据,我们可以使用降维技术,例如PCA和UMAP。尽管PCA将确定所有PC,但我们一次只能绘制两个。相比之下,UMAP将从任何数量的顶级PC中获取信息,以将细胞布置在这个多维空间中。它将在多维空间中获取这些距离,并尝试将其绘制为二维。这样,细胞之间的距离代表表达上的相似性。
为了生成这些可视化,我们需要首先运行PCA和UMAP方法。
# 运行PCA
seurat_integrated <- RunPCA(object = seurat_integrated)
# 绘制PCA图
PCAPlot(seurat_integrated,split.by = "sample")