# 3.5 聚类

## 1 前言

聚类是一种无监督的学习过程，在scRNA分析中，根据表达谱相似性对不同细胞进行分组。

在后面根据marker基因对各个cluster进行注释后，cluster就可以代表不同的细胞类型或状态。因此，聚类是将枯燥的数据与感兴趣生物学意义联系起来的一个桥梁。

不过这里需要明确一下，聚类分的群和细胞类型还不是一回事

* 分群是根据经验计算的结果
* 细胞类型有真正的生物学意义

至于这种问题：“我应该设置多少cluster数量？”，则是个伪命题。因为我们可以根据需要，使用不同方法，得到不同数量的clusters，每个cluster都是高维空间数据的一部分。

不过倒是可以问：“我们的clusters与细胞类型之间的拟合程度如何？”。不过这个问题很难回答，因为取决于自己的分析目标。有的人对于能得到主要的细胞类型就很满意，而有的人还想再进一步，得到细胞亚群；还有的，更想基于细胞亚群，达到不同细胞状态的分辨率（如代谢活性、胁迫反应等）。另外，两种截然不同的分群结果，可能都具有生物学意义，只是关注点不同而已。

对于分群聚类的操作，更像是显微镜的操作，可以通过调整参数（使用不同的镜头组合）来获得不同的分辨率；另外还可以探索多种聚类方法，看到数据的不同面。

### **数据准备之PBMC**

具体的操作在上一章也提到过

```r
# 下载
library(BiocFileCache)
bfc <- BiocFileCache("raw_data", ask = FALSE)
raw.path <- bfcrpath(bfc, file.path("http://cf.10xgenomics.com/samples",
    "cell-exp/2.1.0/pbmc4k/pbmc4k_raw_gene_bc_matrices.tar.gz"))
untar(raw.path, exdir=file.path(tempdir(), "pbmc4k"))

# 读取
suppressMessages(library(DropletUtils))
fname <- file.path(tempdir(), "pbmc4k/raw_gene_bc_matrices/GRCh38")
sce.pbmc <- read10xCounts(fname, col.names=TRUE)
> dim(sce.pbmc)
[1]  33694 737280

# 基因注释之ID整合
library(scater)
rownames(sce.pbmc) <- uniquifyFeatureNames(
    rowData(sce.pbmc)$ID, rowData(sce.pbmc)$Symbol)
# 基因注释之添加位置信息
library(EnsDb.Hsapiens.v86)
location <- mapIds(EnsDb.Hsapiens.v86, keys=rowData(sce.pbmc)$ID, 
    column="SEQNAME", keytype="GENEID")

# 空液滴检测
set.seed(100)
e.out <- emptyDrops(counts(sce.pbmc))
sce.pbmc <- sce.pbmc[,which(e.out$FDR <= 0.001)]

# 质控
stats <- perCellQCMetrics(sce.pbmc, subsets=list(Mito=which(location=="MT")))
high.mito <- isOutlier(stats$subsets_Mito_percent, type="higher")
sce.pbmc <- sce.pbmc[,!high.mito]

# 归一化（去卷积）
library(scran)
set.seed(1000)
clusters <- quickCluster(sce.pbmc)
sce.pbmc <- computeSumFactors(sce.pbmc, cluster=clusters)
sce.pbmc <- logNormCounts(sce.pbmc)

# 利用数据分布来表示技术噪音，并挑选HVGs
set.seed(1001)
dec.pbmc <- modelGeneVarByPoisson(sce.pbmc)
top.pbmc <- getTopHVGs(dec.pbmc, prop=0.1)

# 降维（三种方法）
set.seed(10000)
sce.pbmc <- denoisePCA(sce.pbmc, subset.row=top.pbmc, technical=dec.pbmc)

set.seed(100000)
sce.pbmc <- runTSNE(sce.pbmc, dimred="PCA")

set.seed(1000000)
sce.pbmc <- runUMAP(sce.pbmc, dimred="PCA")

# 看下结果
sce.pbmc
## class: SingleCellExperiment 
## dim: 33694 3922 
## metadata(1): Samples
## assays(2): counts logcounts
## rownames(33694): RP11-34P13.3 FAM138A ... AC213203.1 FAM231B
## rowData names(2): ID Symbol
## colnames(3922): AAACCTGAGAAGGCCT-1 AAACCTGAGACAGACC-1 ...
##   TTTGTCACAGGTCCAC-1 TTTGTCATCCCAAGAT-1
## colData names(3): Sample Barcode sizeFactor
## reducedDimNames(3): PCA TSNE UMAP
## altExpNames(0):
```

