劉小澤寫于2020.6.30 + 7.2
為何取名叫“交響樂”栋猖？因?yàn)閱渭?xì)胞分析就像一個(gè)大樂團(tuán)，需要各個(gè)流程的協(xié)同配合
單細(xì)胞交響樂1-常用的數(shù)據(jù)結(jié)構(gòu)SingleCellExperiment
單細(xì)胞交響樂2-scRNAseq從實(shí)驗(yàn)到下游簡(jiǎn)介
單細(xì)胞交響樂3-細(xì)胞質(zhì)控
單細(xì)胞交響樂4-歸一化
單細(xì)胞交響樂5-挑選高變化基因
單細(xì)胞交響樂6-降維

1 前言

聚類是一種無監(jiān)督的學(xué)習(xí)過程烹骨，在scRNA分析中，根據(jù)表達(dá)譜相似性對(duì)不同細(xì)胞進(jìn)行分組驹马。

在后面根據(jù)marker基因?qū)Ω鱾€(gè)cluster進(jìn)行注釋后卷拘，cluster就可以代表不同的細(xì)胞類型或狀態(tài)仰禽。因此氮墨，聚類是將枯燥的數(shù)據(jù)與感興趣生物學(xué)意義聯(lián)系起來的一個(gè)橋梁。

不過這里需要明確一下坟瓢，聚類分的群和細(xì)胞類型還不是一回事

• 分群是根據(jù)經(jīng)驗(yàn)計(jì)算的結(jié)果
• 細(xì)胞類型有真正的生物學(xué)意義

至于這種問題：“我應(yīng)該設(shè)置多少cluster數(shù)量？”犹撒，則是個(gè)偽命題折联。因?yàn)槲覀兛梢愿鶕?jù)需要，使用不同方法识颊，得到不同數(shù)量的clusters诚镰，每個(gè)cluster都是高維空間數(shù)據(jù)的一部分。

不過倒是可以問：“我們的clusters與細(xì)胞類型之間的擬合程度如何祥款？”清笨。不過這個(gè)問題很難回答，因?yàn)槿Q于自己的分析目標(biāo)刃跛。有的人對(duì)于能得到主要的細(xì)胞類型就很滿意抠艾，而有的人還想再進(jìn)一步，得到細(xì)胞亞群桨昙；還有的检号，更想基于細(xì)胞亞群腌歉，達(dá)到不同細(xì)胞狀態(tài)的分辨率（如代謝活性、脅迫反應(yīng)等）齐苛。另外翘盖，兩種截然不同的分群結(jié)果，可能都具有生物學(xué)意義凹蜂，只是關(guān)注點(diǎn)不同而已馍驯。

對(duì)于分群聚類的操作，更像是顯微鏡的操作玛痊，可以通過調(diào)整參數(shù)（使用不同的鏡頭組合）來獲得不同的分辨率汰瘫；另外還可以探索多種聚類方法，看到數(shù)據(jù)的不同面卿啡。

數(shù)據(jù)準(zhǔn)備之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")

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

# 質(zhì)控
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)

# 利用數(shù)據(jù)分布來表示技術(shù)噪音吟吝，并挑選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")

# 看下結(jié)果
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):

補(bǔ)充知識(shí)點(diǎn)之聚類與分類

• 分類：類別是已知的，通過對(duì)已知分類的數(shù)據(jù)進(jìn)行訓(xùn)練和學(xué)習(xí)颈娜，找到這些不同類的特征剑逃，再對(duì)未分類的數(shù)據(jù)進(jìn)行分類。屬于有監(jiān)督學(xué)習(xí)官辽。

• 聚類：事先不知道數(shù)據(jù)會(huì)分為幾類蛹磺，通過聚類分析將數(shù)據(jù)聚合成幾個(gè)群體。聚類不需要對(duì)數(shù)據(jù)進(jìn)行訓(xùn)練和學(xué)習(xí)同仆。屬于無監(jiān)督學(xué)習(xí)萤捆。

2 基于圖形的聚類 | Graph-based clustering

2.1 是什么？

