# 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细胞


---

# Agent Instructions: Querying This Documentation

If you need additional information that is not directly available in this page, you can query the documentation dynamically by asking a question.

Perform an HTTP GET request on the current page URL with the `ask` query parameter:

```
GET https://jieandze1314.osca.top/03/03-5.md?ask=<question>
```

The question should be specific, self-contained, and written in natural language.
The response will contain a direct answer to the question and relevant excerpts and sources from the documentation.

Use this mechanism when the answer is not explicitly present in the current page, you need clarification or additional context, or you want to retrieve related documentation sections.