clustree ： 聚类可视化利器

我们知道在研究问题时，分组是很重要的，有分组才有故事可讲。比如，两块田一块施肥一块不施肥，可以做比较嘛。在单细胞数据分析中用到较多的数据分组技术是聚类（clustering），这里面有很多的喜怒哀乐，因为聚类是无监督的，而且可以聚成不同的层次，在第一次聚类后，又可以对亚群聚类，真是子子孙孙无穷匮也。这也是单细胞数据分析的魅力所在：不同层次的聚类就像剥洋葱，剥着剥着，说不定就泪流满脸了呢？

如：

Global characterization of T cells in non-small-cell lung cancer by single-cell sequencing

这T细胞分的好细致啊。

从目前大部分的单细胞分析流程来看，单细胞数据分析属于探索性数据分析。因为他用的主要分析方法，降维和聚类，都属于探索性数据分析方法。这也解释了为什么目前大部分的单细胞文章都是描述性质的了。就是把我们所能看到的东西报告给国际同行。一群细胞，有的人能看出十个类（找不同），有的人只能看出两点不同。这里的十个和两点就相当于我们说的resolution.。

探索，对有的人来说是好事；对另一些人来说是坏事，他们需要确定的东西。探索的时候，多表现得贪婪：你帮我做一下看看，我想看看有什么结果，这是正常的，探索嘛。但是数据探索不是无止境的，虽然探索的时候尽量不带有偏见，但是探索的过程就是形成偏见的过程：指导下游分析。

鉴于单细胞数据探索性的特点，这启示我们：

首先，把数据分析流程分为两个阶段：探索与验证。探索时贪婪，多试试；验证时谨慎，慢慢来。二者要二八开，不要在探索中沉浸太久，探索的过程就像走马观花，获得数据概览，总体概貌。这也是不需要学科背景的一部分，所以哪家公司都可以做。他们的流程也就叫做单细胞探索性数据分析流程。而后面的百分之八十，是需要学科背景的，基本上离开课题组哪家公司都不能做，所以不要问我哪家公司单细胞做得好。这就像，你帮生病的住在的国外朋友找医院一样，医生如果不能对病人望闻问切，只能给你走一个标准体检流程。

君不见，有的单细胞数据就像不断转诊的病人，希望找到可以治愈焦虑的医生。因为连探索都没有遇到好的公司，更别说验证了。

第二，对科技服务公司来说，要找准自己的定位：优化探索流程，以提供研究视角。探索性数据分析是可流程化的，数据的分布，质量的检查，降维聚类，这些可以提供数据概览，让客户一看就可以有个切入点。验证这一块，可以配一个专业的团队和客户一起，在学科背景加持下展开分析。总结验证分析的一般原则，转化到探索流程中去。

其实探索性数据分析是一门学问。

探索性数据分析（Exploratory Data Analysis，简称EDA），探索性数据分析是上世纪六十年代从数据分析中提炼出来的，其方法由美国统计学家John Tukey提出。是指在尽量少的先验假定下进行探索，通过作图（可视化）、制表（统计细胞数）、计算特征量（降维），聚类（发现类）等手段探索数据的结构（群）和规律（轨迹）的一种数据分析方法。特别是单细胞数据以高维多高噪声为显著特点，拿到数据后往往不知所措，不知道从哪里开始了解目前拿到手上的数据（excel看不了，或看不全），这时候探索性数据分析就非常有效。

聚类技术广泛应用于大型数据集的分析，将具有相似性质的样本聚类在一起。例如，聚类常用于单细胞rna测序领域，以识别组织样本中存在的不同细胞类型。执行聚类的算法有很多，结果可能有很大差异。特别是，数据集中出现的组的数量通常是未知的，而算法识别的集群数量可以根据所使用的参数而改变。为了探讨和检验不同聚类分辨率的影响，我们使用聚类树（clustree ）可视化显示在多个分辨率下分群之间的关系，允许研究人员看到样本如何随着分群数量的增加而移动。此外，元信息可以覆盖在树中，以告知解决方案的选择，并指导分群的识别。

同样请出我们的数据集：

library(Seurat)
library(SeuratData)
library(clustree)

pbmc3k.final

An object of class Seurat 
13714 features across 2638 samples within 1 assay 
Active assay: RNA (13714 features, 2000 variable features)
 2 dimensional reductions calculated: pca, umap

