比较不同的对单细胞转录组数据聚类的方法

背景介绍

聚类之前必须要对表达矩阵进行normalization,而且要去除一些批次效应等外部因素。通过对表达矩阵的聚类,可以把细胞群体分成不同的状态,解释为什么会有不同的群体。不过从计算的角度来说,聚类还是蛮复杂的,各个细胞并没有预先标记好,而且也没办法事先知道可以聚多少类。尤其是在单细胞转录组数据里面有很高的噪音,基因非常多,意味着的维度很高。

对这样的高维数据,需要首先进行降维,可以选择PCA或者t-SNE方法。聚类的话,一般都是无监督聚类方法,比如:hierarchical clustering, k-means clustering and graph-based clustering。算法略微有一点复杂,略过吧。

这里主要比较6个常见的单细胞转录组数据的聚类包:

  • SINCERA
  • pcaReduce
  • SC3
  • tSNE + k-means
  • SEURAT
  • SNN-Cliq

所以需要安装并且加载一些包,安装代码如下;

install.packages('pcaReduce')
## try http:// if https:// URLs are not supported
source("https://bioconductor.org/biocLite.R") 
biocLite("SC3") 
biocLite("Seurat") 
install.packages("devtools")
library("devtools")
install_github("BPSC","nghiavtr") 
install_github("hemberg-lab/scRNA.seq.funcs")
devtools::install_github("JustinaZ/pcaReduce")

加载代码如下:

library(pcaMethods)
library(pcaReduce)
library(SC3)
library(scater)
library(pheatmap)
set.seed(1234567)

加载测试数据

这里选取的是数据,加载了这个scater包的SCESet对象,包含着一个23730 features, 301 samples 的表达矩阵。

供11已知的种细胞类型,这样聚类的时候就可以跟这个已知信息做对比,看看聚类效果如何。

可以直接用plotPCA来简单PCA并且可视化。

pollen <- readRDS("../pollen/pollen.rds")
pollen
## SCESet (storageMode: lockedEnvironment)
## assayData: 23730 features, 301 samples 
##   element names: exprs, is_exprs, tpm 
## protocolData: none
## phenoData
##   rowNames: Hi_2338_1 Hi_2338_2 ... Hi_GW16_26 (301 total)
##   varLabels: cell_type1 cell_type2 ... is_cell_control (33 total)
##   varMetadata: labelDescription
## featureData
##   featureNames: A1BG A1BG-AS1 ... ZZZ3 (23730 total)
##   fvarLabels: mean_exprs exprs_rank ... feature_symbol (11 total)
##   fvarMetadata: labelDescription
## experimentData: use 'experimentData(object)'
## Annotation:
head(fData(pollen))
##          mean_exprs exprs_rank n_cells_exprs total_feature_exprs
## A1BG     0.56418762      12460            79           169.82048
## A1BG-AS1 0.31265010      10621            37            94.10768
## A1CF     0.05453986       6796            59            16.41650
## A2LD1    0.22572953       9781            28            67.94459
## A2M      0.15087563       8855            21            45.41356
## A2M-AS1  0.02428046       5366             3             7.30842
##          pct_total_exprs pct_dropout total_feature_tpm
## A1BG        1.841606e-03    73.75415            481.37
## A1BG-AS1    1.020544e-03    87.70764            538.18
## A1CF        1.780276e-04    80.39867             13.99
## A2LD1       7.368203e-04    90.69767            350.65
## A2M         4.924842e-04    93.02326           1356.63
## A2M-AS1     7.925564e-05    99.00332             88.61
##          log10_total_feature_tpm pct_total_tpm is_feature_control
## A1BG                    2.683380  1.599256e-04              FALSE
## A1BG-AS1                2.731734  1.787996e-04              FALSE
## A1CF                    1.175802  4.647900e-06              FALSE
## A2LD1                   2.546111  1.164965e-04              FALSE
## A2M                     3.132781  4.507134e-04              FALSE
## A2M-AS1                 1.952356  2.943891e-05              FALSE
##          feature_symbol
## A1BG               A1BG
## A1BG-AS1       A1BG-AS1
## A1CF               A1CF
## A2LD1             A2LD1
## A2M                 A2M
## A2M-AS1         A2M-AS1
table(pData(pollen)$cell_type1)
## 
##   2338   2339     BJ   GW16   GW21 GW21+3  hiPSC   HL60   K562   Kera 
##     22     17     37     26      7     17     24     54     42     40 
##    NPC 
##     15
plotPCA(pollen, colour_by = "cell_type1")

可以看到简单的PCA也是可以区分部分细胞类型的,只不过在某些细胞相似性很高的群体区分力度不够,所以需要开发新的算法来解决这个聚类的问题。

SC聚类

pollen <- sc3_prepare(pollen, ks = 2:5)
pollen <- sc3_estimate_k(pollen)
pollen@sc3$k_estimation
## [1] 11
## 准备 SCESet对象 数据给 SC3方法,先预测能聚多少个类,发现恰好是11个。

