学习目标：

• 构建 QC 指标和相关图便于直观地查看数据质量
• 评估 QC 指标并设置过滤条件去掉低质量细胞

单细胞 RNA-seq：质量控制

image

此工作流程的每一步都有自己的目的和挑战。对于我们原始数据的 QC，它们包括：

目标：

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

挑战：

• 从不太复杂的细胞中划定质量差的细胞
• 选择合适的过滤阈值，保留高质量细胞，而不去除生物学相关细胞类型

建议：

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

生成质控标准

将数据加载到 Seurat 并创建初始对象时，会为矩阵中的每个细胞组装一些基本metadata。要仔细查看此metadata，让我们查看存储在seurat对象meta.data数据：

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

image

共有三列信息：

• orig.ident：通常包含样本标识（已知）。默认是project在加载数据时提供的值
• nCount_RNA：此列表示每个细胞的 UMI 数量
• nFeature_RNA：这一列代表每个细胞检测到的基因数量

为了为质量控制分析创建适当的图，我们需要计算一些额外的指标。这些包括：

• 每个 UMI 检测到的基因数量：这个指标让我们了解数据集的复杂性（每个 UMI 检测到的基因越多，我们的数据就越复杂）
• 线粒体比率：该指标将为我们提供细胞中线粒体基因的read的百分比

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

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

线粒体比率

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

对于我们的分析，我们更愿意使用比率值而不是使用百分比值。因此，我们将通过取输出值并除以 100 来反转函数执行的最后一步。

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

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

其它metadata列

我们现在拥有评估数据所需的质量指标。但是，我们希望在元数据中包含一些有用的附加信息，包括cell ID 和条件信息。

当我们在上面的元数据文件中添加信息列时，我们只是使用$运算符将其直接添加到 Seurat 对象中的元数据中。对于接下来的几列数据，我们可以继续这样做，但我们会将数据帧提取到一个单独的变量中。通过这种方式，我们可以将元数据数据帧作为与 seurat 对象分开的实体来使用，而不会影响存储在对象内的任何其他数据。

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

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

接下来，我们将为细胞标识符添加一个新列。来源于元数据数据框的行名称中。我们将按原样保留行名并将其复制到名为cells的新列中：

# Add cell IDs to metadata
metadata$cells <- rownames(metadata)

可以看到每个细胞 ID 都有一个ctrl_或stim_前缀，正如我们在合并 Seurat 对象时所指定的那样。我们可以使用此前缀创建一个新列，表示每个细胞在哪种条件下分类。我们将此列称为sample：

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

最后，我们将重命名元数据数据框中的一些现有列，以使其更直观：

# Rename columns
metadata <- metadata %>%
        dplyr::rename(seq_folder = orig.ident,
                      nUMI = nCount_RNA,
                      nGene = nFeature_RNA)

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

image

将更新的metadata保存到 Seurat 对象中

在我们评估我们的指标之前，我们要将工作保存到 Seurat 对象。我们可以通过简单地将数据帧分配到meta.data中：

# Add metadata back to Seurat object
merged_seurat@meta.data <- metadata

# Create .RData object to load at any time
save(merged_seurat, file="data/merged_filtered_seurat.RData")

评估质量指标

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

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

双联体细胞呢？在单细胞 RNA 测序实验中，双联体由两个细胞产生。它们通常是由于细胞分选或捕获中的错误引起的，尤其是在涉及数千个细胞的基于液滴的protocol中。当目标是在单细胞水平上表征种群时，双峰显然是不可取的。特别是，它们可能会错误地暗示存在实际上并不存在的中间种群或过渡状态。因此，最好去除双峰库，这样它们就不会影响对结果的解释。

为什么我们不检查双峰？许多工作流程使用 UMI 或基因的最大阈值，认为检测到的读数或基因数量多得多，表明存在多个细胞。虽然这个理由看起来很直观，但并不准确。此外，许多用于检测双峰的工具往往会去除具有中间或连续表型的细胞，尽管它们可能适用于具有非常离散的细胞类型的数据集。Scrublet是一种流行的双峰检测工具，但我们还没有对它进行充分的基准测试。目前，我们建议此时不包括任何阈值。当我们确定了每个簇的标记后，我们建议探索这些标记以确定这些标记是否适用于一种以上的细胞类型。

细胞计数

细胞计数由检测到的唯一细胞barcode的数量决定。对于该实验，预计有 12,000 -13,000 个细胞。

在理想情况下，您希望barcode数量与细胞数量一一对应。然而，情况并非如此，因为细胞的捕获率只是加载内容的一部分。例如，与10X 相比，inDrops 细胞捕获效率更高（70-80%），而10X的效率在50-60% 之间 。

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

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

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

image

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

每个细胞的 UMI 计数（转录本）

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

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

image

我们可以看到两个样本中的大多数细胞都具有 1000 个或更多的 UMI，这很好。

每个细胞检测到的基因

我们对基因检测的期望与 UMI 检测相似，尽管它可能比 UMI 低一点。对于高质量数据，比例直方图应包含代表被封装细胞的单个大峰。如果我们看到主峰左侧的小肩峰（不存在于我们的数据中），或细胞的双峰分布，则可以表明一些情况。可能有一组细胞由于某种原因失败了。也可能存在生物学上不同类型的细胞（即静止细胞群、目标复杂度较低的细胞），或一种类型比另一种小得多（即具有高计数的细胞可能是尺寸较大的细胞）。因此，应使用我们在本课中描述的其他指标来评估此阈值。

