scRNA-seq—质量控制

分享是一种态度

回顾

单细胞RNA-seq分析介绍 单细胞RNA-seq的设计和方法 从原始数据到计数矩阵 差异分析前的准备工作 scRNA-seq—读入数据详解

学习目标

• 构建QC指标和相关的图，以直观地了解数据的质量
• 评估QC指标并设置过滤条件以删除低质量的细胞

scRNA-seq质量控制流程

此工作流程的每个步骤都有自己的目的和挑战。对于原始计数数据的质量控制，包括：

目标

• 筛选数据，使其仅包含高质量的真实细胞，这样当我们对细胞进行聚类时，就更容易识别不同的细胞类群
• 识别任何不合格的样本，并尝试挽救数据或将其从分析中删除，此外，还要尝试了解样本失败的原因

挑战

• 从不太复杂的细胞中划定质量较差的细胞
• 选择合适的阈值进行过滤，以保持高质量的细胞，而不会去除与生物相关的细胞类型

建议

• 在执行QC之前，要求您对细胞类型的期望有一个很好的了解。例如，您是否希望在您的样本中有低复杂性的细胞或具有较高水平的线粒体表达的细胞？如果是这样，那么在评估我们的数据质量时，我们需要考虑到这种生物学因素。

生成质控标准

请记住，Seurat会自动为每个细胞创建一些元数据：

1# Explore merged metadata
2View(merged_seurat@meta.data)

添加的列

• orig.ident ：通常包含样本标识（如果已知），通常默认project为我们为其分配的身份
• nCount_RNA ：每个细胞的UMI数量
• nFeature_RNA ：每个细胞检测到的基因数量

我们需要计算一些用于绘图的其他指标：

• number of genes detected per UMI: 这个度量让我们对数据集的复杂性有了一个概念(每个UMI检测到的基因越多，我们的数据就越复杂)
• mitochondrial ratio: 这个度量将给我们提供一个源自线粒体基因的细胞读数的百分比。

每个细胞的每个UMI的基因数量非常容易计算，我们将对结果进行log10转换，以便更好地在样本之间进行比较。

1# Add number of genes per UMI for each cell to metadata
2merged_seurat$log10GenesPerUMI <- log10(merged_seurat$nFeature_RNA) / log10(merged_seurat$nCount_RNA)

Seurat有一个方便的功能，可以让我们计算映射到线粒体基因的转录本的比例。PercentageFeatureSet()将采用一定模型并搜索基因标识符。此功能可以轻松计算属于每个细胞的可能功能子集的所有计数的百分比。这里的计算只是将属于该集合的要素的计数槽中存在的矩阵的列和除以所有要素的列和，然后乘以100。由于我们要绘制比率值，所以我们将反转这一步，然后除以100

1# Compute percent mito ratio
2merged_seurat$mitoRatio <- PercentageFeatureSet(object = merged_seurat, pattern = "^MT-")
3merged_seurat$mitoRatio <- merged_seurat@meta.data$mitoRatio / 100

注意：提供的模式（“ ^ MT-”）适用于人类基因名称。您可能需要根据您感兴趣的生物进行调整。如果您没有使用基因名称作为基因ID，则此功能将无法使用。我们有代码可用于自行计算该指标。

同时还需要将其他信息添加到QC指标的元数据中，例如单元ID、条件信息和各种指标。虽然使用$操作符将信息直接添加到Seurat对象的元数据槽非常容易，但是我们选择把数据框提取到一个单独的变量中。通过这种方式，我们可以继续插入QC分析所需的其他指标，而不会有影响merded_seurat对象的风险。

我们将通过从seurat对象提取meta.data槽来创建元数据数据框：

1# Create metadata dataframe
2metadata <- merged_seurat@meta.data

可以看到每个细胞ID都有一个ctrl_或stim_前缀，正如我们在合并Seurat对象时指定的那样。这些前缀应该与Orig.ident中列出的样本匹配。让我们首先添加一个包含细胞ID的列，并更改当前列名以使其更直观：