## 这里是并行计算,所以速度还可以
pollen <- sc3(pollen, ks = 11, biology = TRUE)
pollen
## SCESet (storageMode: lockedEnvironment)
## assayData: 23730 features, 301 samples 
##   element names: exprs, is_exprs, tpm 
## protocolData: none
## phenoData
##   rowNames: Hi_2338_1 Hi_2338_2 ... Hi_GW16_26 (301 total)
##   varLabels: cell_type1 cell_type2 ... sc3_11_log2_outlier_score
##     (35 total)
##   varMetadata: labelDescription
## featureData
##   featureNames: A1BG A1BG-AS1 ... ZZZ3 (23730 total)
##   fvarLabels: mean_exprs exprs_rank ... sc3_11_de_padj (16 total)
##   fvarMetadata: labelDescription
## experimentData: use 'experimentData(object)'
## Annotation:
head(fData(pollen))
##          mean_exprs exprs_rank n_cells_exprs total_feature_exprs
## A1BG     0.56418762      12460            79           169.82048
## A1BG-AS1 0.31265010      10621            37            94.10768
## A1CF     0.05453986       6796            59            16.41650
## A2LD1    0.22572953       9781            28            67.94459
## A2M      0.15087563       8855            21            45.41356
## A2M-AS1  0.02428046       5366             3             7.30842
##          pct_total_exprs pct_dropout total_feature_tpm
## A1BG        1.841606e-03    73.75415            481.37
## A1BG-AS1    1.020544e-03    87.70764            538.18
## A1CF        1.780276e-04    80.39867             13.99
## A2LD1       7.368203e-04    90.69767            350.65
## A2M         4.924842e-04    93.02326           1356.63
## A2M-AS1     7.925564e-05    99.00332             88.61
##          log10_total_feature_tpm pct_total_tpm is_feature_control
## A1BG                    2.683380  1.599256e-04              FALSE
## A1BG-AS1                2.731734  1.787996e-04              FALSE
## A1CF                    1.175802  4.647900e-06              FALSE
## A2LD1                   2.546111  1.164965e-04              FALSE
## A2M                     3.132781  4.507134e-04              FALSE
## A2M-AS1                 1.952356  2.943891e-05              FALSE
##          feature_symbol sc3_gene_filter sc3_11_markers_clusts
## A1BG               A1BG            TRUE                     5
## A1BG-AS1       A1BG-AS1            TRUE                     4
## A1CF               A1CF            TRUE                     2
## A2LD1             A2LD1           FALSE                    NA
## A2M                 A2M           FALSE                    NA
## A2M-AS1         A2M-AS1           FALSE                    NA
##          sc3_11_markers_padj sc3_11_markers_auroc sc3_11_de_padj
## A1BG            7.740802e-10            0.8554452   1.648352e-18
## A1BG-AS1        1.120284e-03            0.6507985   5.575777e-03
## A1CF            5.007946e-23            0.8592113   1.162843e-17
## A2LD1                     NA                   NA             NA
## A2M                       NA                   NA             NA
## A2M-AS1                   NA                   NA             NA
## 可以看到SC3方法处理后的SCESet对象的基因信息增加了5列,比较重要的是sc3_gene_filter信息,决定着该基因是否拿去聚类,因为基因太多了,需要挑选
table(fData(pollen)$sc3_gene_filter)
## 
## FALSE  TRUE 
## 11902 11828
### 只有一半的基因被挑选去聚类了

## 后面是一些可视化
sc3_plot_consensus(pollen, k = 11, show_pdata = "cell_type1")
sc3_plot_silhouette(pollen, k = 11)
sc3_plot_expression(pollen, k = 11, show_pdata = "cell_type1")
sc3_plot_markers(pollen, k = 11, show_pdata = "cell_type1")
plotPCA(pollen, colour_by = "sc3_11_clusters")
## 还支持shiny的交互式聚类,暂时不显示
# sc3_interactive(pollen)

很明显可以看到SC3聚类的效果要好于普通的PCA

pcaReduce

# use the same gene filter as in SC3
input <- exprs(pollen[fData(pollen)$sc3_gene_filter, ])
# run pcaReduce 1 time creating hierarchies from 1 to 30 clusters
pca.red <- PCAreduce(t(input), nbt = 1, q = 30, method = 'S')[[1]]
##  这里对2~30种类别的情况都分别对样本进行分组。
## 我们这里取只有11组的时候,这些样本是如何分组的信息来可视化。
pData(pollen)$pcaReduce <- as.character(pca.red[,32 - 11])
plotPCA(pollen, colour_by = "pcaReduce")

tSNE + kmeans

scater包包装了 Rtsne 和 ggplot2 来做tSNE并且可视化。