到底分多少个群是合适的呢？为了规避这个问题，我们多做几次聚类吧。

pbmc3k.final <- FindClusters(
    object = pbmc3k.final,
    resolution = c(seq(0,1.6,.2))
)

用clustree 看看不同resolution 下的相互关系：

clustree(pbmc3k.final@meta.data, prefix = "RNA_snn_res.")

查一下clustree是何方神圣。

原来为Seurat 提供了接口：

clustree(pbmc3k.final, prefix = "RNA_snn_res.", node_colour = "purple", node_size = 10,
         node_alpha = 0.8) # 可以直接输入Seurat对象

然后这是一款专门为探索聚类分析而开发的R包，果然厉害。

我们可以在这个tree上绘制更多信息。

如线粒体信息啦，以及其他你想看的分数值。

clustree(pbmc3k.final@meta.data, prefix = "RNA_snn_res.", node_colour = "percent.mt",
         node_colour_aggr = "mean")

我们可以调整布局：

clustree(pbmc3k.final@meta.data, prefix = "RNA_snn_res.", layout = "sugiyama")

增加节点信息：

clustree(pbmc3k.final, prefix = "RNA_snn_res.", node_label = "RNA_snn_res.")

标出线粒体含量：

clustree(pbmc3k.final, prefix = "RNA_snn_res.", node_label = "percent.mt",node_label_aggr = "median")

如果有先验知识，看看分群与之的匹配程度：

label_position <- function(labels) {
    if (length(unique(labels)) == 1) {
        position <- as.character(unique(labels))
    } else {
        position <- "mixed"
    }
    return(position)
}

clustree(pbmc3k.final@meta.data, prefix = "RNA_snn_res.", node_label = "seurat_annotations",
         node_label_aggr = "label_position")

那些mixed的群可以考察一下，是不是可以再分。

绘制表达量。

clustree(pbmc3k.final, prefix = "RNA_snn_res.",
         node_colour = "CD79A", node_colour_aggr = "median")

这个可以看出clustree对Seurat的支持力度了。

我们来在umap空间绘制不同resolution 的分布情况。

pbmc3k.final<- AddMetaData(pbmc3k.final,pbmc3k.final@reductions$umap@cell.embeddings,
                           col.name = c("UMAP_1","UMAP_2"))
pbmc3k.final<- AddMetaData(pbmc3k.final,pbmc3k.final@reductions$pca@cell.embeddings,
                           col.name = colnames(pbmc3k.final@reductions$pca@cell.embeddings))


clustree_overlay(pbmc3k.final, prefix = "RNA_snn_res.", x_value = "UMAP_1", y_value = "UMAP_2")

clustree_overlay(pbmc3k.final, prefix = "RNA_snn_res.", x_value = "UMAP_1", y_value = "UMAP_2",
                 use_colour = "points", alt_colour = "blue") 

加上标签

clustree_overlay(pbmc3k.final, prefix = "RNA_snn_res.", x_value = "UMAP_1", y_value = "UMAP_2",label_nodes = TRUE)

overlay_list <- clustree_overlay(pbmc3k.final, prefix = "RNA_snn_res.", x_value = "PC_1", y_value = "PC_2",
                                 plot_sides = TRUE)

names(overlay_list)
#> [1] "overlay" "x_side"  "y_side"

overlay_list$x_side

overlay_list$y_side

同时支持SingleCellExperiment格式的数据。

suppressPackageStartupMessages(library("SingleCellExperiment"))
data("sc_example")
names(sc_example)

sce <- SingleCellExperiment(assays = list(counts = sc_example$counts,
                                          logcounts = sc_example$logcounts),
                            colData = sc_example$sc3_clusters,
                            reducedDims = SimpleList(TSNE = sc_example$tsne))

head(colData(sce))

clustree(sce, prefix = "sc3_", suffix = "_clusters")

我们还可以用ggplot语法来修饰。

clustree(pbmc3k.final, prefix = "RNA_snn_res.") +
    guides(edge_colour = FALSE, edge_alpha = FALSE) +scale_color_brewer(palette = "Set1") +
    scale_edge_color_continuous(low = "blue", high = "red")+
    theme(legend.position = "bottom")