1# Add cell IDs to metadata
2metadata$cells <- rownames(metadata)
3
4# Rename columns
5metadata <- metadata %>%
6        dplyr::rename(seq_folder = orig.ident,
7                      nUMI = nCount_RNA,
8                      nGene = nFeature_RNA)

现在，根据细胞前缀获取每个细胞的样本名称

1# Create sample column
2metadata$sample <- NA
3metadata$sample[which(str_detect(metadata$cells, "^ctrl_"))] <- "ctrl"
4metadata$sample[which(str_detect(metadata$cells, "^stim_"))] <- "stim"

现在，您已经设置好了评估数据质量所需的指标！最终的元数据表将包含对应于每个细胞的行，以及包含有关这些细胞的信息的列：

将更新的元数据保存到我们的Seurat对象

在评估指标之前，我们可以把迄今为止完成的所有工作保存回Seurat对象中，这样方便以后调用。我们只需将数据框分配到meta.data插槽即可完成此操作。

1# Add metadata back to Seurat object
2merged_seurat@meta.data <- metadata
3
4# Create .RData object to load at any time
5save(merged_seurat, file="data/merged_filtered_seurat.RData")

评估质量指标

现在我们已经生成了要评估的各种指标，我们可以通过可视化来探索它们。我们将评估各种指标，然后决定哪些细胞质量较低，应该从分析中删除：

• 细胞计数
• 每个细胞的UMI计数
• 每个细胞检测到的基因
• UMI与检测到的基因
• 线粒体比率
• Novelty

What about doublets?在单细胞RNA测序实验中，双胞体是由两个细胞产生的。它们通常是由于细胞分选或捕获中的错误引起的，特别是在涉及数千个细胞的基于液滴的协议中。当目标是描述单细胞水平的群体特征时，双峰显然是不可取的。具体地说，他们可能错误地暗示存在实际并不存在的中间群体或短暂状态。因此，需要移除双峰文库，以便它们不会影响结果的解释。 Why aren’t we checking for doublets? 许多工作流程使用UMI或基因的最大阈值，其想法是检测到的大量读数或基因表明存在多个细胞。尽管此理由似乎很直观，但并不准确。同样，许多用于检测双峰的工具倾向于去除具有中间或连续表型的细胞，尽管它们在具有非常离散的细胞类型的数据集上可能会很好地工作。Scrublet是双态检测的流行工具，但我们尚未对其进行充分的基准测试。目前，我们建议您此时不设置任何阈值。当我们确定了每个簇的标记时，建议您探索这些标记，以确定这些标记是否适用于一种以上的细胞类型。

Cell counts

细胞计数由检测到的唯一细胞条形码的数量确定。对于此实验，预计将有12,000 -13,000个细胞。

在理想的情况下，您希望唯一的细胞条形码的数量与您加载的细胞的数量相对应。但是，情况并非如此，因为细胞的捕获率仅是所装载细胞的比例的一部分。例如，与10x相比inDrops的细胞捕获效率更高(70-80%)，而10x的效率可以下降到50-60%。

注意：如果用于文库制备的细胞浓度不准确，捕获效率可能会显得很低。细胞浓度不应由FACS机器或生物分析仪测定(这些工具对浓度测定不准确)，而应使用血细胞计数仪或自动细胞计数器计算细胞浓度。

细胞编号也可能因方法而异，产生的细胞编号比我们加载的要高得多。例如，在inDrops方法中，细胞条形码存在于水凝胶中，并与单个细胞和裂解/反应混合物封装在液滴中。虽然每个水凝胶都应该有一个与之相关的细胞条形码，但有时一个水凝胶可以有多个细胞条形码。类似地，使用10X协议，有可能只获得乳液液滴(GEM)中的条形码珠子，而没有实际的细胞。这两种情况，加上死亡细胞的存在，都可能导致比细胞更多的细胞条形码。