**补充知识点之聚类与分类**

* 分类：类别是已知的，通过对已知分类的数据进行训练和学习，找到这些不同类的特征，再对未分类的数据进行分类。属于有监督学习。
* 聚类：事先不知道数据会分为几类，通过聚类分析将数据聚合成几个群体。聚类不需要对数据进行训练和学习。**属于无监督学习**。

## 2 基于图形的聚类 |  Graph-based clustering

### **2.1 是什么？**

这种聚类方法在**Seurat中使用最流行**。最基础的想法是：我们首先构建一个图，其中每个节点都是一个细胞，它与高维空间中最邻近的细胞相连。连线基于细胞之间的相似性计算权重，权重越高，表示细胞间关系更密切。

如果一群细胞之间的权重高于另一群细胞，那么这一群细胞就被当做一个群体 “communiity”

这个方法**最大的优点就是：可缩放性**。只需要根据 k-nearest neighbor（KNN）进行搜索，自己去寻找去拓展。相比于k-means、 Gaussian mixture models等方法，不需要预先对cluster形状以及cluster中细胞分布做一个估计。

这种方法**使得每个细胞都被强制连接到一定数量的相邻细胞**，这减少了仅由一两个离群细胞组成的无意义群体的风险

不过这个方法也有缺点：

**缺点一：**

它最后只保留了邻近细胞之间的联系，除此以外没有保存。

**缺点二：**

另外正是由于它主要关注相邻的权重，导致一些噪音也可能聚集在一起，形成一个community，从而影响聚类的分辨率，可能导致数据异质性的过度解读

### **2.2 怎么做？**

在使用时，需要考虑几个问题：

* 构建这个图需要设置几个邻居？
* 用什么方法确定边的权重？
* 用什么检测community的方法来定义cluster？

下面的例子中，就定义了使用一个细胞与10个邻居的关系，进而构建shared nearest neighbor graph（SNNG），如果两个细胞的邻居中有一个是共享的，则通过一条边连接，之后再根据highest average rank of the shared neighbors定义边的权重([Xu and Su 2015](https://pubmed.ncbi.nlm.nih.gov/25805722/))

**首先使用scran包的`buildSNNGraph()`函数，并且使用PCA的前几个PCs；之后利用 igraph包中的Walktrap方法去鉴定clusters**

```r
library(scran)
g <- buildSNNGraph(sce.pbmc, k=10, use.dimred = 'PCA')
clust <- igraph::cluster_walktrap(g)$membership
table(clust)
## clust
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18 
## 585 518 364 458 170 791 295 107  45  46 152  84  40  60 142  16  28  21
```

之后可以把cluster信息放回到`SingleCellExperiment`对象中，保存成一个因子

```r
library(scater)
colLabels(sce.pbmc) <- factor(clust)
plotReducedDim(sce.pbmc, "TSNE", colour_by="label")
```

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-06-30-065754.png)

从上面的计算方法可以看出，有一个很重要的参数是`k` ，含义是：the number of nearest neighbors used to construct the graph。如果k设置越大，得到的图之间联通程度越高，cluster也越大。因此这个参数也是可以不断尝试的

```r
# 比如设置较大的k，就会得到比k=10更少的cluster，但同时每个cluster平均的细胞数更多
g.50 <- buildSNNGraph(sce.pbmc, k=50, use.dimred = 'PCA')
clust.50 <- igraph::cluster_walktrap(g.50)$membership
table(clust.50)
## clust.50
##   1   2   3   4   5   6   7   8 
## 307 729 789 187 516 524 825  45

# 如果设置更小的k，会得到数量更多的cluster
g.5 <- buildSNNGraph(sce.pbmc, k=5, use.dimred = 'PCA')
clust.5 <- igraph::cluster_walktrap(g.5)$membership
table(clust.5)
## clust.5
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
##  81  45 457 296 168 350  79 432 428 897  64 171  68 135  79  26  18  30  21  16 
##  21  22  23 
##  36   9  16
```

这个结果可以利用`force-directed`的布局出图，它与t-SNE、UMAP的降维结果都是紧密相关的

```r
set.seed(2000)
reducedDim(sce.pbmc, "force") <- igraph::layout_with_fr(g)
plotReducedDim(sce.pbmc, colour_by="label", dimred="force")
```

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-06-30-070536.png)