這種聚類方法在Seurat中使用最流行俗批。最基礎(chǔ)的想法是：我們首先構(gòu)建一個(gè)圖俗或，其中每個(gè)節(jié)點(diǎn)都是一個(gè)細(xì)胞，它與高維空間中最鄰近的細(xì)胞相連岁忘。連線基于細(xì)胞之間的相似性計(jì)算權(quán)重辛慰，權(quán)重越高，表示細(xì)胞間關(guān)系更密切干像。

如果一群細(xì)胞之間的權(quán)重高于另一群細(xì)胞帅腌，那么這一群細(xì)胞就被當(dāng)做一個(gè)群體 “communiity”

這個(gè)方法最大的優(yōu)點(diǎn)就是：可縮放性。只需要根據(jù) k-nearest neighbor（KNN）進(jìn)行搜索麻汰，自己去尋找去拓展速客。相比于k-means、 Gaussian mixture models等方法五鲫，不需要預(yù)先對(duì)cluster形狀以及cluster中細(xì)胞分布做一個(gè)估計(jì)溺职。

這種方法使得每個(gè)細(xì)胞都被強(qiáng)制連接到一定數(shù)量的相鄰細(xì)胞，這減少了僅由一兩個(gè)離群細(xì)胞組成的無意義群體的風(fēng)險(xiǎn)

不過這個(gè)方法也有缺點(diǎn)：

缺點(diǎn)一：

它最后只保留了鄰近細(xì)胞之間的聯(lián)系，除此以外沒有保存辅愿。

缺點(diǎn)二：

另外正是由于它主要關(guān)注相鄰的權(quán)重智亮，導(dǎo)致一些噪音也可能聚集在一起，形成一個(gè)community点待，從而影響聚類的分辨率阔蛉，可能導(dǎo)致數(shù)據(jù)異質(zhì)性的過度解讀

2.2 怎么做？

在使用時(shí)癞埠，需要考慮幾個(gè)問題：

• 構(gòu)建這個(gè)圖需要設(shè)置幾個(gè)鄰居状原？
• 用什么方法確定邊的權(quán)重？
• 用什么檢測(cè)community的方法來定義cluster苗踪？

下面的例子中颠区，就定義了使用一個(gè)細(xì)胞與10個(gè)鄰居的關(guān)系，進(jìn)而構(gòu)建shared nearest neighbor graph（SNNG）通铲，如果兩個(gè)細(xì)胞的鄰居中有一個(gè)是共享的毕莱，則通過一條邊連接，之后再根據(jù)highest average rank of the shared neighbors定義邊的權(quán)重(Xu and Su 2015)

首先使用scran包的buildSNNGraph()函數(shù)颅夺，并且使用PCA的前幾個(gè)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對(duì)象中，保存成一個(gè)因子

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

從上面的計(jì)算方法可以看出吧黄，有一個(gè)很重要的參數(shù)是k 部服，含義是：the number of nearest neighbors used to construct the graph。如果k設(shè)置越大拗慨，得到的圖之間聯(lián)通程度越高廓八，cluster也越大。因此這個(gè)參數(shù)也是可以不斷嘗試的

# 比如設(shè)置較大的k赵抢，就會(huì)得到比k=10更少的cluster剧蹂，但同時(shí)每個(gè)cluster平均的細(xì)胞數(shù)更多
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

# 如果設(shè)置更小的k，會(huì)得到數(shù)量更多的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

這個(gè)結(jié)果可以利用force-directed的布局出圖烦却，它與t-SNE宠叼、UMAP的降維結(jié)果都是緊密相關(guān)的

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

2.3 其他的參數(shù)

關(guān)于如何計(jì)算權(quán)重，在buildSNNGraph()函數(shù)中可以設(shè)置：

• type="number" ：根據(jù)近鄰數(shù)量計(jì)算
• type="jaccard" ：根據(jù)Jaccard index of the two sets of neighbors計(jì)算

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

上面得到的各種結(jié)果短绸，都是來自igraph包的graph對(duì)象车吹，因此可以用于igraph包的各種community檢測(cè)方法：

# 之前使用的是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（比較精細(xì)），而fast-greedy結(jié)果相對(duì)較少（比較粗糙）

