6.2 聚類示例
> library(pcaMethods)
> library(SC3)
> library(scater)
> library(SingleCellExperiment)
> library(pheatmap)
> library(mclust)
> set.seed(1234567)
為了說明scRNA-seq數(shù)據(jù)的聚類毕箍，我們使用發(fā)育中的小鼠胚胎細(xì)胞的Deng數(shù)據(jù)集（Deng等，2014）道盏。我們已經(jīng)預(yù)處理了數(shù)據(jù)集并創(chuàng)建了一個(gè)SingleCellExperiment對象而柑。我們還用原始文章中確定的細(xì)胞類型對細(xì)胞進(jìn)行了注釋（colData中的cell_type2）。
6.2.1 Deng數(shù)據(jù)集
下載鏈接：https://singlecellcourse.cog.sanger.ac.uk/index.html?shared=data/
加載數(shù)據(jù)并查看：
> deng <- readRDS("data/deng/deng-reads.rds")
看一下細(xì)胞類型注釋：
> table(colData(deng)$cell_type2)
16cell 4cell 8cell early2cell earlyblast late2cell lateblast
50 14 37 8 43 10 30
mid2cell midblast zy
12 60 4
簡單的PCA分析已經(jīng)分離出一些強(qiáng)大的細(xì)胞類型并為數(shù)據(jù)結(jié)構(gòu)提供了一些見解：
早期細(xì)胞類型（藍(lán)荷逞、橙媒咳、綠）分離得很好，但三個(gè)囊胚時(shí)間點(diǎn)（紫种远、粉涩澡、黃）比較難區(qū)分。
6.2.2 SC3
對Deng數(shù)據(jù)運(yùn)行SC3聚類坠敷。SC3的優(yōu)勢在于它可以直接處理SingleCellExperiment對象妙同。
假設(shè)我們不知道cluster的數(shù)量k。SC3可以進(jìn)行估計(jì)：
> deng <- sc3_estimate_k(deng)
Estimating k...
> metadata(deng)$sc3$k_estimation
[1] 6
有趣的是膝迎，SC3預(yù)測的細(xì)胞類型數(shù)量比原始數(shù)據(jù)注釋中的要少粥帚。
但是，如果將不同細(xì)胞類型的早期限次、中期和晚期組合在一起芒涡，我們將得到6種細(xì)胞類型。我們將合并后的細(xì)胞類型存儲在colData的cell_type1中：
> plotPCA(deng, colour_by = "cell_type1")
現(xiàn)在我們準(zhǔn)備運(yùn)行SC3（同時(shí)要求它計(jì)算cluster的生物學(xué)特性）：
> deng <- sc3(deng, ks = 10, biology = TRUE, n_cores = 1)
Setting SC3 parameters...
Calculating distances between the cells...
Performing transformations and calculating eigenvectors...
Performing k-means clustering...
Calculating consensus matrix...
Calculating biology...
SC3結(jié)果由幾種不同的輸出組成（更多詳細(xì)信息請參閱（Kiselev等卖漫，2017）和SC3 vignette）费尽。這里我們展示其中的一些：
一致性矩陣：
輪廓圖：
表達(dá)矩陣熱圖：
識別marker基因：
SC3聚類的PCA圖：
將SC3聚類結(jié)果與原文細(xì)胞類型標(biāo)簽進(jìn)行比較：
> adjustedRandIndex(colData(deng)$cell_type2, colData(deng)$sc3_10_clusters)
[1] 0.7804189
SC3也可以在交互式Shiny會話中運(yùn)行：
> sc3_interactive(deng)
此命令將在瀏覽器中打開SC3。
注意：由于直接計(jì)算距離羊始，當(dāng)細(xì)胞數(shù)量>5000時(shí)旱幼，SC3會變得非常慢。對于數(shù)量級在個(gè)細(xì)胞的大型數(shù)據(jù)集店枣，我們建議使用
Seurat速警。
6.2.3 tSNE + kmeans
我們之前使用scater包看到的tSNE圖是使用Rtsne和ggplot2包制作的。在這里我們將做同樣的操作：
> deng <- runTSNE(deng, rand_seed = 1)
> plotTSNE(deng)
注意鸯两，上圖所有點(diǎn)都是黑色的闷旧。這與我們之前看到的不同，當(dāng)時(shí)細(xì)胞的顏色是基于注釋的钧唐。這里我們沒有任何注釋忙灼，所有細(xì)胞都來自同一批，因此所有點(diǎn)都是黑色的。
現(xiàn)在我們對tSNE圖上的點(diǎn)應(yīng)用k均值聚類算法该园。你看到了多少組酸舍？
我們將從k=8開始：
> colData(deng)$tSNE_kmeans <- as.character(kmeans(reducedDim(deng, "TSNE"), centers = 8)$clust)
> plotTSNE(deng, colour_by = "tSNE_kmeans")
你可能已經(jīng)注意到，tSNE+kmeans是隨機(jī)的双妨，每次運(yùn)行時(shí)都會產(chǎn)生不同的結(jié)果淮阐。為了達(dá)到更好的效果，我們需要多次運(yùn)行刁品。SC3也是隨機(jī)的泣特，但由于一致性步驟，它更加穩(wěn)健挑随，不太可能產(chǎn)生不同的結(jié)果状您。
6.2.4 SINCERA
如上一章所述，SINCERA基于層次聚類兜挨，它在進(jìn)行聚類之前會執(zhí)行z-score轉(zhuǎn)換：
> input <- logcounts(deng[rowData(deng)$sc3_gene_filter, ])
# perform gene-by-gene per-sample z-score transformation
> dat <- apply(input, 1, function(y) scRNA.seq.funcs::z.transform.helper(y))
# hierarchical clustering
> dd <- as.dist((1 - cor(t(dat), method = "pearson"))/2)
> hc <- hclust(dd, method = "average")
如果不知道聚類的數(shù)量膏孟，SINCERA可以將生成不超過特定數(shù)量的單例聚類（僅包含1個(gè)細(xì)胞的聚類）的層次樹的最小高度設(shè)為k：
> num.singleton <- 0
> kk <- 1
> for (i in 2:dim(dat)[2]) {
clusters <- cutree(hc, k = i)
clustersizes <- as.data.frame(table(clusters))
singleton.clusters <- which(clustersizes$Freq < 2)
if (length(singleton.clusters) <= num.singleton) {
kk <- i
} else {
break;
}
}
> cat(kk)
6
將SINCERA結(jié)果可視化為熱圖：
往期內(nèi)容：
重生之我在劍橋大學(xué)學(xué)習(xí)單細(xì)胞RNA-seq分析——5. scRNA-seq數(shù)據(jù)的基本質(zhì)量控制 (QC) 和探索（5）
重生之我在劍橋大學(xué)學(xué)習(xí)單細(xì)胞RNA-seq分析——5. scRNA-seq數(shù)據(jù)的基本質(zhì)量控制 (QC) 和探索（6）
重生之我在劍橋大學(xué)學(xué)習(xí)單細(xì)胞RNA-seq分析——6. 生物學(xué)分析（1）