### **2.3 其他的参数**

**关于如何计算权重**，在`buildSNNGraph()`函数中可以设置：

* `type="number"` ：根据近邻数量计算
* `type="jaccard"` ：根据Jaccard index of the two sets of neighbors计算

```r
g.num <- buildSNNGraph(sce.pbmc, use.dimred="PCA", type="number")
g.jaccard <- buildSNNGraph(sce.pbmc, use.dimred="PCA", type="jaccard")
```

上面得到的各种结果，都是来自igraph包的graph对象，因此可以用于igraph包的**各种community检测方法：**

```r
# 之前使用的是Walktrap聚类方法，还有这么多种
clust.louvain <- igraph::cluster_louvain(g)$membership
clust.infomap <- igraph::cluster_infomap(g)$membership
clust.fast <- igraph::cluster_fast_greedy(g)$membership
clust.labprop <- igraph::cluster_label_prop(g)$membership
clust.eigen <- igraph::cluster_leading_eigen(g)$membership
```

最后就可以**比较两种不同的聚类方法差异**

```r
library(pheatmap)

# 比较Infomap 和 Walktrap
# Using a large pseudo-count for a smoother color transition
# between 0 and 1 cell in each 'tab'.
tab <- table(paste("Infomap", clust.infomap), 
    paste("Walktrap", clust))
ivw <- pheatmap(log10(tab+10), main="Infomap vs Walktrap",
    color=viridis::viridis(100), silent=TRUE)

# Fast-greedy 和 Walktrap
tab <- table(paste("Fast", clust.fast), 
    paste("Walktrap", clust))
fvw <- pheatmap(log10(tab+10), main="Fast-greedy vs Walktrap",
    color=viridis::viridis(100), silent=TRUE)

gridExtra::grid.arrange(ivw[[4]], fvw[[4]])
```

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-06-30-071422.png)

可以看出，Infomap比Walktrap得到了更多的cluster（比较精细），而fast-greedy结果相对较少（比较粗糙）

**可以通过`cut_at()`指定分群的数量**

```r
# 指定5群
community.walktrap <- igraph::cluster_walktrap(g)
table(igraph::cut_at(community.walktrap, n=5))
## 
##    1    2    3    4    5 
## 3612  198   45   46   21

# 指定20群
table(igraph::cut_at(community.walktrap, n=20))
## 
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
## 533 364 458 442 170 791 295 107  45  46 152  84  40  60 142  76  52  16  28  21
```

scran默认使用rank-based的权重计算方法，Seurat使用Jaccard-based的方法（后面辅助Louvain聚类）

### **2.4 分群的评价**

当使用graph-based方法时，模块化（modularity）是一个评价cluster的指标。**modularity值越大，表示cluster内部的细胞与细胞之间连线（edge）数最多，内部分散效果越好**，从而避免了不同cluster之间邻近的细胞形成连线。