可以通過cut_at()指定分群的數(shù)量

# 指定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默認(rèn)使用rank-based的權(quán)重計(jì)算方法朝卒，Seurat使用Jaccard-based的方法（后面輔助Louvain聚類）

2.4 分群的評(píng)價(jià)

當(dāng)使用graph-based方法時(shí)证逻，模塊化（modularity）是一個(gè)評(píng)價(jià)cluster的指標(biāo)。modularity值越大抗斤，表示cluster內(nèi)部的細(xì)胞與細(xì)胞之間連線（edge）數(shù)最多囚企，內(nèi)部分散效果越好丈咐，從而避免了不同cluster之間鄰近的細(xì)胞形成連線。

之前通過cluster_walktrap的方法得到了18個(gè)cluster龙宏，下面就看看這個(gè)方法做的如何：會(huì)用到一個(gè)函數(shù)clusterModularity() 棵逊，并且最好設(shè)置一個(gè)參數(shù)as.ratio=TRUE ，表示兩兩cluster之間比較银酗，觀察到的權(quán)重與預(yù)期權(quán)重之比辆影，這個(gè)比值（ratio）越大，表示正在比較的兩個(gè)cluster之間的細(xì)胞間連線數(shù)越多黍特，聯(lián)系越密切蛙讥。

下面看一個(gè)基于ratio的熱圖

# g是利用buildSNNGraph得到的，clust是基于g利用igraph::cluster_walktrap得到的
ratio <- clusterModularity(g, clust, as.ratio=TRUE)
dim(ratio)
## [1] 18 18
# 這是一個(gè)矩陣灭衷，每一行/列都表示一個(gè)cluster次慢，中間的每個(gè)ratio值都表示觀察到的權(quán)重與預(yù)期權(quán)重之比（可以衡量實(shí)際中細(xì)胞間連線的數(shù)量）。然后我們把ratio畫出來看看
library(pheatmap)
pheatmap(log2(ratio+1), cluster_rows=FALSE, cluster_cols=FALSE,
    color=colorRampPalette(c("white", "blue"))(100))

從圖中可以看到翔曲，深色區(qū)域（ratio較大的值）基本都集中在對(duì)角線迫像，也就是同一個(gè)cluter之中。因此部默，同一個(gè)cluster中會(huì)存在最多的細(xì)胞間連線數(shù)量侵蒙，而不是與其他的cluster，這不也正是設(shè)置cluster的目的么傅蹂？一個(gè)cluster中的細(xì)胞的相關(guān)性就是應(yīng)該高于其他cluster中的細(xì)胞

對(duì)角線說明了大部分的cluster劃分合理纷闺，當(dāng)然也不排除有一些對(duì)角線以外的第二高ratio值，表示那個(gè)cluster與其他的cluster之間也有細(xì)胞間連線份蝴。

下面看一個(gè)基于ratio的點(diǎn)線圖

就像這種犁功，和之前使用的graph不同，之前的點(diǎn)（node）是細(xì)胞婚夫，這里的node是cluster浸卦，看像不像分化軌跡做的”擬時(shí)序分析“？其實(shí)就是把細(xì)胞聚類后案糙，再聚一次限嫌，找到它們的前后相互關(guān)系

這個(gè)圖是這么做的：

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

3 k-means聚類

3.1 是什么

顧名思義时捌，就是通過不斷迭代找出k 個(gè)中心(centroid)怒医，然后將各個(gè)細(xì)胞分配至最近的中心點(diǎn)∩萏郑看似晦澀稚叹，來看看https://blog.csdn.net/huangfei711/article/details/78480078這里的解讀：

這個(gè)方法算法簡(jiǎn)單、易于實(shí)現(xiàn)，具有速度優(yōu)勢(shì)扒袖。但是存在幾個(gè)嚴(yán)重的缺點(diǎn)塞茅，可能會(huì)導(dǎo)致分配的cluster不清不楚：

