单细胞交响乐7-聚类分群
原书名为：Orchestrating Single-Cell Analysis
我给它取名叫“单细胞交响乐”
因为单细胞分析就像一个大乐团，需要各个流程的协同配合
1 前言
聚类是一种无监督的学习过程，在scRNA分析中，根据表达谱相似性对不同细胞进行分组。
在后面根据marker基因对各个cluster进行注释后，cluster就可以代表不同的细胞类型或状态。因此，聚类是将枯燥的数据与感兴趣生物学意义联系起来的一个桥梁。
不过这里需要明确一下，聚类分的群和细胞类型还不是一回事
分群是根据经验计算的结果
细胞类型有真正的生物学意义
至于这种问题：“我应该设置多少cluster数量？”，则是个伪命题。因为我们可以根据需要，使用不同方法，得到不同数量的clusters，每个cluster都是高维空间数据的一部分。
不过倒是可以问：“我们的clusters与细胞类型之间的拟合程度如何？”。不过这个问题很难回答，因为取决于自己的分析目标。有的人对于能得到主要的细胞类型就很满意，而有的人还想再进一步，得到细胞亚群；还有的，更想基于细胞亚群，达到不同细胞状态的分辨率（如代谢活性、胁迫反应等）。另外，两种截然不同的分群结果，可能都具有生物学意义，只是关注点不同而已。
对于分群聚类的操作，更像是显微镜的操作，可以通过调整参数（使用不同的镜头组合）来获得不同的分辨率；另外还可以探索多种聚类方法，看到数据的不同面。
数据准备之PBMC
具体的操作在上一章也提到过
# 下载
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)
首先使用scran包的buildSNNGraph()函数，并且使用PCA的前几个PCs；之后利用 igraph包中的Walktrap方法去鉴定clusters
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对象中，保存成一个因子
library(scater)
colLabels(sce.pbmc) <- factor(clust)
plotReducedDim(sce.pbmc, "TSNE", colour_by="label")
从上面的计算方法可以看出，有一个很重要的参数是k ，含义是：the number of nearest neighbors used to construct the graph。如果k设置越大，得到的图之间联通程度越高，cluster也越大。因此这个参数也是可以不断尝试的
# 比如设置较大的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的降维结果都是紧密相关的
set.seed(2000)
reducedDim(sce.pbmc, "force") <- igraph::layout_with_fr(g)
plotReducedDim(sce.pbmc, colour_by="label", dimred="force")
2.3 其他的参数
关于如何计算权重，在buildSNNGraph()函数中可以设置：
type="number"：根据近邻数量计算type="jaccard"：根据Jaccard index of the two sets of neighbors计算
g.num <- buildSNNGraph(sce.pbmc, use.dimred="PCA", type="number")
g.jaccard <- buildSNNGraph(sce.pbmc, use.dimred="PCA", type="jaccard")
上面得到的各种结果，都是来自igraph包的graph对象，因此可以用于igraph包的各种community检测方法：
# 之前使用的是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
最后就可以比较两种不同的聚类方法差异
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]])
可以看出，Infomap比Walktrap得到了更多的cluster（比较精细），而fast-greedy结果相对较少（比较粗糙）
可以通过cut_at()指定分群的数量
# 指定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的热图
# 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))
从图中可以看到，深色区域（ratio较大的值）基本都集中在对角线，也就是同一个cluter之中。因此，同一个cluster中会存在最多的细胞间连线数量，而不是与其他的cluster，这不也正是设置cluster的目的么？一个cluster中的细胞的相关性就是应该高于其他cluster中的细胞
对角线说明了大部分的cluster划分合理，当然也不排除有一些对角线以外的第二高ratio值，表示那个cluster与其他的cluster之间也有细胞间连线。
下面看一个基于ratio的点线图
就像这种，和之前使用的graph不同，之前的点（node）是细胞，这里的node是cluster，看像不像分化轨迹做的”拟时序分析“？其实就是把细胞聚类后，再聚一次，找到它们的前后相互关系
这个图是这么做的：
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)，然后将各个细胞分配至最近的中心点。看似晦涩，来看看https://blog.csdn.net/huangfei711/article/details/78480078这里的解读：
这个方法算法简单、易于实现，具有速度优势。但是存在几个严重的缺点，可能会导致分配的cluster不清不楚：
k值需要提前定义。假如真实有10种细胞类型，但我现在只定义k=9，那么有两个细胞类型就会硬生生地”合二为一“。但其他方法，例如上面的基于图形的，即使分辨率设置的很低（意思就是看的很模糊），也是能看出这两种不是同一种类型
根据它的计算方法，需要不断重复多运行几次，保证得到的cluster结果是稳定的
既然k-means是找k个中心，那么有中心就会有区域半径，就会限定数据的范围。但实际上，很多分组并没有常规的数据大小和结构，也就是说，一组数据不会这么”乖乖地“落在我们画好的范围中
看了上面的优缺点后，k-means依然成为了最好用的样本数据压缩方法之一。另外它可以不受细胞密度的影响，保证了即使是数量最多的细胞类型也不会主导下游分析
3.2 怎么做？
基础包中就有函数kmeans() ，我们这里基于PCA的结果，指定目标中心数为10，另外也是需要设置随机种子
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")
当然，可以在一开始就设置k大一点（多找一些”核心人物“），避免发生细胞亚群重叠的情况。
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")
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群，乍一看上去这个群颜色很单一，但放大看，就会发现其中除了蓝色还有灰色）
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的大小来衡量亚群的异质性。另外它并没有义务帮我们保留cluster的相对位置，因此我们不能根据点的位置远近来确定cluster之间的相似程度，每运行一次点的位置都是变动的。
那么可能会想，如果不通过t-SNE的结果来判断cluster之间的关系，那怎么看呢？
可以通过层次聚类
cent.tree <- hclust(dist(clust.kmeans2$centers), "ward.D2")
plot(cent.tree)
同样看到，19群格格不入，与之前t-SNE的可视化结果，以及rms结果都相吻合
3.4 k-means的更高级用法
之前提到过，k-means可以将邻近的细胞靠一个”中心大佬“来体现，用一个中心(centroid)来表示其余多个点，达到数据压缩的目的。因此，它可以用到其他更复杂的聚类方法中，作为铺垫先压缩一下数据。
例如scran包中有一个函数clusterSNNGraph() ，就是整合了k-means和graph-based
先利用k-means得到一些关键的中心点（centroid）
之后将中心点进行graph-based聚类
结果再把每个中心点（centroid）周围的”小弟们“接过来，放在中心点周围
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))
这种方法相比于单纯的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这个数据集，毕竟一二百个细胞还是可以算的
## 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计算细胞间的距离
dist.416b <- dist(reducedDim(sce.416b, "PCA"))
然后使用Ward‘s方法进行聚类
tree.416b <- hclust(dist.416b, "ward.D2")
准备画复杂一点的聚类树
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)
图中可以看到，处理前后确实形成了红、蓝分界
另外更推荐使用Ward’s方法，它在合并时受到的cluster差异影响最小（毕竟人家的出发点就是每次合并时的信息损失最小化，所以如果差异很大，那我就不合并呗）。
为了得到更清楚地分群结果，可以对树进行”修建“
最简单的是
cutree()，例如cutree(hc, 4)手动切分数量复杂一点但更好用的是
cutreeDynamic()，来自dynamicTreeCut这个包
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)
一般这个结果与t-SNE的结果是吻合的
colLabels(sce.416b) <- factor(clust.416b)
plotReducedDim(sce.416b, "TSNE", colour_by="label")
不过这里注意到，t-SNE中cluster2拆分开了，值得怀疑：是真的聚类出了问题，还是t-SNE本身的显示问题，毕竟我们知道t-SNE有时真的会把一个亚群拆开画。于是继续下面👇的探索评价
4.3 评价
可以借助”轮廓图“（silhouette width）检查分群的质量。会对每个细胞都计算一个silhouette width值，如果一个细胞的width值为正并且越大，表示相对于其他亚群的细胞，这个细胞和它所在亚群中的细胞更接近，分群效果越好；如果width为负，就表示这个亚群的这个细胞和其他亚群的细胞更接近，即分群效果不太理想。
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)
图中可以反映的问题：
这个图是对上面得到的各个cluster中的细胞做的barplot，每个cluster是一种颜色。
这里看到cluster2都是正值并且横坐标width数值还是最大的，因此说明了分群没有问题，之前t-SNE的结果的确是它本身的问题
如果发现width值较小，表示分群结果一般，还有可能是分群多度，本来属于一个群的，又被拆分成小群
利用这个图，我们就能调整之前的参数，来调整分群效果，不过也不需要太过纠结完美的分群结果。因为即使图上看似合理的分群，可能实际上也不会得到更多的生物信息
5 评价分群的稳定性
什么叫做分群的稳定性？
就是分群之前操作的小小波动，也不会影响到最后的分群结果；高稳定性的分群结果其实也方便了别人进行重复。scran利用自助法（bootstrapping）进行评价
在统计学中，自助法（Bootstrap Method，Bootstrapping，或自助抽樣法）是一种从给定训练集中有放回的均匀抽样，也就是说，每当选中一个样本，它等可能地被再次选中并被再次添加到训练集中。自助法由Bradley Efron于1979年在《Annals of Statistics》上發表。
细胞有放回的抽样，得到一个”自助重复“数据集，然后用他进行聚类，看看是否能得到相同的亚群
# 例如使用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）
# 借助热图判断
pheatmap(coassign, cluster_row=FALSE, cluster_col=FALSE,
color=rev(viridis::magma(100)))
需要注意的是
这个自助法（Bootstrapping）是一个通用的评价cluster稳定性的方法，适用于各种聚类算法
高稳定性的亚群不一定是高质量的，因为如果一个亚群质量一直很差，它也很稳定
6 数据集取小再聚类 | Subclustering
首先通过常规操作得到分群结果，找到某些感兴趣的亚群，然后在某一个cluster中再一次进行特征基因挑选、聚类操作，从而增加分析的精确度。
举个例子：
首先利用整个PBMC数据集，根据T细胞的marker基因筛出可能是T细胞的亚群，然后推测其中某个亚群可能是记忆T细胞
# 看到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)))
现在推测cluster6可能是记忆T细胞
得到亚群cluster6的CD4+和CD8+ 的亚亚群
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)))
这种数据集取小再聚类的方法，可以简化亚群的解读。只需要在某一个亚群的背景下继续挖掘，就可能会获得重要信息，一层层递进。比如上面的思路：一堆细胞=》分18群=》鉴定T细胞的亚群=》取出某个T细胞亚群再去鉴定记忆T细胞
点击底部的“阅读原文”，获得更好的阅读体验哦😻
如果你对单细胞转录组研究感兴趣，但又不知道如何入门，也许你可以关注一下下面的课程
生信爆款入门-第7期（线上直播4周，马拉松式陪伴，带你入门）你的生物信息入门课
数据挖掘学习班第5期（线上直播3周，马拉松式陪伴，带你入门） 医学生/医生首选技能提高课
生信技能树的2019年终总结 你的生物信息成长宝藏
看完记得顺手点个“在看”哦！
长按扫码可关注