之前通过`cluster_walktrap`的方法得到了18个cluster，下面就看看这个方法做的如何：会用到一个函数`clusterModularity()` ，并且最好设置一个参数`as.ratio=TRUE` ，表示两两cluster之间比较，观察到的权重与预期权重之比，这个比值（ratio）越大，表示正在比较的两个cluster之间的细胞间连线数越多，联系越密切。

**下面看一个基于ratio的热图**

```r
# g是利用buildSNNGraph得到的，clust是基于g利用igraph::cluster_walktrap得到的
ratio <- clusterModularity(g, clust, as.ratio=TRUE)
dim(ratio)
## [1] 18 18
# 这是一个矩阵，每一行/列都表示一个cluster，中间的每个ratio值都表示观察到的权重与预期权重之比（可以衡量实际中细胞间连线的数量）。然后我们把ratio画出来看看
library(pheatmap)
pheatmap(log2(ratio+1), cluster_rows=FALSE, cluster_cols=FALSE,
    color=colorRampPalette(c("white", "blue"))(100))
```

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-01-151552.png)

从图中可以看到，深色区域（ratio较大的值）**基本都集中在对角线**，也就是同一个cluter之中。因此，**同一个cluster中会存在最多的细胞间连线数量，而不是与其他的cluster**，这不也正是设置cluster的目的么？一个cluster中的细胞的相关性就是应该高于其他cluster中的细胞

对角线说明了大部分的cluster划分合理，当然**也不排除有一些对角线以外的第二高ratio值**，表示那个cluster与其他的cluster之间也有细胞间连线。

**下面看一个基于ratio的点线图**

就像这种，和之前使用的graph不同，之前的点（node）是细胞，这里的node是cluster，看像不像分化轨迹做的”拟时序分析“？其实就是把细胞聚类后，再聚一次，找到它们的前后相互关系

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-01-152620.png)

这个图是这么做的：

```r
cluster.gr <- igraph::graph_from_adjacency_matrix(log2(ratio+1), 
    mode="upper", weighted=TRUE, diag=FALSE)
set.seed(11001010)
# 权重（weight）越大，线越明显
plot(cluster.gr, edge.width=igraph::E(cluster.gr)$weight*5,
    layout=igraph::layout_with_lgl)
```

## 3 k-means聚类

### **3.1 是什么**

顾名思义，就是通过**不断迭代找出k 个中心(centroid)，然后将各个细胞分配至最近的中心点**。看似晦涩，来看看[k-means 聚类解读](https://blog.csdn.net/huangfei711/article/details/78480078)

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-02-042204.png)

这个方法算法简单、易于实现，具有速度优势。但是存在几个严重的缺点，可能会导致分配的cluster不清不楚：

* k值需要提前定义。假如真实有10种细胞类型，但我现在只定义k=9，那么有两个细胞类型就会硬生生地”合二为一“。但其他方法，例如上面的基于图形的，即使分辨率设置的很低（意思就是看的很模糊），也是能看出这两种不是同一种类型
* 根据它的计算方法，需要不断重复多运行几次，保证得到的cluster结果是稳定的
* 既然k-means是找k个中心，那么有中心就会有区域半径，就会限定数据的范围。但实际上，很多分组并没有常规的数据大小和结构，也就是说，一组数据不会这么”乖乖地“落在我们画好的范围中

看了上面的优缺点后，k-means依然成为了最好用的样本数据压缩方法之一。另外它可以不受细胞密度的影响，保证了即使是数量最多的细胞类型也不会主导下游分析

### **3.2 怎么做？**

基础包中就有函数`kmeans()` ，我们这里基于PCA的结果，指定目标中心数为10，另外也是需要设置随机种子

```r
set.seed(100)
clust.kmeans <- kmeans(reducedDim(sce.pbmc, "PCA"), centers=10)
table(clust.kmeans$cluster)
## 
##   1   2   3   4   5   6   7   8   9  10 
## 472 200  46 515  90 320 241 865 735 438

# 作图
colLabels(sce.pbmc) <- factor(clust.kmeans$cluster)
plotReducedDim(sce.pbmc, "TSNE", colour_by="label")
```

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-01-155839.png)