# Visualize the distribution of genes detected per cell via histogram
metadata %>% 
    ggplot(aes(color=sample, x=nGene, fill= sample)) + 
    geom_density(alpha = 0.2) + 
    theme_classic() +
    scale_x_log10() + 
    geom_vline(xintercept = 300)

# Visualize the distribution of genes detected per cell via boxplot
metadata %>% 
    ggplot(aes(x=sample, y=log10(nGene), fill=sample)) + 
    geom_boxplot() + 
    theme_classic() +
    theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
    theme(plot.title = element_text(hjust=0.5, face="bold")) +
    ggtitle("NCells vs NGenes")

image

image

UMIs 与基因检测

通常一起评估的两个指标是 UMI 的数量和每个细胞检测到的基因数量。在这里，我们绘制了基因数量与由线粒体read分数着色的 UMI 数量的关系图。线粒体read分数仅在检测到的基因很少（深色数据点）的特别低计数细胞中很高。这可能表明受损/垂死的细胞的细胞质 mRNA 通过破损的膜泄漏，因此，只有位于线粒体中的 mRNA 仍然是保守的。这些细胞被我们的计数和基因数量阈值过滤掉。联合可视化计数和基因阈值显示联合过滤效果。

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

通过这个图，我们还评估了线的斜率，以及图右下象限中数据点的任何散点。这些细胞具有大量的 UMI，但只有少数基因。这些可能是垂死的细胞，但也可能代表低复杂性细胞类型（即红细胞）的群体。

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

image

线粒体计数比

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

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

image

复杂性

我们可以通过使用称为新颖性评分的衡量标准，根据 RNA 种类的复杂程度来评估每个细胞。新颖性分数是通过 nGenes 与 nUMI 的比率来计算的。如果有很多捕获的转录本（高 nUMI）和在细胞中检测到的基因数量很少，这可能意味着您只捕获了少量基因，并且只是一遍又一遍地对这些较少数量的基因的转录本进行测序。这些低复杂性（低新颖性）细胞可能代表特定的细胞类型（即缺乏典型转录组的红细胞），或者可能是由于一些其他奇怪的人工制品或污染。通常，我们预计优质细胞的新颖性得分高于 0.80。

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

image

注意： 每个细胞的reads是另一个可用于探索的指标；但是，使用的工作流程需要保存这些信息以进行评估。通常，使用此指标，您希望看到在每个细胞 10,000 到 100,000 个读数之间的相对相同位置具有峰值的所有样本。

过滤

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

细胞水平过滤

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

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

在 Seurat 对象中使用该subset()函数过滤：

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

基因水平过滤

在我们的数据中，我们将有许多零计数的基因。这些基因可以明显的降低细胞的平均表达，因此我们将从我们的数据中删除它们。我们将首先确定每个细胞中哪些基因的计数为零：

# Extract counts
counts <- GetAssayData(object = filtered_seurat, slot = "counts")

# Output a logical matrix specifying for each gene on whether or not there are more than zero counts per cell
nonzero <- counts > 0

现在，我们将按prevalence进行过滤。如果一个基因只在少数细胞中表达，它并不是特别有意义，因为它仍然会降低所有其他细胞的平均值。对于我们的数据，我们选择只保留在 10 个或更多表达的基因细胞。通过使用这个过滤器，所有细胞中计数为零的基因将被有效地去除。

# Sums all TRUE values and returns TRUE if more than 10 TRUE values per gene
keep_genes <- Matrix::rowSums(nonzero) >= 10

# Only keeping those genes expressed in more than 10 cells
filtered_counts <- counts[keep_genes, ]

最后，获取这些过滤后的计数并创建一个新的 Seurat 对象用于下游分析。

# Reassign to filtered Seurat object
filtered_seurat <- CreateSeuratObject(filtered_counts, meta.data = filtered_seurat@meta.data)

重新评估 QC 指标

执行过滤后，建议查看相关指标以确保您的数据符合您的期望并且有利于下游分析。

练习

1. 使用下面提供的代码从过滤后的 Seurat 对象中提取新的元数据：

```
 # Save filtered subset to new metadata
 metadata_clean <- filtered_seurat@meta.data

```

2. 使用过滤后的数据执行所有相同的 QC 绘图。

3. 报告每个样本剩余的细胞数，并评论移除的细胞数是高还是低。你能解释为什么这个数字仍然不是~12K（这是为实验加载了多少个细胞）？

4. 每个细胞的 nGene 过滤后，您仍应观察到主峰右侧的小肩峰。这个肩膀代表什么？

5. 当根据 nUMI 绘制 nGene 时，您是否观察到图右下象限中的任何数据点？你对这些被移除的细胞有什么想法？

保存过滤的细胞

根据这些 QC 指标，我们将识别任何不合格的样本，并继续处理我们过滤的细胞。通常我们会使用不同的过滤标准来遍历 QC 指标；它不一定是一个线性过程。当满足过滤条件时，我们将保存过滤后的细胞对象以进行聚类和标记识别。

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

注意：我们使用的数据质量非常好。如果您有兴趣了解在执行 QC 时“不良”数据可能是什么样子，我们在此处链接了一些材料，我们在其中探索了劣质样品的类似 QC 指标。