3.5 聚类

刘小泽写于2020.6.30 + 7.2

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),然后将各个细胞分配至最近的中心点。看似晦涩,来看看k-means 聚类解读

这个方法算法简单、易于实现,具有速度优势。但是存在几个严重的缺点,可能会导致分配的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细胞