当然，可以在一开始就设置k大一点（多找一些”核心人物“），避免发生细胞亚群重叠的情况。

```r
set.seed(100)
clust.kmeans2 <- kmeans(reducedDim(sce.pbmc, "PCA"), centers=20)
table(clust.kmeans2$cluster)
## 
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
## 153 172  47 254 125 207 160 334 204 442 163  68 192 271 113 168 124 420  45 260

# 作图 
colLabels(sce.pbmc) <- factor(clust.kmeans2$cluster)
plotTSNE(sce.pbmc, colour_by="label", text_by="label")
```

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-01-160738.png)

### **3.3 分群的评价**

之前了解了k-means的计算方法（计算每个小弟与各个中心大佬的距离，离哪个近就跟哪个），那么就可以用within-cluster sum of squares (**WCSS**) 方法对每个亚群进行评价。顾名思义，这个方法就是寻找最近的距离。

另外，基于WCSS又可以计算root-mean-squared deviation (**RMSD**)，反映的是偏离平均位置的程度。当然，它的计算也很简单，就是用每个亚群的WCSS值除以该亚群的细胞数量然后再开方。

**RMSD值低的好。如果RMSD值低，就表示和其他亚群的混杂程度更小，更像是一个纯粹的亚群**（例如rms值较小的8、10群，它们当中混杂的其他亚群数量很少）。反观rms值大的群，如16、19群，都混杂了其他亚群的细胞（尤其是19群，乍一看上去这个群颜色很单一，但放大看，就会发现其中除了蓝色还有灰色）

```r
ncells <- tabulate(clust.kmeans2$cluster)
tab <- data.frame(wcss=clust.kmeans2$withinss, ncells=ncells)
tab$rms <- sqrt(tab$wcss/tab$ncells)
tab
##         wcss ncells      rms
## 1   2872.081    153 4.332641
## 2   4204.124    172 4.943944
## 3   1443.475     47 5.541862
## 4   4275.279    254 4.102659
## 5   1711.482    125 3.700251
## 6   3054.815    207 3.841557
## 7   1632.526    160 3.194259
## 8   2159.987    334 2.543035
## 9   7026.528    204 5.868881
## 10  2858.314    442 2.542985
## 11  4259.492    163 5.111932
## 12  2288.852     68 5.801688
## 13  2111.912    192 3.316555
## 14  6391.151    271 4.856293
## 15  1603.372    113 3.766847
## 16 12399.822    168 8.591185
## 17  1823.082    124 3.834354
## 18  7026.462    420 4.090192
## 19  2590.561     45 7.587359
## 20  3639.745    260 3.741526
```

再结合上一张t-SNE的图，看到rms指标与t-SNE图上反映的cluster大小基本没有相关性（例如第3群比第12群细胞数量少，同样第12群rms大；第3群比第1群细胞数量少，却是第3群rms大）。这里也侧面反映了之前介绍t-SNE时说的：**t-SNE就是一个可视化的工具，不要试图量化t-SNE结果上的cluster**。**上面的点不能说明什么问题**，t -SNE的算法会让稠密的点变膨胀，让原本稀疏的点变紧凑，**不可以用cluster的大小来衡量亚群的异质性。**&#x53E6;外它并没有义务帮我们保留cluster的相对位置，因此我们**不能根据点的位置远近来确定cluster之间的相似程度，每运行一次点的位置都是变动的。**

> 那么可能会想，如果不通过t-SNE的结果来判断cluster之间的关系，那怎么看呢？

可以通过层次聚类

```r
cent.tree <- hclust(dist(clust.kmeans2$centers), "ward.D2")
plot(cent.tree)
```

同样看到，19群格格不入，与之前t-SNE的可视化结果，以及rms结果都相吻合

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-02-040843.png)

### **3.4 k-means的更高级用法**

之前提到过，k-means可以将邻近的细胞靠一个”中心大佬“来体现，用一个中心(centroid)来表示其余多个点，达到数据压缩的目的。因此，它可以用到其他更复杂的聚类方法中，作为铺垫先压缩一下数据。

