本文来源于哈佛大学的单细胞课程系列,在此做一些学习,不当之处请指正。
scRNA-seq/04_SC_quality_control.md at master · hbctraining/scRNA-seq · GitHub
https://github.com/hbctraining/scRNA-seq/blob/master/lessons/04_SC_quality_control.md
构造质量控制指标和相关图以可视化方式探索数据质量
评估质量控制指标并过滤以删除低质量的细胞
此工作流程的每个步骤都有其自己的目标和挑战。对于我们的原始计数数据的质量控制,它们包括:
目标:
要过滤数据以仅包括高质量的真实细胞,以便在对细胞进行聚类时,更容易识别不同的细胞类型种群
要确定任何不合格的样品,要么试图挽救数据或分析除去,除了,试图理解为什么样品不合格
挑战:
从复杂程度较低的细胞中划定质量较差的细胞
选择适当的过滤阈值,以便在不去除生物学相关细胞类型的情况下保持高质量细胞
建议:
在执行质量控制之前,对您对存在的细胞类型的期望有个好主意。例如,您是否期望样品中具有低复杂度的细胞或具有较高线粒体表达水平的细胞?如果是这样,那么在评估我们的数据质量时,我们需要考虑这种生物学。
请记住,Seurat为每个细胞自动创建一些元数据:
#探索合并的元数据
View(merged_seurat@meta.data)
添加的列包括:
orig.ident
:通常包含样本名称(如果已知),但默认project
为我们为其分配的身份
nCount_RNA
:每个细胞的UMI数量
nFeature_RNA
:每个细胞检测到的基因数量
我们需要计算一些其他度量标准以进行绘图:
每个UMI检测到的基因数量:此度量标准使我们对数据集的复杂性有所了解(每个UMI检测到的基因更多,我们的数据更复杂)
线粒体比率:该指标将为我们提供一定百分比的源自线粒体基因的细胞读数
每个细胞每个UMI的基因数量非常容易计算,我们将对结果进行log10转换以更好地进行样本之间的比较。
#将每个细胞的每个UMI的基因数量添加到元数据
merged_seurat$log10GenesPerUMI <- log10(merged_seurat$nFeature_RNA) / log10(merged_seurat$nCount_RNA)
Seurat具有便利的功能,可以让我们计算映射到线粒体基因的转录本的比例。在PercentageFeatureSet()
将采取的模式和搜索基因标识。对于每个列(细胞),它将使用属于该集合的要素的计数槽之和,除以所有要素的列之和,然后乘以100。由于我们需要比率值进行绘制,因此我们将这一步反向进行然后除以100。
#计算线粒体比率
merged_seurat$mitoRatio <- PercentageFeatureSet(object = merged_seurat, pattern = "^MT-")
merged_seurat$mitoRatio <- merged_seurat@meta.data$mitoRatio / 100
注意:提供的模式(“ ^ MT-”)适用于人类基因名称。您可能需要根据您感兴趣的生物体进行调整。如果您没有使用基因名称作为基因ID,则此功能将无法使用。我们有代码可用于自行计算该指标。
我们需要为我们的质量控制指标向元数据中添加其他信息,例如细胞ID,条件信息和各种指标。尽管使用$
操作符将信息直接添加到Seurat对象的元数据槽中非常容易,但我们将数据框提取到一个单独的变量中。这样,我们可以继续插入进行质量控制分析所需的其他指标,而不会影响merged_seurat
对象。
我们将通过从Seurat对象中提取meta.data
来创建元数据数据框:
#创建元数据
metadata <- merged_seurat@meta.data
您应该看到每个细胞ID在合并Seurat对象时都指定了一个ctrl_
或stim_
前缀。这些前缀应与中列出的示例匹配orig.ident
。让我们开始添加带有细胞ID的列,并更改当前列的名称以使其更加直观:
# 添加cell IDs 到 metadata
metadata$cells <- rownames(metadata)
# 重命名列名
metadata <- metadata %>%
dplyr::rename(seq_folder = orig.ident,
nUMI = nCount_RNA,
nGene = nFeature_RNA)
现在,让我们根据细胞前缀获取每个细胞的样本名称:
#创建样品列
metadata$sample <- NA
metadata$sample[which(str_detect(metadata$cells, "^ctrl_"))] <- "ctrl"
metadata$sample[which(str_detect(metadata$cells, "^stim_"))] <- "stim"
现在,您已经完成了评估数据质量所需的度量标准的设置!最终的元数据表将具有与每个细胞相对应的行,以及包含有关这些细胞信息的列:
在评估指标之前,我们将把迄今为止完成的所有工作保存到Seurat对象中。我们可以通过简单地将数据帧分配到meta.data
插槽中来做到这一点:
# 将metadata添加回Seurat对象
merged_seurat@meta.data <- metadata
# 保存.RData object
save(merged_seurat, file="data/merged_filtered_seurat.RData")
现在,我们已经生成了各种要评估的指标,我们可以通过可视化对其进行探索。我们将评估各种指标,然后决定哪些电池质量较低,应从分析中将其删除:
细胞计数
每个细胞的UMI计数
每个细胞检测到的基因
UMI与检测到的基因
线粒体计数比
新奇
那双细胞呢?在单细胞RNA测序实验中,两个细胞产生了双峰。它们通常是由于细胞分选或捕获中的错误而产生的,尤其是在涉及数千个细胞的基于液滴的协议中。当目标是在单细胞水平上表征种群时,双峰显然是不可取的。特别是,它们可能错误地建议存在实际上不存在的中间人口或过渡状态。因此,希望删除双峰文库,以使它们不影响结果的解释。
我们为什么不检查双峰呢?许多工作流程使用UMI或基因的最大阈值,其想法是检测到的大量读数或基因表明存在多个细胞。尽管此理由似乎很直观,但并不准确。同样,许多用于检测双峰的工具倾向于去除具有中间或连续表型的细胞,尽管它们在具有非常离散的细胞类型的数据集上可以很好地工作。Scrublet是用于双态检测的流行工具,但我们尚未对其进行充分的基准测试。目前,我们建议您此时不包括任何阈值。当我们确定了每个簇的标记后,我们建议您探索这些标记,以确定这些标记是否适用于一种以上的细胞类型。
细胞计数由检测到的独特细胞条形码的数量确定。对于此实验,预计将有12,000 -13,000个细胞。
在理想的世界中,您希望唯一的细胞条形码的数量与您加载的细胞的数量相对应。但是,情况并非如此,因为细胞的捕获率仅是所装载细胞的比例的一部分。例如,与10X相比,inDrops细胞的捕获效率更高(70-80%),后者可能下降50-60%。
注意:如果用于文库制备的细胞浓度不准确,捕获效率可能会低得多。不应通过FACS机器或生物分析仪确定细胞浓度(这些工具无法准确确定浓度),而应使用血细胞计数器或自动细胞计数器来计算细胞浓度。
信元号也可以根据协议而变化,从而产生比我们加载的要高得多的信元号。例如,在inDrops协议中,细胞条形码存在于水凝胶中,并与单个细胞和裂解/反应混合物一起包裹在液滴中。尽管每个水凝胶应具有与其关联的单个细胞条形码,但有时水凝胶可以具有多个细胞条形码。类似地,使用10X协议,有可能仅在乳液液滴(GEM)中获得条形码条形码的珠,而没有实际的细胞。除存在即将死亡的细胞外,这两种方法均可导致比细胞数量更多的细胞条形码。
#可视化每个样本
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")
我们看到每个样本超过15,000个细胞,这大大超过了预期的12-13,000。显然,我们可能存在一些垃圾“细胞”。
每个细胞的UMI计数通常应高于500,这是我们预期的下限。如果UMI计数在500-1000计数之间,则可以使用,但可能应该对细胞进行更深的测序。
#可视化每个细胞
# 可视化每个细胞的UMIs/transcripts
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)
我们可以看到,两个样本中的大多数细胞都具有1000个UMI或更高,这非常好。
我们对基因检测的期望与对UMI检测的期望相似,尽管它可能比UMI更低。对于高质量数据,比例直方图应包含一个单个大峰,该峰代表封装的细胞。如果我们在主要峰的右侧看到一个小的肩膀(在我们的数据中不存在),或者细胞的双峰分布,则可能表明有两件事。可能是因为某些原因而导致一组细胞失败。也可能是有生物学上不同类型的细胞(例如,静态细胞群,不太复杂的目标细胞)和/或一种类型比另一种小得多(即,计数高的细胞可能是尺寸较大的细胞)。因此,应该使用本课中介绍的其他指标来评估此阈值。
# 通过密度图可视化每个细胞检测到的基因的分布
metadata %>%
ggplot(aes(color=sample, x=nGene, fill= sample)) +
geom_density(alpha = 0.2) +
theme_classic() +
scale_x_log10() +
geom_vline(xintercept = 300)
# 通过箱线图可视化每个细胞检测到的基因的分布
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")
通常一起评估的两个指标是UMI的数量和每个细胞检测到的基因数量。在这里,我们以线粒体读段的比例绘制了基因数量与UMI数量的关系图。线粒体阅读分数在具有很少检测到的基因的特别低计数的细胞中仅高(浅蓝色)。这可能表明细胞质的mRNA已经通过破损的膜漏出了受损/垂死的细胞,因此,仅位于线粒体中的mRNA仍然是保守的。这些细胞被我们的计数和基因数量阈值滤除。联合可视化计数和基因阈值可显示联合过滤效果。
细胞是质量差的可能有低的基因,而每个细胞的UMI,并对应于数据点在情节的左下象限。好的细胞通常会同时显示每个细胞更多的基因和更高数量的UMI。
通过此图,我们还可以评估线的斜率以及该图右下象限中数据点的任何散布。这些细胞具有大量的UMI,但只有少数基因。这些可能是垂死的细胞,但也可能代表低复杂度细胞类型(即红细胞)的种群。
#可视化检测到的基因与UMI数量之间的相关性,并确定是否存在大量基因/ UMI
元数据较少的细胞% >%
# 可视化检测到的基因与UMI数量之间的相关性,并确定是否存在大量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)
该度量可以识别死细胞或垂死细胞中是否存在大量的线粒体污染。我们将线粒体计数质量差的样品定义为超过0.2的线粒体比率标记的细胞,除非您当然希望样品中有这种情况。
# 可视化每个细胞元数据中检测到的线粒体基因表达的分布
metadata %>%
ggplot(aes(color=sample, x=mitoRatio, fill=sample)) +
geom_density(alpha = 0.2) +
scale_x_log10() +
theme_classic() +
geom_vline(xintercept = 0.2)
我们可以看到对每个细胞进行较少测序的样品具有较高的总体复杂性,这是因为我们尚未开始为这些样品的任何给定基因饱和测序。这些样品中的异常细胞可能是具有比其他细胞复杂的RNA种类的细胞。有时,我们可以通过此指标检测低复杂度的细胞类型(如红细胞)的污染。通常,我们希望新奇分数高于0.80。
#通过可视化每个UMI元数据检测到的基因来可视化基因表达的整体复杂性
# Filter out low quality reads using selected thresholds - these will change with experiment
metadata %>%
ggplot(aes(x=log10GenesPerUMI, color = sample, fill=sample)) +
geom_density(alpha = 0.2) +
theme_classic() +
geom_vline(xintercept = 0.8)
注意: 每个细胞的读取数是另一个有用的度量标准;但是,所使用的工作流程将需要保存此信息以进行评估。通常,您希望使用此度量标准来查看所有样本,每个样本的峰值在相对相同的位置,每个细胞的读数介于10,000和100,000之间。
总之,孤立地考虑任何这些QC指标都可能导致对细胞信号的误解。例如,线粒体计数相对较高的细胞可能参与呼吸过程,并且可能是您想要保留的细胞。同样,其他指标可以具有其他生物学解释。因此,在设置阈值时,请始终考虑这些指标的共同影响,并将其设置为尽可能宽松,以避免无意间滤除可行的细胞群。
现在我们已经可视化了各种指标,我们可以决定要应用的阈值,这将导致去除低质量的细胞。通常,前面提到的建议是一个粗略的指导原则,具体的实验需要告知所选择的确切阈值。我们将使用以下阈值:
nUMI > 500
nGene > 250
log10GenesPerUMI > 0.8
mitoRatio < 0.2
要进行过滤,我们将回到Seurat对象并使用以下subset()
函数:
#筛选出低质量的读出使用所选择的阈值-这些将与实验改变
# Filter out low quality reads 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))
在我们的数据中,我们将有许多计数为零的基因。这些基因可以显着降低细胞的平均表达,因此我们将其从数据中删除。首先,我们将删除所有细胞中零表达的基因。此外,我们将按流行程度执行一些过滤。如果一个基因仅在少数几个细胞中表达,那么它就没有什么特别的意义,因为它仍会降低未在其中表达的所有其他细胞的平均值。对于我们的数据,我们选择仅保留在10个或更多的基因中表达的基因细胞。
#为每个基因输出一个逻辑向量,确定每个细胞是否计数超过零
#提取计数
counts <- GetAssayData(object = filtered_seurat, slot = "counts")
#为每个基因输出一个逻辑向量,以了解每个细胞的计数是否大于零,
nonzero <- counts > 0
# 总结TRUE TRUE所有值和回报,如果每个基因超过10个TRUE值
keep_genes <- Matrix::rowSums(nonzero) >= 10
# 只有保持在超过10个细胞中表达的那些基因
filtered_counts <- counts[keep_genes, ]
# 重新分配给已过滤的Seurat对象
filtered_seurat <- CreateSeuratObject(filtered_counts, meta.data = filtered_seurat@meta.data)
执行过滤之后,建议回顾一下指标,以确保您的数据符合您的期望并适合进行下游分析。
练习题
从过滤的Seurat对象中提取新的元数据,以进行与未过滤的数据相同的绘图
# 保存过滤子集,以新的元数据
metadata_clean <- filtered_seurat@meta.data
执行与未过滤数据相同的所有图,并确定所使用的阈值是否合适。
基于这些QC指标,我们将识别出任何失败的样本,并继续进行过滤后的细胞。通常,我们使用不同的过滤条件来遍历QC指标。它不一定是线性过程。满足过滤条件后,我们将保存过滤后的细胞对象以进行聚类和标记识别。
# Create .RData object to load at any time
save(filtered_seurat, file="data/seurat_filtered.RData")
注意:我们还有其他材料可用于探索质量较差的样本的指标。
本课程由哈佛大学生物信息学核心(HBC)的教学团队成员开发。这些是根据知识共享署名许可(CC BY 4.0)的条款分发的开放获取材料,只要原始作者和出处均受到认可,就可以在任何介质中不受限制地使用,分发和复制。
######################################
####04_SC_quality_control.md
# Explore merged metadata
View(merged_seurat@meta.data)
# Add number of genes per UMI for each cell to metadata
merged_seurat$log10GenesPerUMI <- log10(merged_seurat$nFeature_RNA) / log10(merged_seurat$nCount_RNA)
# Compute percent mito ratio
merged_seurat$mitoRatio <- PercentageFeatureSet(object = merged_seurat, pattern = "^MT-")
merged_seurat$mitoRatio <- merged_seurat@meta.data$mitoRatio / 100
# Create metadata dataframe
metadata <- merged_seurat@meta.data
# Add cell IDs to metadata
metadata$cells <- rownames(metadata)
# Rename columns
metadata <- metadata %>%
dplyr::rename(seq_folder = orig.ident,
nUMI = nCount_RNA,
nGene = nFeature_RNA)
# Create sample column
metadata$sample <- NA
metadata$sample[which(str_detect(metadata$cells, "^ctrl_"))] <- "ctrl"
metadata$sample[which(str_detect(metadata$cells, "^stim_"))] <- "stim"
# 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")
# 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")
# 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)
# 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")
# 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)
# 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)
# 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)
# Filter out low quality reads 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))
dim(filtered_seurat)
# Output a logical vector for every gene on whether the more than zero counts per cell
# Extract counts
counts <- GetAssayData(object = filtered_seurat, slot = "counts")
# Output a logical vector for every gene on whether the more than zero counts per cell
nonzero <- counts > 0
# 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, ]
# Reassign to filtered Seurat object
filtered_seurat <- CreateSeuratObject(filtered_counts, meta.data = filtered_seurat@meta.data)
# Save filtered subset to new metadata
metadata_clean <- filtered_seurat@meta.data
# Create .RData object to load at any time
save(filtered_seurat, file="data/seurat_filtered.RData")