1# Visualize the number of cell counts per sample
2metadata %>% 
3      ggplot(aes(x=sample, fill=sample)) + 
4      geom_bar() +
5      theme_classic() +
6      theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
7      theme(plot.title = element_text(hjust=0.5, face="bold")) +
8      ggtitle("NCells")

我们看到每个样本都超过15，000个细胞，这比预期的12-13,000个要多得多。很明显，我们现在可能有一些垃圾“细胞”。

UMI counts (transcripts) per cell

每个细胞的UMI计数通常应高于500，这是我们预期的下限。如果UMI计数在500-1000计数之间，则可以使用，但可能应该对细胞进行更深的测序。

1# Visualize the number UMIs/transcripts per cell
2metadata %>% 
3      ggplot(aes(color=sample, x=nUMI, fill= sample)) + 
4      geom_density(alpha = 0.2) + 
5      scale_x_log10() + 
6      theme_classic() +
7      ylab("Cell density") +
8      geom_vline(xintercept = 500)

我们可以看到，两个样本中的大多数单元都具有1000个UMI或更高，这非常好。

Genes detected per cell

我们对基因检测的期望值与UMI检测的期望值相似，尽管可能比UMIs略低。对于高质量数据，比例直方图应包含表示封装的细胞的单个大峰值。如果我们看到主峰右边的一个小肩膀(我们的数据中没有出现)，或者细胞的双峰分布，这可能表明了一些事情。可能是因为某些原因导致一组细胞发生错误。也可能是存在生物上不同类型的细胞(即，静止的细胞群体、不太复杂的目标细胞)，还有就是可能是一种类型比另一种小得多(即，高计数的细胞可能是尺寸较大的细胞)。因此，应该使用我们在本课中描述的其他指标来评估此阈值。

UMIs vs. genes detected

通常一起评估的两个指标是UMI的数量和每个细胞检测到的基因数量。在这里，我们绘制了基因数量与线粒体reads所占比例UMI数量的关系图。线粒体reads部分只在检测到很少基因的特别低计数的细胞中才高(浅蓝色)。这可能是损伤/死亡的细胞，其细胞质的mRNA已经通过破裂的膜泄漏出来，因此，只有位于线粒体的mRNA仍然是保守的。这些细胞被我们的计数和基因数量阈值过滤掉。联合可视化计数和基因阈值可显示联合过滤效果。

质量差的细胞很可能每个细胞的基因和UMI都很低，并且与图左下象限的数据点相对应。好的细胞通常会表现为每个细胞有更多的基因和更高数量的UMI。

通过此图，我们还评估了线的斜率，以及图的右下角象限中数据点的任何散布情况。这些细胞有大量的UMI，但只有几个基因。这些可能是濒临死亡的细胞，但也可能代表一个低复杂性细胞类型的群体(即红细胞)。

 1# Visualize the correlation between genes detected and number of UMIs and determine whether strong presence of cells with low numbers of genes/UMIs
 2metadata %>% 
 3      ggplot(aes(x=nUMI, y=nGene, color=mitoRatio)) + 
 4      geom_point() + 
 5    scale_colour_gradient(low = "gray90", high = "black") +
 6      stat_smooth(method=lm) +
 7      scale_x_log10() + 
 8      scale_y_log10() + 
 9      theme_classic() +
10      geom_vline(xintercept = 500) +
11      geom_hline(yintercept = 250) +
12      facet_wrap(~sample)

Mitochondrial counts ratio

这一指标可以确定死亡或濒临死亡的细胞是否存在大量的线粒体污染。我们将线粒体计数质量差的样品定义为超过0.2线粒体比率标记的细胞，除非您希望样品中有这种情况。

1# Visualize the distribution of mitochondrial gene expression detected per cell
2metadata %>% 
3      ggplot(aes(color=sample, x=mitoRatio, fill=sample)) + 
4      geom_density(alpha = 0.2) + 
5      scale_x_log10() + 
6      theme_classic() +
7      geom_vline(xintercept = 0.2)

Complexity