例如scran包中有一个函数`clusterSNNGraph()` ，就是整合了k-means和graph-based

* 先利用k-means得到一些关键的中心点（centroid）
* 之后将中心点进行graph-based聚类
* 结果再把每个中心点（centroid）周围的”小弟们“接过来，放在中心点周围

```r
set.seed(0101010)
# 参数kmeans.centers是给第一步k-means用的；参数k是给第二步graph-based聚类用的
kgraph.clusters <- clusterSNNGraph(sce.pbmc, use.dimred="PCA", 
    use.kmeans=TRUE, kmeans.centers=1000, k=5)
table(kgraph.clusters)
## kgraph.clusters
##   1   2   3   4   5   6   7   8   9  10  11 
## 840 137 550 517 220 528 829  46 127  83  45

plotTSNE(sce.pbmc, colour_by=I(kgraph.clusters))
```

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-02-060743.png)

这种方法相比于单纯的graph-based聚类，速度大大提升。（回顾一下之前说的graph-based方法）因为不用去为每个细胞找邻居，从而构建一个图。另外之前说graph-based这种方法使得每个细胞都被强制连接到一定数量的相邻细胞，但k-means的”核心大佬“加入减轻了这个问题。

注意到，上面的代码使用了一个参数`kmeans.centers`， 指的是k-means获得cluster的数量，它的设置取决于是要处理速度还是要数据还原度。设置越大越能表示细胞真实分布，但也为第二步聚类造成了计算量激增。

## 4 层次聚类

### **4.1 是什么**

这是一种历史悠久的聚类方法，最后会给出样本的聚类树（dendrogram）。思路是这样的：**贪婪地将样本聚集成簇，然后将这些簇聚集成更大的簇，依此类推，直到所有样本都属于一个簇为止**。

**层次聚类包含的算法也有很多种，差别在于如何层层聚集**。例如complete linkage方法（`hclust`的默认方法）是保证合并簇之间样本距离最小， Ward’s方法是保证一个簇内部在聚合的同时方差增加的最少，也就是将每次合并时的信息损失最小化。

在实际使用中，层次聚类运行很慢，因此**只适用于小型scRNA数据**。因为需要计算一个细胞间的距离矩阵，如果动辄上千细胞，那么这个计算量会很大。

### **4.2 怎么做**

之前用的PBMC数据集中细胞数太多，于是切换到sce.416b这个数据集，毕竟一二百个细胞还是可以算的

```r
## class: SingleCellExperiment 
## dim: 46604 185 
## metadata(0):
## assays(3): counts logcounts corrected
## rownames(46604): 4933401J01Rik Gm26206 ... CAAA01147332.1
##   CBFB-MYH11-mcherry
## rowData names(4): Length ENSEMBL SYMBOL SEQNAME
## colnames(185): SLX-9555.N701_S502.C89V9ANXX.s_1.r_1
##   SLX-9555.N701_S503.C89V9ANXX.s_1.r_1 ...
##   SLX-11312.N712_S507.H5H5YBBXX.s_8.r_1
##   SLX-11312.N712_S517.H5H5YBBXX.s_8.r_1
## colData names(10): Source Name cell line ... block sizeFactor
## reducedDimNames(2): PCA TSNE
## altExpNames(2): ERCC SIRV
```

首先利用前几个PCs计算细胞间的距离

```r
dist.416b <- dist(reducedDim(sce.416b, "PCA"))
```

然后使用Ward‘s方法进行聚类

```r
tree.416b <- hclust(dist.416b, "ward.D2")
```

准备画复杂一点的聚类树