• k值需要提前定義。假如真實(shí)有10種細(xì)胞類型季率，但我現(xiàn)在只定義k=9野瘦，那么有兩個(gè)細(xì)胞類型就會(huì)硬生生地”合二為一“。但其他方法飒泻，例如上面的基于圖形的缅刽，即使分辨率設(shè)置的很低（意思就是看的很模糊），也是能看出這兩種不是同一種類型
• 根據(jù)它的計(jì)算方法蠢络，需要不斷重復(fù)多運(yùn)行幾次衰猛，保證得到的cluster結(jié)果是穩(wěn)定的
• 既然k-means是找k個(gè)中心，那么有中心就會(huì)有區(qū)域半徑刹孔，就會(huì)限定數(shù)據(jù)的范圍啡省。但實(shí)際上，很多分組并沒有常規(guī)的數(shù)據(jù)大小和結(jié)構(gòu)髓霞，也就是說卦睹，一組數(shù)據(jù)不會(huì)這么”乖乖地“落在我們畫好的范圍中

看了上面的優(yōu)缺點(diǎn)后，k-means依然成為了最好用的樣本數(shù)據(jù)壓縮方法之一方库。另外它可以不受細(xì)胞密度的影響结序，保證了即使是數(shù)量最多的細(xì)胞類型也不會(huì)主導(dǎo)下游分析

3.2 怎么做？

基礎(chǔ)包中就有函數(shù)kmeans() 纵潦，我們這里基于PCA的結(jié)果徐鹤，指定目標(biāo)中心數(shù)為10，另外也是需要設(shè)置隨機(jī)種子

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")

當(dāng)然邀层，可以在一開始就設(shè)置k大一點(diǎn)（多找一些”核心人物“）返敬，避免發(fā)生細(xì)胞亞群重疊的情況。

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 分群的評(píng)價(jià)

之前了解了k-means的計(jì)算方法（計(jì)算每個(gè)小弟與各個(gè)中心大佬的距離寥院，離哪個(gè)近就跟哪個(gè)）劲赠，那么就可以用within-cluster sum of squares (WCSS) 方法對(duì)每個(gè)亞群進(jìn)行評(píng)價(jià)。顧名思義秸谢，這個(gè)方法就是尋找最近的距離凛澎。

另外，基于WCSS又可以計(jì)算root-mean-squared deviation (RMSD)估蹄，反映的是偏離平均位置的程度塑煎。當(dāng)然，它的計(jì)算也很簡(jiǎn)單元媚，就是用每個(gè)亞群的WCSS值除以該亞群的細(xì)胞數(shù)量然后再開方轧叽。

RMSD值低的好贼陶。如果RMSD值低右锨，就表示和其他亞群的混雜程度更小岂丘，更像是一個(gè)純粹的亞群（例如rms值較小的8般妙、10群梳玫，它們當(dāng)中混雜的其他亞群數(shù)量很少）导梆。反觀rms值大的群贩绕，如16艳汽、19群嗤无，都混雜了其他亞群的細(xì)胞（尤其是19群震束，乍一看上去這個(gè)群顏色很單一，但放大看当犯，就會(huì)發(fā)現(xiàn)其中除了藍(lán)色還有灰色）

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

再結(jié)合上一張t-SNE的圖垢村，看到rms指標(biāo)與t-SNE圖上反映的cluster大小基本沒有相關(guān)性（例如第3群比第12群細(xì)胞數(shù)量少，同樣第12群rms大嚎卫；第3群比第1群細(xì)胞數(shù)量少嘉栓，卻是第3群rms大）。這里也側(cè)面反映了之前介紹t-SNE時(shí)說的：t-SNE就是一個(gè)可視化的工具拓诸，不要試圖量化t-SNE結(jié)果上的cluster侵佃。上面的點(diǎn)不能說明什么問題，t -SNE的算法會(huì)讓稠密的點(diǎn)變膨脹奠支，讓原本稀疏的點(diǎn)變緊湊馋辈，不可以用cluster的大小來衡量亞群的異質(zhì)性。另外它并沒有義務(wù)幫我們保留cluster的相對(duì)位置倍谜，因此我們不能根據(jù)點(diǎn)的位置遠(yuǎn)近來確定cluster之間的相似程度迈螟，每運(yùn)行一次點(diǎn)的位置都是變動(dòng)的。

