scRNA-seq Clustering quality control
scRNA-seq Clustering quality control
- 评估是否存在clustering artifacts
- 使用PCA，tSNE和UMAP图确定聚类的质量，并了解何时重新聚类
- 评估已知的细胞类型标记以假设群集的细胞类型同一性
目标
- **生成特定于细胞类型的簇，**并使用已知的标记来确定簇的身份。
- 确定分群是否代表真实的细胞类型或由于生物或技术差异而形成的群集 ，例如处于细胞周期S期的细胞群集、特定批次的群集或高线粒体含量的细胞。
挑战
- 识别每个分群的 细胞类型
- 保持耐心，因为这可能是聚类和标记识别之间高度迭代的过程(有时甚至回到QC过滤或标准化)
建议
- 对存在的细胞类型和这些细胞类型的几个标记基因有一个很好的预期。了解你所期望的是低复杂性的细胞类型还是高线粒体含量的细胞类型，以及这些细胞是否正在分化。
- 确定要删除的任何垃圾群集。垃圾群集可能是包含高线粒体含量和低UMIs/基因的那些。
- 如果 未将所有细胞类型检测为单独的群 ，请尝试更改UMAP分辨率、用于分群的PC数量或使用的可变基因数量
探讨质控指标
为了确定我们的分群是否可能是由于细胞周期阶段或线粒体表达等人工因素造成的，可视化探索这些指标以查看是否有任何簇表现出富集或与其他簇不同，这是很有用的。然而，如果观察到特定簇的富集或差异，它可以用细胞类型来解释，那就可以不必担忧。
为了探索和可视化各种质量指标，我们将使用来自Seurat的通用的 DimPlot() 和 FeaturePlot() 函数。
按样本划分群集
通过样本分离群集，我们可以从探索每个样本中每个群集的细胞分布开始：
# Extract identity and sample information from seurat object to determine the number of cells per cluster per sample
n_cells <- FetchData(seurat_integrated,
vars = c("ident", "orig.ident")) %>%
dplyr::count(ident, orig.ident) %>%
tidyr::spread(ident, n)
# View table
View(n_cells)
我们可以使用UMAP可视化每个样本的每个群集的细胞：
# UMAP of cells in each cluster by sampleDimPlot(seurat_integrated,
label = TRUE,
split.by = "sample") + NoLegend()
按细胞周期阶段划分群集
接下来，我们将探讨细胞是否会因不同的细胞周期阶段聚集。当我们对无意义的变异源进行SCTransform归一化和回归时，并没有因为细胞周期阶段而使变异消退。如果我们的细胞簇在线粒体表达上表现出很大的差异，这预示着我们要重新运行SCTransform，并将 S.Score 和 G2M.Score 添加到我们的变量中以进行回归，然后重新运行其余步骤。
# Explore whether clusters segregate by cell cycle phase
DimPlot(seurat_integrated,
label = TRUE,
split.by = "Phase") + NoLegend()
根据细胞周期划分，我们没有看到太多的聚集，所以我们可以继续进行QC。
按各种无意义的变异源进行的群集分离
接下来，我们将探索其他指标，例如每个细胞的UMI和基因数量，S期和G2M期标记，以及通过UMAP进行的线粒体基因表达。查看各个S和G2M分数可以为我们提供更多信息来检查细胞周期因素的影响，就像我们之前所做的那样。
# Determine metrics to plot present in seurat_integrated@meta.data
metrics <- c("nUMI", "nGene", "S.Score", "G2M.Score", "mitoRatio")
FeaturePlot(seurat_integrated,
reduction = "umap",
features = metrics,
pt.size = 0.4,
sort.cell = TRUE,
min.cutoff = 'q10',
label = TRUE)
注意：
sort.cell参数将在低表达细胞上方绘制高表达细胞，而min.cutoff参数将确定着色的阈值。min.cutoff的q10意味着基因表达最低的10%的细胞，不会表现出任何紫色阴影（完全灰色）
指标在整个群集中相对均匀，但nUMI和nGene在群集3、9、14和15中显示出更高的值，也许还包括群集17。我们将密切关注这些群集，看看细胞类型是否可以解释这种增加。
按PCs分群的探索
我们还可以探索群集通过不同的PC分开的情况。我们希望定义的PC能够很好地区分细胞类型。要可视化此信息，我们需要提取细胞的UMAP坐标信息及其对应的分数，以便通过UMAP查看每个PC。
首先，我们确定要从Seurat对象提取的信息，然后，可以使用 FetchData() 函数提取它。
# Defining the information in the seurat object of interest
columns <- c(paste0("PC_", 1:16),
"ident",
"UMAP_1", "UMAP_2")
# Extracting this data from the seurat object
pc_data <- FetchData(seurat_integrated,
vars = columns)
注意：我们如何在知道
FetchData()函数能提取UMAP_1以获得UMAP坐标？Seurat速查表(https://satijalab.org/seurat/essential_commands.html)描述了该函数能够从表达式矩阵，细胞嵌入或元数据中提取任何数据。例如，如果您浏览seurat_Integrated@Reductions列表对象，第一个组件是用于PCA的，并且包括一个用于cell.embedding的槽。我们可以使用列名(PC_1、PC_2、PC_3等)。以调出对应于每个PC的每个细胞的坐标或PC分数。 我们可以对UMAP做同样的事情：
# Extract the UMAP coordinates for the first 10 cells
seurat_integrated@reductions$umap@cell.embeddings[1:10, 1:2]
FetchData()函数仅使我们能更轻松地提取数据。
在下面的UMAP图中，细胞根据各自主成分的PC得分进行着色。
让我们快速了解一下前16个PC：
# Adding cluster label to center of cluster on UMAP
umap_label <- FetchData(seurat_integrated,
vars = c("ident", "UMAP_1", "UMAP_2")) %>%
group_by(ident) %>%
summarise(x=mean(UMAP_1), y=mean(UMAP_2))
# Plotting a UMAP plot for each of the PCs
map(paste0("PC_", 1:16), function(pc){
ggplot(pc_data,
aes(UMAP_1, UMAP_2)) +
geom_point(aes_string(color=pc),
alpha = 0.7) +
scale_color_gradient(guide = FALSE,
low = "grey90",
high = "blue") +
geom_text(data=umap_label,
aes(label=ident, x, y)) +
ggtitle(pc)
}) %>%
plot_grid(plotlist = .)
我们可以看到在不同的PC下的群集展示。例如，驱动PC_2基因在簇6、11和17中表现出更高的表达（在15中也可能更高）。我们可以回顾一下驱动此PC的基因，以了解细胞类型可能是什么：
# Examine PCA results
print(seurat_integrated[["pca"]], dims = 1:5, nfeatures = 5)
以CD79A基因和HLA基因作为PC_2的阳性标志物，可以推测第6、11、17簇对应于B细胞。这只是暗示了集群身份可能是什么，集群的身份是通过PC的组合来确定的。
为了真正确定群集的身份和 resolution 是否合适，探索几个预期细胞类型的已知标记是很有帮助的。
未完待续......
注：以上内容来自哈佛大学生物信息中心(HBC)_的教学团队的生物信息学培训课程。原文链接：https://hbctraining.github.io/scRNA-seq/schedule/ 点击 “阅读原文” 可直达