```r
library(dendextend)
# seq_along根据tree.416b$labels，也就是185个细胞的ID，产生了1-185个数字，结果就是把字符串变数字
tree.416b$labels <- seq_along(tree.416b$labels)
dend <- as.dendrogram(tree.416b, hang=0.1)
# 合并两种批次信息
combined.fac <- paste0(sce.416b$block, ".", 
    sub(" .*", "", sce.416b$phenotype))
> table(combined.fac)
combined.fac
20160113.induced    20160113.wild 20160325.induced    20160325.wild 
              46               47               47               45

# 首先按照不同处理设置蓝色和红色，然后按照批次设置淡蓝色和淡红色
labels_colors(dend) <- c(
    `20160113.wild`="blue",
    `20160113.induced`="red",
    `20160325.wild`="dodgerblue",
    `20160325.induced`="salmon"
)[combined.fac][order.dendrogram(dend)]

plot(dend)
```

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-02-063949.png)

图中可以看到，处理前后确实形成了红、蓝分界

另外更推荐使用Ward’s方法，它在合并时受到的cluster差异影响最小（毕竟人家的出发点就是每次合并时的信息损失最小化，所以如果差异很大，那我就不合并呗）。

**为了得到更清楚地分群结果，可以对树进行”修建“**

* 最简单的是`cutree()` ，例如`cutree(hc, 4)` 手动切分数量
* 复杂一点但更好用的是`cutreeDynamic()`，来自`dynamicTreeCut` 这个包

```r
library(dynamicTreeCut)
clust.416b <- cutreeDynamic(tree.416b, distM=as.matrix(dist.416b),
    minClusterSize=10, deepSplit=1)
table(clust.416b)
## clust.416b
##  1  2  3  4 
## 78 69 24 14

labels_colors(dend) <- clust.416b[order.dendrogram(dend)]
plot(dend)
```

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-02-070851.png)

一般这个结果与t-SNE的结果是吻合的

```r
colLabels(sce.416b) <- factor(clust.416b)
plotReducedDim(sce.416b, "TSNE", colour_by="label")
```

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-02-071938.png)

不过这里注意到，t-SNE中cluster2拆分开了，值得怀疑：是真的聚类出了问题，还是t-SNE本身的显示问题，毕竟我们知道t-SNE有时真的会把一个亚群拆开画。于是继续下面👇的探索评价

### **4.3 评价**

可以借助”轮廓图“（silhouette width）检查分群的质量。会对每个细胞都计算一个silhouette width值，如果一个细胞的**width值为正并且越大**，表示相对于其他亚群的细胞，这个细胞和它所在亚群中的细胞更接近，**分群效果越好**；如果width为负，就表示这个亚群的这个细胞和其他亚群的细胞更接近，即分群效果不太理想。

```r
library(cluster)
sil <- silhouette(clust.416b, dist = dist.416b)
> colnames(sil)
[1] "cluster"   "neighbor"  "sil_width"

# 设置每个cluster的颜色
clust.col <- scater:::.get_palette("tableau10medium") # hidden scater colours
# 这一句的意思是：如果sil_width>0，就属于和自己接近，否则属于和其他亚群接近
sil.cols <- clust.col[ifelse(sil[,3] > 0, sil[,1], sil[,2])]
sil.cols <- sil.cols[order(-sil[,1], sil[,3])]

plot(sil, main = paste(length(unique(clust.416b)), "clusters"), 
     border=sil.cols, col=sil.cols, do.col.sort=FALSE)
```

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-02-092004.png)

图中可以反映的问题：

* 这个图是对上面得到的各个cluster中的细胞做的barplot，每个cluster是一种颜色。
* 这里看到cluster2都是正值并且横坐标width数值还是最大的，因此说明了分群没有问题，之前t-SNE的结果的确是它本身的问题
* 如果发现width值较小，表示分群结果一般，还有可能是分群多度，本来属于一个群的，又被拆分成小群
* 利用这个图，我们就能调整之前的参数，来调整分群效果，不过也不需要太过纠结完美的分群结果。因为即使图上看似合理的分群，可能实际上也不会得到更多的生物信息

## 5 评价分群的稳定性

### **什么叫做分群的稳定性？**

就是分群之前操作的小小波动，也不会影响到最后的分群结果；高稳定性的分群结果其实也方便了别人进行重复。scran利用自助法（bootstrapping）进行评价