pollen <- plotTSNE(pollen, rand_seed = 1, return_SCESet = TRUE)
## 上面的tSNE的结果,下面用kmeans的方法进行聚类,假定是8类细胞类型。
pData(pollen)$tSNE_kmeans <- as.character(kmeans(pollen@reducedDimension, centers = 8)$clust)
plotTSNE(pollen, rand_seed = 1, colour_by = "tSNE_kmeans")

SNN-Cliq

这个有一点难用,算了吧。

distan <- "euclidean"
par.k <- 3
par.r <- 0.7
par.m <- 0.5
# construct a graph
scRNA.seq.funcs::SNN(
    data = t(input),
    outfile = "snn-cliq.txt",
    k = par.k,
    distance = distan
)
# find clusters in the graph
snn.res <- 
    system(
        paste0(
            "python snn-cliq/Cliq.py ", 
            "-i snn-cliq.txt ",
            "-o res-snn-cliq.txt ",
            "-r ", par.r,
            " -m ", par.m
        ),
        intern = TRUE
    )
cat(paste(snn.res, collapse = "\n"))
snn.res <- read.table("res-snn-cliq.txt")
# remove files that were created during the analysis
system("rm snn-cliq.txt res-snn-cliq.txt")
pData(pollen)$SNNCliq <- as.character(snn.res[,1])
plotPCA(pollen, colour_by = "SNNCliq")

SINCERA

至少是在这个数据集上面表现不咋地

# 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")
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)
## 14
pheatmap(
    t(dat),
    cluster_cols = hc,
    cutree_cols = 14,
    kmeans_k = 100,
    show_rownames = FALSE
)

SEURAT

library(Seurat)
pollen_seurat <- new("seurat", raw.data = get_exprs(pollen, exprs_values = "tpm"))
pollen_seurat <- Setup(pollen_seurat, project = "Pollen")
pollen_seurat <- MeanVarPlot(pollen_seurat)
pollen_seurat <- RegressOut(pollen_seurat, latent.vars = c("nUMI"), 
                            genes.regress = pollen_seurat@var.genes)
pollen_seurat <- PCAFast(pollen_seurat)
pollen_seurat <- RunTSNE(pollen_seurat)
pollen_seurat <- FindClusters(pollen_seurat)
TSNEPlot(pollen_seurat, do.label = T)
pData(pollen)$SEURAT <- as.character(pollen_seurat@ident)
sc3_plot_expression(pollen, k = 11, show_pdata = "SEURAT")
markers <- FindMarkers(pollen_seurat, 2)
FeaturePlot(pollen_seurat, 
            head(rownames(markers)), 
            cols.use = c("lightgrey", "blue"), 
            nCol = 3)

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2018-01-19

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏ATYUN订阅号

【实践】伪造名人的脸—做一个小示例了解生成式对抗网络

生成式对抗网络(GAN)的概念由Ian Goodfellow提出。Goodfellow使用了艺术评论家和艺术家的比喻来描述这两个模型比喻发生器和鉴别,它们组成了...

3484
来自专栏生信技能树

第4篇:对ATAC-Seq/ChIP-seq的质量评估(一)——phantompeakqualtools

在下游分析前,最好是先对peak calling 后的ChIP-Seq数据进行质量评估。

743
来自专栏机器之心

教程 | 从头开始了解PyTorch的简单实现

选自GitHub 机器之心编译 参与:路 本教程展示了如何从了解张量开始到使用 PyTorch 训练简单的神经网络,是非常基础的 PyTorch 入门资源。Py...

3455
来自专栏机器之心

教程 | 利用TensorFlow和神经网络来处理文本分类问题

3197
来自专栏量子位

如何捕获一只彩色卓别林?黑白照片AI上色教程很友好 | 哈佛大触

1112
来自专栏CreateAMind

论文:生成模型采样-类比学习应用 代码

之前发的这篇文章(之前内容在文章底部)介绍了生成模型的高效采样及隐变量空间特征特点,最近的How to Train a GAN? Tips and tricks...

722
来自专栏崔庆才的专栏

从头开始了解PyTorch的简单实现

1565
来自专栏AI科技大本营的专栏

如何让渣画质图片达到逼真效果,试试GAN吧

翻译 | 梁红丽 编辑 | Just 【AI科技大本营导读】在最终视觉呈现效果上,现有的用于极限学习图片压缩的算法似乎都不尽人意,本文作者则使用了 GAN,允许...

3479
来自专栏人工智能头条

Hype:组合机器学习和超参数优化

1728
来自专栏AILearning

【Scikit-Learn 中文文档】分解成分中的信号(矩阵分解问题) - 无监督学习 - 用户指南 | ApacheCN

2.5. 分解成分中的信号(矩阵分解问题) 2.5.1. 主成分分析(PCA) 2.5.1.1. 准确的PCA和概率解释(Exact PCA and p...

2267

扫描关注云+社区