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

背景介绍

聚类之前必须要对表达矩阵进行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 条评论
登录 后参与评论

相关文章

来自专栏WD学习记录

n-gram

N-Gram是大词汇连续语音识别中常用的一种语言模型,对中文而言,我们称之为汉语语言模型(CLM, Chinese Language Model)。汉语语言模型...

773
来自专栏小小挖掘机

DQN三大改进(三)-Dueling Network

论文地址:https://arxiv.org/pdf/1511.06581.pdf 代码地址:https://github.com/princewen/tens...

1K6
来自专栏Ldpe2G的个人博客

图像素描风格生成

论文链接:Combining Sketch and Tone for Pencil Drawing Production

2417
来自专栏Python小屋

Python使用tensorflow中梯度下降算法求解变量最优值

TensorFlow是一个用于人工智能的开源神器,是一个采用数据流图(data flow graphs)用于数值计算的开源软件库。数据流图使用节点(nodes)...

3348
来自专栏人工智能头条

SVM大解密(附代码和公式)

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

如何快速优化机器学习的模型参数

【导读】一般来说机器学习模型的优化没什么捷径可循。用什么架构,选择什么优化算法和参数既取决于我们对数据集的理解,也要不断地试错和修正。所以快速构建和测试模型的能...

722
来自专栏Android点滴积累

Android XML shape 标签使用详解(apk瘦身,减少内存好帮手)

  一个android开发者肯定懂得使用 xml 定义一个 Drawable,比如定义一个 rect 或者 circle 作为一个 View 的背景。但是,也肯...

910
来自专栏Android点滴积累

Android XML shape 标签使用详解(apk瘦身,减少内存好帮手)

Android XML shape 标签使用详解   一个android开发者肯定懂得使用 xml 定义一个 Drawable,比如定义一个 rect 或者 c...

2687
来自专栏数据派THU

一文读懂支持向量机SVM(附实现代码、公式)

支持向量机(SVM),一个神秘而众知的名字,在其出来就受到了莫大的追捧,号称最优秀的分类算法之一,以其简单的理论构造了复杂的算法,又以其简单的用法实现了复杂的问...

1963
来自专栏真皮专栏

聚类算法

p=2时就说平时计算的几何距离,当p趋向于正无穷的时候,其实求的就不是x,y的距离了,而是求x y中最长的一个了。因为如果x大于y,在指数增长下x回远大于y,所...

1182

扫码关注云+社区