SC3 : 单细胞转录组聚类分析R包

Kiselev VY, Kirschner K, Schaub MT, Andrews T, Yiu A, Chandra T, Natarajan KN, Reik W, Barahona M, Green AR, Hemberg M (2017). “SC3 - consensus clustering of single-cell RNA-Seq data.” Nature Methods.

一开始我觉得这么火的单细胞数据分析，一定很高大上吧！当我知道单细胞一个主要的功能是揭示细胞异质性，而这个异质性是有聚类算法来实现的。一听聚类，我乐了。关于聚类我能说出一本书啊：

所谓聚类, 就是将一个数据单位的集合分割成几个称为簇或类别的子集 , 每个类中的数据都有相似性,它的划分依据就是“物以类聚”。数据聚类分析是根据事物本身的特性, 研究对被聚类的对象进行类别划分的方法 。聚类分析依据的原则是使同一聚簇中的对象具有尽可能大的相似性, 而不同聚簇中的对象具有尽可能大的相异性, 聚类分析主要解决的问题就是如何在没有先验知识的前提下, 实现满足这种要求的聚簇的聚合。聚类分析称为无监督学习 (Unsuper-vised Study),主要体现在聚类学习的数据对象没有类别标记,需要由聚类学习算法自动计算 。

在我做过的上百场的培训中，每每讲到聚类，我都会提到一句话：说到聚类，有一个概念一定要像思想钢印一样刻在脑子里：距离。不管哪种算法，都需要一个统计量来判断对象之间的远近（或相似性）关系，这个就是距离。同样的对象，不同的距离来度量亲疏不同，就像现在我站你面前，距离很近，但是血缘关系来讲呢，就没有那么近了。

今天我们来认识一种聚类算法（其实是一种组合的聚类算法， Nature Methods，2017）单细胞RNA-seq能够基于转录组特征对细胞类型进行定量表征。我们提出了单细胞一致性聚类(single cell consensus clustering, SC3)，这是一种用户友好的无监督聚类工具，它通过一致方法将多个聚类解决方案组合在一起，从而得到高精度和鲁棒性的分群结果(http://bioconductor.org/packages/SC3)。这中聚类算法允许用户自定义聚类的个数，这个自由不是任何人都消受得起的：一些人根本不知道聚多少类是合适自己的。

我是SC3聚类的一般流程

(a)使用SC3框架进行聚类的概述(参见方法)。用Treutlein数据举例说明了一致步骤。
(b)用于设置SC3参数的已发布数据集。N是数据集中的细胞数;k为作者最初确定的簇数;单位:RPKM是每千碱基每百万次读的转录本，RPM是每百万次读的转录本，FPKM是每千碱基每百万次读的转录本片段，TPM是每百万次读的转录本片段。
(c) ARI>处d值的直方图。金标准数据集达到95。黑色竖线表示细胞总数N的d = 4-7%，分类准确率高。
(d) (b)中所示数据集的SC3聚类的100个实现。条对应点的中位数。红色和灰色分别对应有一致步长和无一致步长时的聚类。这条黑线对应的ARI=0。8。黑色虚线分隔了金和银标准数据集。

文献用到的数据集包含了420个细胞中20000余个基因的表达信息，直接进行聚类资源消耗较大且由于数据噪声过大会导致效果不佳。其实，先降维（特征选择）再聚类已经是通识了。SC3算法的第一步也是要进行基因过滤，文中将表达率低于6%和高于94%的基因都从数据集中去除，将数据集的规模减少了50%。不同的是，SC3第二步会计算细胞与细胞之间的“距离”，为了算法的普适性，这里一共计算了三种距离，分别是欧氏距离，皮尔森相关系数与斯皮尔曼相关系数，最终得到了三个420维的距离矩阵。基于距离矩阵来做pca降维，而后再用k-means聚类。所以SC3与其说优化了聚类算法，不如说是调整了降维的算法。

实现起来是很简单的，因为已经有成熟的R包了：SC3。下载安装就可以用来解锁细胞异质性啦！

卑微小王就不再一步一步跑官网的教程了，SC3也是基于SingleCellExperiment对象，如果是10X的数据可以直接用Seurat过滤QC一下，转化为SingleCellExperiment即可分析。

library(Seurat)
library(SingleCellExperiment)
library(SC3,lib.loc = "D:/R-3.5.1/library")
library(scater,lib.loc = "D:/R-3.5.1/library")
pbmc <- readRDS(file = "D:\\Users\\Administrator\\Desktop\\Novo周运来\\SingleCell\\scrna_tools/pbmc3k_final.rds")

sce <- SingleCellExperiment(
  assays = list(
    counts = as.matrix(pbmc@assays$RNA@counts),
    logcounts = as.matrix(pbmc@assays$RNA@data)
  ), 
  colData = pbmc@meta.data
)

当然，也可以用seurat直接转：

sce <- as.SingleCellExperiment(pbmc)

library(future)
# check the current active plan
plan()
# change the current plan to access parallelization
plan("multiprocess", workers = 4)
plan()
sce <- sc3(sce, ks = 2:10, biology = TRUE) # too time

sc3_plot_expression(
    sce, k = 3, 
    show_pdata = c(
        "cell_type1", 
        "log10_total_features",
        "sc3_3_clusters", 
        "sc3_3_log2_outlier_score"
    )
)

我只是官网示例

SC3聚类算法实现初探
vignettes||SC3
SC3: consensus clustering of single-cell RNA-seq data