6.2 聚类示例
> library(pcaMethods)
> library(SC3)
> library(scater)
> library(SingleCellExperiment)
> library(pheatmap)
> library(mclust)
> set.seed(1234567)
为了说明scRNA-seq数据的聚类，我们使用发育中的小鼠胚胎细胞的Deng数据集（Deng等，2014）。我们已经预处理了数据集并创建了一个SingleCellExperiment对象。我们还用原始文章中确定的细胞类型对细胞进行了注释（colData中的cell_type2）。
6.2.1 Deng数据集
下载链接：https://singlecellcourse.cog.sanger.ac.uk/index.html?shared=data/
加载数据并查看：
> deng <- readRDS("data/deng/deng-reads.rds")
看一下细胞类型注释：
> 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分析已经分离出一些强大的细胞类型并为数据结构提供了一些见解：
早期细胞类型（蓝、橙、绿）分离得很好，但三个囊胚时间点（紫、粉、黄）比较难区分。
6.2.2 SC3
对Deng数据运行SC3聚类。SC3的优势在于它可以直接处理SingleCellExperiment对象。
假设我们不知道cluster的数量k。SC3可以进行估计：
> deng <- sc3_estimate_k(deng)
Estimating k...
> metadata(deng)$sc3$k_estimation
[1] 6
有趣的是，SC3预测的细胞类型数量比原始数据注释中的要少。
但是，如果将不同细胞类型的早期、中期和晚期组合在一起，我们将得到6种细胞类型。我们将合并后的细胞类型存储在colData的cell_type1中：
> plotPCA(deng, colour_by = "cell_type1")
现在我们准备运行SC3（同时要求它计算cluster的生物学特性）：
> 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结果由几种不同的输出组成（更多详细信息请参阅（Kiselev等，2017）和SC3 vignette）。这里我们展示其中的一些：
一致性矩阵：
轮廓图：
表达矩阵热图：
识别marker基因：
SC3聚类的PCA图：
将SC3聚类结果与原文细胞类型标签进行比较：
> adjustedRandIndex(colData(deng)$cell_type2, colData(deng)$sc3_10_clusters)
[1] 0.7804189
SC3也可以在交互式Shiny会话中运行：
> sc3_interactive(deng)
此命令将在浏览器中打开SC3。
注意：由于直接计算距离，当细胞数量>5000时，SC3会变得非常慢。对于数量级在个细胞的大型数据集，我们建议使用Seurat。
6.2.3 tSNE + kmeans
我们之前使用scater包看到的tSNE图是使用Rtsne和ggplot2包制作的。在这里我们将做同样的操作：
> deng <- runTSNE(deng, rand_seed = 1)
> plotTSNE(deng)
患者数据的tSNE图
注意，上图所有点都是黑色的。这与我们之前看到的不同，当时细胞的颜色是基于注释的。这里我们没有任何注释，所有细胞都来自同一批，因此所有点都是黑色的。
现在我们对tSNE图上的点应用k均值聚类算法。你看到了多少组？
我们将从k=8开始：
> colData(deng)$tSNE_kmeans <- as.character(kmeans(reducedDim(deng, "TSNE"), centers = 8)$clust)
> plotTSNE(deng, colour_by = "tSNE_kmeans")
患者数据的tSNE图，包含8个聚类，由k均值聚类算法识别。
你可能已经注意到，tSNE+kmeans是随机的，每次运行时都会产生不同的结果。为了达到更好的效果，我们需要多次运行。SC3也是随机的，但由于一致性步骤，它更加稳健，不太可能产生不同的结果。
6.2.4 SINCERA
如上一章所述，SINCERA基于层次聚类，它在进行聚类之前会执行z-score转换：
> 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")
如果不知道聚类的数量，SINCERA可以将生成不超过特定数量的单例聚类（仅包含1个细胞的聚类）的层次树的最小高度设为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结果可视化为热图：
使用k的SINCERA方法聚类结果
往期内容：
重生之我在剑桥大学学习单细胞RNA-seq分析——5. scRNA-seq数据的基本质量控制 (QC) 和探索（5）
重生之我在剑桥大学学习单细胞RNA-seq分析——5. scRNA-seq数据的基本质量控制 (QC) 和探索（6）
重生之我在剑桥大学学习单细胞RNA-seq分析——6. 生物学分析（1）