那么可能會(huì)想尔崔，如果不通過t-SNE的結(jié)果來判斷cluster之間的關(guān)系井联，那怎么看呢？

可以通過層次聚類

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

同樣看到您旁，19群格格不入烙常，與之前t-SNE的可視化結(jié)果，以及rms結(jié)果都相吻合

3.4 k-means的更高級(jí)用法

之前提到過鹤盒，k-means可以將鄰近的細(xì)胞靠一個(gè)”中心大佬“來體現(xiàn)蚕脏，用一個(gè)中心(centroid)來表示其余多個(gè)點(diǎn)，達(dá)到數(shù)據(jù)壓縮的目的侦锯。因此驼鞭，它可以用到其他更復(fù)雜的聚類方法中，作為鋪墊先壓縮一下數(shù)據(jù)尺碰。

例如scran包中有一個(gè)函數(shù)clusterSNNGraph() 挣棕，就是整合了k-means和graph-based

• 先利用k-means得到一些關(guān)鍵的中心點(diǎn)（centroid）
• 之后將中心點(diǎn)進(jìn)行g(shù)raph-based聚類
• 結(jié)果再把每個(gè)中心點(diǎn)（centroid）周圍的”小弟們“接過來译隘，放在中心點(diǎn)周圍

set.seed(0101010)
# 參數(shù)kmeans.centers是給第一步k-means用的；參數(shù)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方法）因?yàn)椴挥萌槊總€(gè)細(xì)胞找鄰居，從而構(gòu)建一個(gè)圖词身。另外之前說graph-based這種方法使得每個(gè)細(xì)胞都被強(qiáng)制連接到一定數(shù)量的相鄰細(xì)胞厅目，但k-means的”核心大佬“加入減輕了這個(gè)問題。

注意到法严，上面的代碼使用了一個(gè)參數(shù)kmeans.centers损敷， 指的是k-means獲得cluster的數(shù)量，它的設(shè)置取決于是要處理速度還是要數(shù)據(jù)還原度深啤。設(shè)置越大越能表示細(xì)胞真實(shí)分布拗馒，但也為第二步聚類造成了計(jì)算量激增。

4 層次聚類

4.1 是什么

這是一種歷史悠久的聚類方法溯街，最后會(huì)給出樣本的聚類樹（dendrogram）瘟忱。思路是這樣的：貪婪地將樣本聚集成簇，然后將這些簇聚集成更大的簇苫幢，依此類推访诱，直到所有樣本都屬于一個(gè)簇為止。

層次聚類包含的算法也有很多種韩肝，差別在于如何層層聚集触菜。例如complete linkage方法（hclust的默認(rèn)方法）是保證合并簇之間樣本距離最小， Ward’s方法是保證一個(gè)簇內(nèi)部在聚合的同時(shí)方差增加的最少哀峻，也就是將每次合并時(shí)的信息損失最小化涡相。

在實(shí)際使用中，層次聚類運(yùn)行很慢剩蟀，因此只適用于小型scRNA數(shù)據(jù)催蝗。因?yàn)樾枰?jì)算一個(gè)細(xì)胞間的距離矩陣，如果動(dòng)輒上千細(xì)胞育特，那么這個(gè)計(jì)算量會(huì)很大丙号。

4.2 怎么做

之前用的PBMC數(shù)據(jù)集中細(xì)胞數(shù)太多，于是切換到sce.416b這個(gè)數(shù)據(jù)集缰冤，畢竟一二百個(gè)細(xì)胞還是可以算的

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

首先利用前幾個(gè)PCs計(jì)算細(xì)胞間的距離

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

然后使用Ward‘s方法進(jìn)行聚類

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

準(zhǔn)備畫復(fù)雜一點(diǎn)的聚類樹

library(dendextend)
# seq_along根據(jù)tree.416b$labels犬缨，也就是185個(gè)細(xì)胞的ID，產(chǎn)生了1-185個(gè)數(shù)字棉浸，結(jié)果就是把字符串變數(shù)字
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

