学习目标：
-构建质量控制指标，可视化地评估数据的质量
-应用适当的过滤器去除低质量的细胞
单细胞RNA-seq:质量控制
单细胞分析流程
这个工作流程的每一步都有它自己的目标和挑战。对于我们原始计数数据的QC，包括：
目标：
-过滤数据，只包含高质量的真正细胞，这样当我们聚集细胞时，更容易识别不同的细胞类型群体
-识别任何失败的样本、试图挽救数据或从分析中删除，此外，试图理解样本失败的原因
挑战：
-描述不太复杂的细胞中质量较差的细胞
-选择适当的筛选阈值，以保持高质量的细胞，而不去除生物学上相关的细胞类型
推荐:
-在执行QC之前，要充分了解您对出现的细胞类型的期望。例如，你希望在你的样本中有低复杂度的细胞还是线粒体表达水平较高的细胞吗?如果是这样，那么在评估我们的数据质量时，我们需要考虑解释这种生物学意义。
生成质量指标
当数据加载到Seurat并创建初始对象时，计数矩阵中的每个细胞都组装了一些基本元数据。为了仔细查看这个元数据，让我们查看存储在元数据中的数据帧。merged_seurat对象的数据槽:
# Explore merged metadata
View(merged_seurat@meta.data)
image.png
-orig.ident: 如果已知，此列将包含示例标识。它将默认为我们在加载数据时为project参数提供的值。
-nCount_RNA: 这一列表示每个细胞的UMIs数量
-nFeature_RNA: 这一栏表示每个细胞中检测到的基因数量
为了为质量控制分析创建适当的图，我们需要计算一些额外的指标。这些包括:
-number of genes detected per UMI: 这个度量让我们了解了我们数据集的复杂性(每个UMI检测到的基因越多，我们的数据就越复杂)。
-mitochondrial ratio: 这个指标将给我们一个源自线粒体基因的细胞reads百分比。
每UMI检测到的基因数
This value is quite easy to calculate, as we take the log10 of the number of genes detected per cell and the log10 of the number of UMIs per cell, then divide the log10 number of genes by the log10 number of UMIs.
基因数去log10除以每个细胞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()函数接受一个pattern参数，并在数据集中的所有基因标识符中搜索该pattern。因为我们寻找线粒体基因，寻找的是任何以" MT- "模式开始的基因标识符。对于每个细胞，该函数取属于“Mt-”集的所有基因(特征)的总数，然后除以所有基因(特征)的总数。将该值乘以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-”)适用于人类基因名称。您可能需要根据您感兴趣的组织来调整pattern参数。此外，如果你没有使用基因名称作为基因ID，那么这个函数就不能像我们在上面使用的那样工作，pattern是能起作用的。由于使用此函数有一些注意事项，建议手动计算此指标。如果你感兴趣，我们有代码可用来计算你自己的度量。
额外的元数据列
我们现在已经具备了评估数据所需的质量指标。但是，我们希望在元数据（metadata）中包含一些有用的附加信息，包括细胞id和环境信息。
当我们在上面的元数据文件中添加信息列时，我们只需使用$操作符将其直接添加到Seurat对象中的元数据槽中。我们可以在接下来的几列数据中继续这样做，但是我们将把这个数据帧提取到一个单独的变量中。通过这种方式，我们可以将元数据帧作为seurat对象的独立实体来处理，而不会影响对象中存储的任何其他数据。
让我们从提取元数据开始创建元数据数据框架。Seurat对象的数据槽:
# 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)
现在，您已经掌握了评估数据质量所需的指标!你最终的元数据表将有对应于每个细胞的行，以及包含这些细胞信息的列:
将更新的元数据保存到我们的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")
评估质量指标
现在我们已经生成了要评估的各种指标，我们可以用可视化的方式来研究它们。我们将评估各种指标，然后决定哪些细胞是低质量的，应该从分析中删除:
-Cell counts
-UMI counts per cell
-Genes detected per cell
-UMIs vs. genes detected
-Mitochondrial counts ratio
-Novelty
-双联体 在单细胞RNA测序实验中，由两个细胞产生双联体。它们通常是由于细胞分类或捕获错误而产生的，特别是在涉及数千个细胞的基于液滴的协议中。当目的是在单细胞水平上刻画种群特征时，双体显然是不可取的。特别是，它们可能错误地提示中间群体或实际上并不存在的过渡状态的存在。因此，双连体文库，以便他们不会危害的结果的解释。
-我们为什么不检查一下有没有双联体? 许多工作流使用UMIs或基因的最大阈值，其想法是，更多的读取或检测到的基因表明有多个细胞。虽然这一理论似乎很直观，但并不准确。此外，许多用于检测双重基因的工具倾向于去除具有中间或连续表型的细胞，尽管它们可能在非常离散的细胞类型的数据集上工作得很好。Scrublet是一种流行的双重检测工具，但是我们还没有对它进行充分的基准测试。目前，我们建议此时不包含任何阈值。当我们确定了每个簇的标记后，我们建议探索这些标记，以确定这些标记是否适用于不止一种细胞类型。
细胞计数
细胞计数是由检测到的唯一细胞barcodes的数量决定的。对于这个实验，预计需要12000 - 13000个细胞。
在理想的情况下，您希望唯一的细胞barcodes码的数量与您装载的细胞数量相对应。然而，情况并非如此，因为细胞的捕获率仅占所加载内容的一部分。例如，inDrops细胞的捕获效率更高(70-80%)，而10X在50-60%之间。
*注意:如果用于文库制备的细胞浓度不准确，则捕获效率可能会低得多。细胞浓度不应由FACS机器或生物分析仪测定(这些工具对浓度测定不准确)，而应使用血细胞计或自动细胞计数器来计算细胞浓度。
细胞数量也可以根据方案而变化，从而产生比我们加载的细胞数量高得多的细胞数量。例如，在inDrops方案中，细胞barcodes存在于水凝胶中，水凝胶被封装在具有单个细胞和裂解/反应混合物的液滴中。虽然每个水凝胶都应该有一个单独的细胞barcodes，但有时一个水凝胶可以有多个细胞barcodes。类似地，使用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.png
我们在每个样本中发现了超过15000个细胞，远远超过了预期的12- 13000个。很明显，我们可能有一些垃圾“细胞”存在。
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.png
我们可以看到，两个样本中的大多数细胞都有1000个或更多的UMIs，这很好。
每个细胞检测到的基因
我们对基因检测的期望与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.png
image.png
UMIs vs.检测到的基因
两个经常一起评估的指标是UMIs的数量和每个细胞检测到的基因数量。在这里，我们已经绘制了基因的数量与线粒体读取部分着色的UMIs的数量。线粒体读取分数只有在检测到很少基因的特别低的细胞中才高(深色的数据点)。这可能表明，受损或死亡的细胞，其胞质mRNA已通过破膜泄漏，因此，只有位于线粒体中的mRNA仍然是保守的。这些细胞被我们的计数和基因数量阈值过滤掉了。联合可视化计数和基因阈值显示联合过滤效应。
质量差的细胞很可能有低基因和每个细胞的UMIs，并与图中左下方象限的数据点相对应。良好的细胞通常表现为每个细胞有较多的基因和较多的UMIs。
通过这个图，我们还可以计算直线的斜率，以及图右下方象限中任意散点的数据。这些细胞有大量的UMIs，但只有少量的基因。这些细胞可能是正在死亡的细胞，但也可能代表一种低复杂度的细胞类型(即红细胞)。
# 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.png
线粒体数量比例
这个指标可以确定是否有大量的线粒体污染死亡或死亡细胞。我们将线粒体计数的低质量样本定义为超过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.png
复杂度
我们可以根据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.png
注意:每个细胞读取量是另一个有用的度量指标;然而，使用的工作流将需要保存该信息以进行评估。通常，使用这个度量标准，您希望看到所有的样本在每个细胞10,000到100,000个读数之间的相对相同位置上都有峰值。
过滤
总之，孤立地考虑这些QC度量中的任何一个都可能导致对细胞信号的误解。例如，线粒体计数相对较高的细胞可能参与呼吸过程，也可能是你想要保留的细胞。同样，其他指标也可以有其他的生物学解释。因此，在设置阈值时，始终考虑这些指标的联合影响，并将它们设置为尽可能宽松的，以避免无意中过滤掉活细胞群。
细胞层级过滤
现在我们已经可视化了各种度量，我们可以决定要应用的阈值，这将导致移除低质量的细胞。前面提到的建议通常只是一个粗略的指导方针，具体的实验需要告知所选择的确切阈值。我们将使用以下阈值:
-nUMI > 500
-nGene > 250
-log10GenesPerUMI > 0.8
-mitoRatio < 0.2
# 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))
基因层级过滤
在我们的数据中，我们将有许多零计数的基因。这些基因可以显著降低细胞的平均表达，所以我们将把它们从我们的数据中删除。我们将从每个细胞中哪些基因的计数为0开始:
# 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
现在，我们将根据流行程度执行一些过滤。如果一个基因只在少数细胞中表达，它就没有什么特别的意义，因为它仍然会降低所有其他没有表达的细胞的平均水平。对于我们的数据，我们选择只保留在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)
评估质量控制指标
在执行过滤之后，建议回顾这些指标，以确保您的数据符合您的期望，并且适合下游分析。(再做一遍）
练习
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")
网友评论