我们可以看到，我们对每个细胞测序较少的样本具有更高的整体复杂性，这是因为我们还没有开始对这些样本的任何给定基因进行饱和测序。这些样本中的异常值细胞可能是RNA种类比其他细胞简单的细胞。有时，我们可以通过这个指标检测到低复杂性细胞类型的污染，比如红细胞。一般来说，我们预计novelty得分在0.80以上。

1# Visualize the overall complexity of the gene expression by visualizing the genes detected per UMI
2metadata %>%
3      ggplot(aes(x=log10GenesPerUMI, color = sample, fill=sample)) +
4      geom_density(alpha = 0.2) +
5      theme_classic() +
6      geom_vline(xintercept = 0.8)

注意：Reads per cell是另一个值得研究的指标；但是，使用的工作流需要保存此信息以供评估。通常，可以使用此度量标准来查看所有样本，每个样本的峰值在相对相同的位置，每个细胞的读数介于10,000和100,000之间。

过滤

总之，孤立地考虑这些质量控制指标中的任何一个都可能导致对细胞信号的误解。例如，线粒体计数相对较高的细胞可能参与呼吸过程，可能是您想要保留的细胞。同样，其他指标也可以有其他生物学解释。因此，在设置阈值时，请始终考虑这些指标的共同影响，并将其设置为尽可能宽松，以避免无意中过滤掉可行的细胞群体。

Cell-level filtering

现在我们已经可视化了各种指标，我们可以决定要应用的阈值，这将导致删除低质量的细胞。通常，前面提到的建议是一个粗略的指导原则，而具体的实验需要设定所选的确切阈值。我们将使用以下阈值：

• nUMI > 500
• nGene > 250
• log10GenesPerUMI > 0.8
• mitoRatio < 0.2

进行过滤，使用以下subset()函数：

1# Filter out low quality reads using selected thresholds - these will change with experiment
2filtered_seurat <- subset(x = merged_seurat, 
3                         subset= (nUMI >= 500) & 
4                           (nGene >= 250) & 
5                           (log10GenesPerUMI > 0.80) & 
6                           (mitoRatio < 0.20))

Gene-level filtering

在我们的数据中，我们将有许多零计数的基因。这些基因可以极大地降低细胞的平均表达量，所以我们将把它们从我们的数据中删除。首先，我们将删除所有细胞中零表达的基因。此外，我们还将根据prevalence执行一些过滤。如果一个基因仅在少数细胞中表达，那么它就没有什么特别的意义，因为它仍会降低未在其中表达的所有其他细胞的平均值。对于我们的数据，我们选择只保留在10个或更多细胞中表达的基因。

 1# Output a logical vector for every gene on whether the more than zero counts per cell
 2# Extract counts
 3counts <- GetAssayData(object = filtered_seurat, slot = "counts")
 4
 5# Output a logical vector for every gene on whether the more than zero counts per cell
 6nonzero <- counts > 0
 7
 8# Sums all TRUE values and returns TRUE if more than 10 TRUE values per gene
 9keep_genes <- Matrix::rowSums(nonzero) >= 10
10
11# Only keeping those genes expressed in more than 10 cells
12filtered_counts <- counts[keep_genes, ]
13
14# Reassign to filtered Seurat object
15filtered_seurat <- CreateSeuratObject(filtered_counts, meta.data = filtered_seurat@meta.data)

重新评估QC指标

执行过滤之后，建议回顾一下指标，以确保您的数据符合您的期望并适合进行下游分析。

保存过滤后的细胞

基于这些QC指标，我们将识别出任何不合格的样本，并继续进行过滤后的单元格。通常，我们使用不同的过滤条件来迭代QC指标。它不一定是线性过程。满足过滤条件后，我们将保存过滤后的细胞对象以进行聚类和标记识别。

1# Create .RData object to load at any time
2save(filtered_seurat, file="data/seurat_filtered.RData")

文中链接

[1] 代码可用于自行计算该指标: https://github.com/hbctraining/scRNA-seq/blob/master/lessons/mitoRatio.md [2] Scrublet: https://github.com/AllonKleinLab/scrublet