R中单细胞RNA-seq分析教程 (4)

数据科学工厂

引言

本系列开启R中单细胞RNA-seq数据分析教程，持续更新，欢迎关注，转发！

8. 细胞聚类

分析 scRNA-seq 数据时，绘制标记基因的特征图通常是一个良好的起点。然而，要深入理解数据中的异质性，需要通过无偏的方式将细胞分组，这就是聚类的作用。理论上，任何聚类方法都可以应用于 scRNA-seq 数据，包括层次聚类和 k-means 等常用于 bulk RNA-seq 数据的方法。但由于 scRNA-seq 数据样本量通常极大（如一次 10x 实验可能包含数千个细胞），这些方法运行速度非常慢。此外，由于 scRNA-seq 数据本身的稀疏性，即便通过 PCA 等降维处理去噪，不同细胞之间的差异也难以像 bulk RNA-seq 那样精准量化。因此，scRNA-seq 数据分析中更常采用基于图的社区识别算法。这里的“图”是一个数学概念，表示一组对象及其之间的关系；通俗地说，就是细胞间构建的网络。

具体来说，首先会生成一个细胞的 k-最近邻（KNN）网络。每个细胞根据其主成分（PC）值，与距离最近的细胞建立连接，只有互为邻居的细胞才会被视为连接对。接着，计算每对细胞之间共享邻居的比例，用以表示两细胞之间连接的强度，并将较弱的连接剪除。这便形成了最终的共享最近邻（SNN）网络。

seurat <- FindNeighbors(seurat, dims = 1:20)

构建好网络后，可以使用 Louvain 社区识别算法来划分网络中的社区，也就是细胞群体。在同一个群体内，细胞之间的连接较为密集，而不同群体之间的连接则相对稀疏。

seurat <- FindClusters(seurat, resolution = 1)

resolution 参数用来决定聚类结果的颗粒度：是返回较大的细胞群体（如主要细胞类型），还是更小、更精细的细胞群（如细胞亚型）。分辨率一般设置在 0.1 到 1 之间，具体选择取决于分析目标。在这里，使用较高的分辨率以实现更细致的聚类。可以多次运行 FindClusters 函数，尝试不同的分辨率值。最新的聚类结果可以通过 Idents(seurat) 或 [email]seurat@active.ident[/email] 获取，而之前的聚类结果则保存在 meta.data 槽（[email]seurat@meta.data[/email]）中的不同列里。

然后，可以利用之前生成的 tSNE 和 UMAP 嵌入来可视化聚类结果。

plot1 <- DimPlot(seurat, reduction = "tsne", label = TRUE)
plot2 <- DimPlot(seurat, reduction = "umap", label = TRUE)
plot1 + plot2

9. 为细胞簇添加注释

给细胞分组能为每个细胞分配一个识别标签，通常认为拥有相同标签的细胞在特性上是相似的，因此可以归为同一类型的细胞或细胞状态。但要确定这些细胞簇具体代表哪种细胞类型或状态，并非易事，通常也没有绝对正确的答案。解决这个问题有几种策略可以尝试，比如：

1. 检视这些细胞簇中典型细胞类型和状态的标记基因表达情况；
2. 找出每个已辨识细胞簇的标志性或标记性基因。通过这些标记基因，可以查阅文献、进行富集分析或进行实验（或咨询同行）来确定细胞簇的类型；
3. 将每个细胞簇的基因表达图谱与已知的参考数据进行对比。

显然，第一种方法需要对研究系统有一定的了解。研究者需要掌握一份被广泛认可的标记基因清单。对于所举例的数据集（大脑类器官）而言，一些标记基因已经提及。一个更长的列表（永远不可能完整）包括：

• NES / SOX2：神经前体细胞的标志
• NHLH1：新生神经元的标志
• DCX / MAP2 / MAPT：神经元的标志
• FOXG1：大脑皮层的标志
• EMX1 / EMX2：大脑皮层背侧的标志

要直观展示细胞簇中标记基因的表达情况，热图可能是最简便的方法。

ct_markers <- c("MKI67","NES","DCX","FOXG1", # G2M, NPC, neuron, telencephalon
                "DLX2","DLX5","ISL1","SIX3","NKX2.1","SOX6","NR2F2", # ventral telencephalon related
                "EMX1","PAX6","GLI3","EOMES","NEUROD6", # dorsal telencephalon related
                "RSPO3","OTX2","LHX9","TFAP2A","RELN","HOXB2","HOXB5") # non-telencephalon related
DoHeatmap(seurat, features = ct_markers) + NoLegend()

接下来，为了使细胞注释过程更加客观，首先需要为每个已识别的细胞簇确定其特定的标记基因。在Seurat软件包中，可以使用 FindAllMarkers功能来完成这一任务。这个功能的作用是针对特定的细胞簇，通过Wilcoxon秩和检验方法，分析该簇内的细胞与簇外细胞在基因表达上的差异。

cl_markers <- FindAllMarkers(seurat, only.pos = TRUE, min.pct = 0.25, logfc.threshold = log(1.2))
library(dplyr)
cl_markers %>% group_by(cluster) %>% top_n(n = 2, wt = avg_logFC)

鉴于单细胞RNA测序（scRNA-seq）数据中样本量巨大（每个细胞作为一个样本），在分析时，不仅需要关注p值，还应该考虑基因在特定细胞簇中的检测频率（pct）以及该基因在簇内外细胞表达的相对变化量（logfc）。因此，相关函数提供了 min.pct和 logfc.threshold参数，用以设定效应大小的阈值。