# 首先按照不同處理設(shè)置藍(lán)色和紅色怀薛，然后按照批次設(shè)置淡藍(lán)色和淡紅色
labels_colors(dend) <- c(
    `20160113.wild`="blue",
    `20160113.induced`="red",
    `20160325.wild`="dodgerblue",
    `20160325.induced`="salmon"
)[combined.fac][order.dendrogram(dend)]

plot(dend)

圖中可以看到，處理前后確實(shí)形成了紅迷郑、藍(lán)分界

另外更推薦使用Ward’s方法枝恋，它在合并時(shí)受到的cluster差異影響最写淳蟆（畢竟人家的出發(fā)點(diǎn)就是每次合并時(shí)的信息損失最小化，所以如果差異很大焚碌，那我就不合并唄）畦攘。

為了得到更清楚地分群結(jié)果，可以對(duì)樹進(jìn)行”修建“

• 最簡(jiǎn)單的是cutree() 呐能，例如cutree(hc, 4) 手動(dòng)切分?jǐn)?shù)量
• 復(fù)雜一點(diǎn)但更好用的是cutreeDynamic()，來自dynamicTreeCut 這個(gè)包

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)

一般這個(gè)結(jié)果與t-SNE的結(jié)果是吻合的

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

不過這里注意到抑堡，t-SNE中cluster2拆分開了摆出，值得懷疑：是真的聚類出了問題，還是t-SNE本身的顯示問題首妖，畢竟我們知道t-SNE有時(shí)真的會(huì)把一個(gè)亞群拆開畫偎漫。于是繼續(xù)下面??的探索評(píng)價(jià)

4.3 評(píng)價(jià)

可以借助”輪廓圖“（silhouette width）檢查分群的質(zhì)量。會(huì)對(duì)每個(gè)細(xì)胞都計(jì)算一個(gè)silhouette width值有缆，如果一個(gè)細(xì)胞的width值為正并且越大象踊，表示相對(duì)于其他亞群的細(xì)胞，這個(gè)細(xì)胞和它所在亞群中的細(xì)胞更接近棚壁，分群效果越好杯矩；如果width為負(fù)，就表示這個(gè)亞群的這個(gè)細(xì)胞和其他亞群的細(xì)胞更接近袖外，即分群效果不太理想史隆。

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

# 設(shè)置每個(gè)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) 

圖中可以反映的問題：

• 這個(gè)圖是對(duì)上面得到的各個(gè)cluster中的細(xì)胞做的barplot泌射，每個(gè)cluster是一種顏色。
• 這里看到cluster2都是正值并且橫坐標(biāo)width數(shù)值還是最大的鬓照，因此說明了分群沒有問題熔酷，之前t-SNE的結(jié)果的確是它本身的問題
• 如果發(fā)現(xiàn)width值較小，表示分群結(jié)果一般豺裆，還有可能是分群多度拒秘，本來屬于一個(gè)群的，又被拆分成小群
• 利用這個(gè)圖臭猜，我們就能調(diào)整之前的參數(shù)翼抠，來調(diào)整分群效果，不過也不需要太過糾結(jié)完美的分群結(jié)果获讳。因?yàn)榧词箞D上看似合理的分群阴颖，可能實(shí)際上也不會(huì)得到更多的生物信息

5 評(píng)價(jià)分群的穩(wěn)定性

什么叫做分群的穩(wěn)定性？

就是分群之前操作的小小波動(dòng)丐膝，也不會(huì)影響到最后的分群結(jié)果量愧；高穩(wěn)定性的分群結(jié)果其實(shí)也方便了別人進(jìn)行重復(fù)钾菊。scran利用自助法（bootstrapping）進(jìn)行評(píng)價(jià)

在統(tǒng)計(jì)學(xué)中，自助法（Bootstrap Method偎肃，Bootstrapping煞烫，或自助抽樣法）是一種從給定訓(xùn)練集中有放回的均勻抽樣，也就是說累颂，每當(dāng)選中一個(gè)樣本滞详，它等可能地被再次選中并被再次添加到訓(xùn)練集中。 自助法由Bradley Efron于1979年在《Annals of Statistics》上發(fā)表紊馏。

