重生之我在剑桥大学学习单细胞RNA-seq分析——6. 生物学分析（2）

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）