> 在统计学中，自助法（Bootstrap Method，Bootstrapping，或自助抽样法）是一种从**给定训练集中有放回的均匀抽样**，也就是说，**每当选中一个样本，它等可能地被再次选中并被再次添加到训练集中**。 自助法由Bradley Efron于1979年在《Annals of Statistics》上发表。

细胞有放回的抽样，得到一个”自助重复“数据集，然后用他进行聚类，看看是否能得到相同的亚群

```r
# 例如使用graph-based聚类
myClusterFUN <- function(x) {
    g <- buildSNNGraph(x, use.dimred="PCA", type="jaccard")
    igraph::cluster_louvain(g)$membership
}

originals <- myClusterFUN(sce.pbmc)

# 进行自助法判断
set.seed(0010010100)
coassign <- bootstrapCluster(sce.pbmc, FUN=myClusterFUN, clusters=originals)
dim(coassign)
## [1] 19 19
```

结果返回一个coassignment的矩阵，表示随机从clusterX与clusterY中挑一个细胞，这两个细胞在经历”自助法“判断后，依然落在同一个cluster的概率。我们希望落在对角线上的概率更高（因为对角线表示自己和自己比较，说明多次重复后细胞还是属于这个cluster）

```r
# 借助热图判断
pheatmap(coassign, cluster_row=FALSE, cluster_col=FALSE,
    color=rev(viridis::magma(100)))
```

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-02-082643.png)

**需要注意的是**

* 这个自助法（Bootstrapping）是一个通用的评价cluster稳定性的方法，适用于各种聚类算法
* 高稳定性的亚群不一定是高质量的，因为如果一个亚群质量一直很差，它也很稳定

## 6 数据集取小再聚类 | Subclustering

首先通过常规操作得到分群结果，找到某些感兴趣的亚群，然后在某一个cluster中再一次进行特征基因挑选、聚类操作，从而增加分析的精确度。

举个例子：

首先利用整个PBMC数据集，根据T细胞的marker基因筛出可能是T细胞的亚群，然后推测其中某个亚群可能是记忆T细胞

```r
# 看到buildSNNGraph就知道是graph-based 画图，然后用igraph去鉴定
g.full <- buildSNNGraph(sce.pbmc, use.dimred = 'PCA')
clust.full <- igraph::cluster_walktrap(g.full)$membership

# 画几个marker基因的表达量，使用的是log-normalized表达量
plotExpression(sce.pbmc, features=c("CD3E", "CCR7", "CD69", "CD44"),
    x=I(factor(clust.full)), colour_by=I(factor(clust.full)))
```

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-02-084532.png)

**现在推测cluster6可能是记忆T细胞**

得到亚群cluster6的CD4+和CD8+ 的亚亚群

```r
memory <- 6 
sce.memory <- sce.pbmc[,clust.full==memory]
# 对它重新挑一遍特征基因，进行PCA
dec.memory <- modelGeneVar(sce.memory)
sce.memory <- denoisePCA(sce.memory, technical=dec.memory,
    subset.row=getTopHVGs(dec.memory, prop=0.1))

# 再走一遍聚类
g.memory <- buildSNNGraph(sce.memory, use.dimred="PCA")
clust.memory <- igraph::cluster_walktrap(g.memory)$membership

# 最后画图（最后就聚了2类，变成因子型放在x轴上），并且用CD8A和CD4进行验证
plotExpression(sce.memory, features=c("CD8A", "CD4"),
    x=I(factor(clust.memory)))
```

![](https://jieandze1314-1255603621.cos.ap-guangzhou.myqcloud.com/blog/2020-07-02-093404.png)

这种数据集取小再聚类的方法，可以**简化亚群的解读**。只需要在某一个亚群的背景下继续挖掘，就可能会获得重要信息，一层层递进。比如上面的思路：一堆细胞=》分18群=》鉴定T细胞的亚群=》取出某个T细胞亚群再去鉴定记忆T细胞