細(xì)胞有放回的抽樣料饥，得到一個(gè)”自助重復(fù)“數(shù)據(jù)集，然后用他進(jìn)行聚類朱监，看看是否能得到相同的亞群

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

originals <- myClusterFUN(sce.pbmc)

# 進(jìn)行自助法判斷
set.seed(0010010100)
coassign <- bootstrapCluster(sce.pbmc, FUN=myClusterFUN, clusters=originals)
dim(coassign)
## [1] 19 19

結(jié)果返回一個(gè)coassignment的矩陣岸啡，表示隨機(jī)從clusterX與clusterY中挑一個(gè)細(xì)胞，這兩個(gè)細(xì)胞在經(jīng)歷”自助法“判斷后赫编，依然落在同一個(gè)cluster的概率巡蘸。我們希望落在對(duì)角線上的概率更高（因?yàn)閷?duì)角線表示自己和自己比較，說明多次重復(fù)后細(xì)胞還是屬于這個(gè)cluster）

# 借助熱圖判斷
pheatmap(coassign, cluster_row=FALSE, cluster_col=FALSE,
    color=rev(viridis::magma(100)))

需要注意的是

• 這個(gè)自助法（Bootstrapping）是一個(gè)通用的評(píng)價(jià)cluster穩(wěn)定性的方法擂送，適用于各種聚類算法
• 高穩(wěn)定性的亞群不一定是高質(zhì)量的悦荒，因?yàn)槿绻粋€(gè)亞群質(zhì)量一直很差，它也很穩(wěn)定

6 數(shù)據(jù)集取小再聚類 | Subclustering

首先通過常規(guī)操作得到分群結(jié)果嘹吨，找到某些感興趣的亞群逾冬，然后在某一個(gè)cluster中再一次進(jìn)行特征基因挑選、聚類操作躺苦，從而增加分析的精確度身腻。

舉個(gè)例子：

首先利用整個(gè)PBMC數(shù)據(jù)集，根據(jù)T細(xì)胞的marker基因篩出可能是T細(xì)胞的亞群匹厘，然后推測(cè)其中某個(gè)亞群可能是記憶T細(xì)胞

# 看到buildSNNGraph就知道是graph-based 畫圖嘀趟，然后用igraph去鑒定
g.full <- buildSNNGraph(sce.pbmc, use.dimred = 'PCA')
clust.full <- igraph::cluster_walktrap(g.full)$membership

# 畫幾個(gè)marker基因的表達(dá)量，使用的是log-normalized表達(dá)量
plotExpression(sce.pbmc, features=c("CD3E", "CCR7", "CD69", "CD44"),
    x=I(factor(clust.full)), colour_by=I(factor(clust.full)))

現(xiàn)在推測(cè)cluster6可能是記憶T細(xì)胞

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

memory <- 6 
sce.memory <- sce.pbmc[,clust.full==memory]
# 對(duì)它重新挑一遍特征基因愈诚，進(jìn)行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進(jìn)行驗(yàn)證
plotExpression(sce.memory, features=c("CD8A", "CD4"),
    x=I(factor(clust.memory)))

這種數(shù)據(jù)集取小再聚類的方法炕柔，可以簡(jiǎn)化亞群的解讀酌泰。只需要在某一個(gè)亞群的背景下繼續(xù)挖掘，就可能會(huì)獲得重要信息匕累，一層層遞進(jìn)陵刹。比如上面的思路：一堆細(xì)胞=》分18群=》鑒定T細(xì)胞的亞群=》取出某個(gè)T細(xì)胞亞群再去鑒定記憶T細(xì)胞

歡迎關(guān)注我們的公眾號(hào)～_～　　
我們是兩個(gè)農(nóng)轉(zhuǎn)生信的小碩，打造生信星球欢嘿，想讓它成為一個(gè)不拽術(shù)語衰琐、通俗易懂的生信知識(shí)平臺(tái)也糊。需要幫助或提出意見請(qǐng)后臺(tái)留言或發(fā)送郵件到jieandze1314@gmail.com

Welcome to our bioinfoplanet!