你可能已经意识到这个过程耗时较长。实际上，有一个名为“presto”的包提供了更快速的解决方案。

library(presto)
cl_markers_presto <- wilcoxauc(seurat)
cl_markers_presto %>%
    filter(logFC > log(1.2) & pct_in > 20 & padj < 0.05) %>%
    group_by(group) %>%
    arrange(desc(logFC), .by_group=T) %>%
    top_n(n = 2, wt = logFC) %>%
    print(n = 40, width = Inf)

presto 的输出结果与 Seurat 的原生方法类似，但提供了更多的统计指标。

无论选择哪种方法，最终鉴定出的主要聚类标记基因都可以通过绘制热图来直观展示。

top10_cl_markers <- cl_markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
DoHeatmap(seurat, features = top10_cl_markers$gene) + NoLegend()

可以通过特征图或小提琴图更详细地查看不同细胞群的标记基因。例如，NEUROD2 和 NEUROD6 是第 2 群的显著标记基因，可以进一步分析它们的表达情况。

plot1 <- FeaturePlot(seurat, c("NEUROD2","NEUROD6"), ncol = 1)
plot2 <- VlnPlot(seurat, features = c("NEUROD2","NEUROD6"), pt.size = 0)
plot1 + plot2 + plot_layout(widths = c(1, 2))

NEUROD2 和 NEUROD6 不仅在聚类 2 中高度表达，而且在聚类 6 中也有较高的表达。如果你记得这些聚类在 tSNE/UMAP 图中的位置，你会发现聚类 2 和聚类 6 紧邻在一起，表明这两个聚类可能代表了相互关联的细胞类型，它们都具有强烈的背侧端脑特征。它们的分离可能反映了它们的成熟状态。聚类 6 中的神经元可能较为幼稚，因为它们与聚类 0 相关，而聚类 0 很可能是背侧端脑的神经祖细胞（NPCs）。此外，聚类 6 表达了高水平的 EOMES，这是一种背侧端脑中间祖细胞（IP）的标记基因。综合分析，可以有较高的信心认为，聚类 0、6 和 2 都是背侧端脑细胞，聚类 0 是祖细胞，聚类 6 是中间祖细胞，聚类 2 是神经元。

如果看聚类 0 另一侧，会发现它与聚类 5 相连，接着是聚类 10。再次查看热图后，可以看到这些聚类虽然具有各自独特的标记基因和表达模式，但它们都呈现出典型的背侧端脑神经祖细胞特征。此外，聚类 5 和 10 的细胞表现出较高的细胞周期 G2M 相位标记基因的表达，这表明它们也是背侧端脑的神经祖细胞，它们与聚类 0 的分离很可能与细胞周期的不同阶段有关。

有趣的是，聚类 10、5、0、6 和 2 中的细胞在 UMAP 图上形成了一个类似轨迹的结构，这可能反映了从神经祖细胞到神经元成熟的分化过程。稍后会进一步探讨这个问题。

这就是细胞聚类注释的常见方法。你可能觉得这种方法过于主观，依赖个人判断。如果是这样，也有更客观、无偏的方式来进行自动化或半自动化注释。目前有一些新的工具，比如 Cole Trapnell 实验室开发的 Garnett 和 Samantha Morris 实验室开发的 Capybara。这些工具的工作原理类似：它们首先对现有的 scRNA-seq 数据进行细胞类型注释的标准化，利用这些数据训练预测模型，然后将模型应用到新的数据集进行自动注释。然而，这些工具目前有一定局限性，它们通常只适用于主要细胞类型的注释，且其表现很大程度上取决于训练数据集的质量。

值得注意的是，细胞聚类注释并不总是需要复杂的机器学习模型。计算 scRNA-seq 数据中细胞或细胞群体的基因表达谱与已知数据（如原料数据）的相关性，也可以提供很有价值的信息。例如，团队开发的 VoxHunt 工具，通过将细胞或细胞群体的表达谱与 Allen Brain Atlas 中小鼠大脑发育过程的原位杂交图谱相关联，来帮助注释。这个方法对于注释大脑类器官样本的 scRNA-seq 数据尤其有用。

library(voxhunt)
load_aba_data('ABA_data')
genes_use <- variable_genes('E13', 300)$gene
vox_map <- voxel_map(seurat, genes_use=genes_use)
plot_map(vox_map)

从投影结果来看，也可以得出相似的结论：聚类 10、5、0、6 和 2 都属于背侧端脑。最后，可以对所有聚类进行初步的注释。

可以选择用注释替换细胞聚类标签，但这一步是可选的。

new_ident <- setNames(c("Dorsal telen. NPC",
                        "Midbrain-hindbrain boundary neuron",
                        "Dorsal telen. neuron",
                        "Dien. and midbrain excitatory neuron",
                        "MGE-like neuron","G2M dorsal telen. NPC",
                        "Dorsal telen. IP","Dien. and midbrain NPC",
                        "Dien. and midbrain IP and excitatory early neuron",
                        "G2M Dien. and midbrain NPC",
                        "G2M dorsal telen. NPC",
                        "Dien. and midbrain inhibitory neuron",
                        "Dien. and midbrain IP and early inhibitory neuron",
                        "Ventral telen. neuron",
                        "Unknown 1",
                        "Unknown 2"),
                      levels(seurat))
seurat <- RenameIdents(seurat, new_ident)
DimPlot(seurat, reduction = "umap", label = TRUE) + NoLegend()

微信扫一扫分享文章