我从Science中，偷学到这个聚类分析技能！真舍不得分享

东西二王 2021-03-29

展开全文

大家好，我是风。欢迎来到风风的从零开始单细胞系列。前面我们已经学习了数据下载、构建分析对象和数据质控。如果你的scater出现了一些警告内容，提示函数被替代，那也不用着急。正如其他内容一样，scater可以进行数据质控，那么Seurat一样可以。为了方便大家学习，今天我们先接着scater的流程往下分析，Seurat进行数据质控的内容我们在后面介绍Seurat流程时再统一学习。

今天我们继续新的内容——聚类cluster。话不多说，直接来看一篇2014年发表在Science杂志上的经典文章：Single-Cell RNA-Seq Reveals Dynamic, Random Monoallelic Gene Expression in Mammalian Cells，老规矩，原文还是在后台回复关键词获取。（咦？好像这次Science之后，风风系列就有了Nature, Cell和Science了(＃°Д°)）

我们来看看文章的这个图：

这个图是单细胞数据PCA分析的可视化，那今天我们就从这个图开始拓展，看看单细胞分析中的聚类和PCA可以玩出什么花样吧，比如，像下面这个图：

期刊信息

读取数据

这个步骤基于scater对象，关于如何构建scater对象的内容可以看前面的推文，我们不再重复，后台数据已经使用原文数据构建scater对象完毕，这里直接加载进来分析：

options(stringsAsFactors = FALSE)set.seed(20210322) # 设置随机种子，方便重复结果# 加载所需要的R包，如果没装可以直接使用Biocmanager安装library(pcaMethods)library(scater)library(SingleCellExperiment)library(pheatmap)library(mclust)# 加载数据scRNA <- readRDS('data.rds')scRNA # 查看数据整体情况## class: SingleCellExperiment ## dim: 22431 268 ## metadata(0):## assays(2): counts logcounts## rownames(22431): Hvcn1 Gbp7 ... Sox5 Alg11## rowData names(10): feature_symbol is_feature_control ... total_counts##   log10_total_counts## colnames(268): 16cell 16cell.1 ... zy.2 zy.3## colData names(30): cell_type2 cell_type1 ... pct_counts_ERCC##   is_cell_control## reducedDimNames(0):## spikeNames(1): ERCC## altExpNames(0):# 查看细胞类型的注释Cell <- table(colData(scRNA)$cell_type2)Cell # 总共有10种细胞类型## ##     16cell      4cell      8cell early2cell earlyblast  late2cell  lateblast ##         50         14         37          8         43         10         30 ##   mid2cell   midblast         zy ##         12         60          4# 我们直接使用plotPCA函数进行PCA分析plotPCA(scRNA, colour_by = 'cell_type2')## Warning: call 'runPCA' explicitly to compute results

好了，现在一个简单的PCA分析已经完成，我们可以看到早期的细胞全部聚在一起，聚类结果较好，但是有一些其他类型的细胞没有很好地聚类,那该怎么办呢？我们顺着思路往下走，既然原始数据的PCA结果不好，我们给它重新进行聚类分析，然后再PCA进行降维，看看结果是否符合我们要求，是不是就可以了呢？好像是这样没错，但是紧接着第二个问题又来了，该怎么进行聚类分析？

单细胞数据的聚类分析跟二代测序数据的聚类分析一致，有kmeans、PAM等方法，使用的R包像常见的ConsensusCluster或者NMF都可以，这里我们用一个专属于单细胞数据分析的R包——SC3包。SC3包进行聚类非常方便，常规使用Kmeans我们需要自己判断数据可以区分为几个亚类，但是SC3包的sc3_estimate_k函数可以帮我们自动聚类，选取最佳cluster数量：

library(SC3)scRNA1 <- sc3_estimate_k(scRNA) # 对构建的scater对象进行kmeas聚类## Estimating k...metadata(scRNA1)$sc3$k_estimation # 查看最佳聚类的数目，发现为6类## [1] 6# 结果发现SC3聚类识别到的聚类数目和原始数据的细胞类型大相径庭，SC3认为6类就是最佳聚类，这个数量比我们前面看到的10类要少得多，那到底哪个更加准确呢？口说无凭，我们进行PCA分析看一下结果：plotPCA(scRNA1, colour_by = 'cell_type1')## Warning: call 'runPCA' explicitly to compute results

# 从PCA分析的结果来看就很明显了，SC3进行kmeas聚类后的结果明显更好，把各时期的细胞都很好地区分出来# 那如果我不需要SC3自动获取最佳聚类，或者SC3自动聚类的结果不好，我们想要自己获得最佳聚类该怎么办呢？我们可以直接给SC3指定聚类的数目，这次用的是sc3函数：scRNA2 <- sc3(scRNA, # 分析对象             ks = 10, # 聚类数目             biology = TRUE,  # 定义是否计算差异表达基因，标记基因和细胞异常值             gene_filter = TRUE, # 选择使用基因过滤             pct_dropout_min = 10, # 去除百分比小于10%的基因             pct_dropout_max = 90, # 去除百分比大于90%的基因             n_cores = 4 # 定义使用的电脑核数，我的电脑是8核，这里使用4核来加快运行速度             )# sc3运行需要一小段时间，如果电脑配置稍微差点，需要的时间会更多，运行完毕后，我们来查看下结果：# 查看Consensus matrixsc3_plot_consensus(scRNA2, k = 10, show_pdata = 'cell_type2')

# 发现确实分为10类，但是有的类特别少# 轮廓图sc3_plot_silhouette(scRNA2, k = 10)

# 展示聚类矩阵热图sc3_plot_expression(scRNA2, k = 10, show_pdata = 'cell_type2')

# 识别标记基因sc3_plot_markers(scRNA2, k = 10, show_pdata = 'cell_type2')

# 再次进行PCA分析plotPCA(scRNA2, colour_by = 'sc3_10_clusters')## Warning: call 'runPCA' explicitly to compute results

# 对比原始的PCA图plotPCA(scRNA, colour_by = 'cell_type2')## Warning: call 'runPCA' explicitly to compute results

# 我们也可以直接提取SC3聚类结果sc3_10_cluster <- colData(scRNA2)$sc3_10_clusterssc3_10_cluster##   [1] 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  4  1  7  1  1  1  1  1  1  1 ##  [26] 1  1  1  1  1  1  1  4  1  1  1  1  1  1  1  1  1  1  1  1  4  1  1  1  1 ##  [51] 5  5  5  5  5  5  5  5  5  5  5  5  5  5  1  1  1  1  1  1  1  1  1  1  1 ##  [76] 1  1  1  1  1  1  1  1  1  1  1  1  4  1  1  1  1  4  4  4  4  4  4  4  4 ## [101] 4  6  6  6  6  6  6  6  6  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8 ## [126] 7  7  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8  8 ## [151] 8  8  3  3  3  3  3  3  3  3  3  3  2  2  2  2  2  2  2  2  2  2  10 2  10## [176] 2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  2  3  3  3  3  3  3  3  3 ## [201] 3  3  3  3  10 9  9  10 10 9  9  9  10 9  10 10 10 9  10 10 9  10 10 10 10## [226] 10 10 10 7  9  10 10 10 10 9  10 9  10 10 10 9  10 9  9  10 10 9  9  10 9 ## [251] 10 9  10 10 9  9  9  9  9  9  9  9  9  9  6  6  6  6 ## Levels: 1 2 3 4 5 6 7 8 9 10# 再看一下原始数据的类别rawtype <- colData(scRNA)$cell_type2rawtype##   [1] 16cell     16cell     16cell     16cell     16cell     16cell    ##   [7] 16cell     16cell     16cell     16cell     16cell     16cell    ##  [13] 16cell     16cell     16cell     16cell     16cell     16cell    ##  [19] 16cell     16cell     16cell     16cell     16cell     16cell    ##  [25] 16cell     16cell     16cell     16cell     16cell     16cell    ##  [31] 16cell     16cell     16cell     16cell     16cell     16cell    ##  [37] 16cell     16cell     16cell     16cell     16cell     16cell    ##  [43] 16cell     16cell     16cell     16cell     16cell     16cell    ##  [49] 16cell     16cell     4cell      4cell      4cell      4cell     ##  [55] 4cell      4cell      4cell      4cell      4cell      4cell     ##  [61] 4cell      4cell      4cell      4cell      8cell      8cell     ##  [67] 8cell      8cell      8cell      8cell      8cell      8cell     ##  [73] 8cell      8cell      8cell      8cell      8cell      8cell     ##  [79] 8cell      8cell      8cell      8cell      8cell      8cell     ##  [85] 8cell      8cell      8cell      8cell      8cell      8cell     ##  [91] 8cell      8cell      8cell      8cell      8cell      8cell     ##  [97] 8cell      8cell      8cell      8cell      8cell      early2cell## [103] early2cell early2cell early2cell early2cell early2cell early2cell## [109] early2cell earlyblast earlyblast earlyblast earlyblast earlyblast## [115] earlyblast earlyblast earlyblast earlyblast earlyblast earlyblast## [121] earlyblast earlyblast earlyblast earlyblast earlyblast earlyblast## [127] earlyblast earlyblast earlyblast earlyblast earlyblast earlyblast## [133] earlyblast earlyblast earlyblast earlyblast earlyblast earlyblast## [139] earlyblast earlyblast earlyblast earlyblast earlyblast earlyblast## [145] earlyblast earlyblast earlyblast earlyblast earlyblast earlyblast## [151] earlyblast earlyblast late2cell  late2cell  late2cell  late2cell ## [157] late2cell  late2cell  late2cell  late2cell  late2cell  late2cell ## [163] lateblast  lateblast  lateblast  lateblast  lateblast  lateblast ## [169] lateblast  lateblast  lateblast  lateblast  lateblast  lateblast ## [175] lateblast  lateblast  lateblast  lateblast  lateblast  lateblast ## [181] lateblast  lateblast  lateblast  lateblast  lateblast  lateblast ## [187] lateblast  lateblast  lateblast  lateblast  lateblast  lateblast ## [193] mid2cell   mid2cell   mid2cell   mid2cell   mid2cell   mid2cell  ## [199] mid2cell   mid2cell   mid2cell   mid2cell   mid2cell   mid2cell  ## [205] midblast   midblast   midblast   midblast   midblast   midblast  ## [211] midblast   midblast   midblast   midblast   midblast   midblast  ## [217] midblast   midblast   midblast   midblast   midblast   midblast  ## [223] midblast   midblast   midblast   midblast   midblast   midblast  ## [229] midblast   midblast   midblast   midblast   midblast   midblast  ## [235] midblast   midblast   midblast   midblast   midblast   midblast  ## [241] midblast   midblast   midblast   midblast   midblast   midblast  ## [247] midblast   midblast   midblast   midblast   midblast   midblast  ## [253] midblast   midblast   midblast   midblast   midblast   midblast  ## [259] midblast   midblast   midblast   midblast   midblast   midblast  ## [265] zy         zy         zy         zy        ## 10 Levels: 16cell 4cell 8cell early2cell earlyblast late2cell ... zy# 比较两种聚类结果是否一致，这里使用rand指数adjustedRandIndex(sc3_10_cluster, rawtype)## [1] 0.6620724# 结果约为0.66，说明两种聚类具有一定的相似度# 当然如果你觉得这样还是太麻烦，那我们也可以直接交互式处理，在网站上进行点击sc3_interactive(scRNA2)

打开网页后看到如下页面：

左边可以勾选热图的注释信息，中间是可视化热图，右边则是描述的内容，可以作为文章的Figure Legend，当然还有其他图，比如轮廓图：

识别了Marker的热图：

基本可以靠点击鼠标完成个性化分析，并且还直接出来Figure legend，连写作都搞定了，是不是觉得非常贴心！别急，还没结束，有了PCA，我们再将大名鼎鼎的tSNE和层次聚类结合：

# 进行tSNE分析scRNA3 <- runTSNE(scRNA,                 scale = TRUE,                 perplexity = 30,                 feature_set = metadata(scRNA3)$hvg_genes,                 set.seed = 1)# 可视化plotReducedDim(scRNA3, 'TSNE', colour_by = 'cell_type2')

# 使用以前推文中学过的层次聚类distance <- dist(t(logcounts(scRNA)))ward <- hclust(distance,method = 'ward.D')# 跟前面一样直接设置为10个亚群hclust_10_cluster <- cutree(ward, k = 10)colData(scRNA3)$hclust_10_cluster <- factor(hclust_10_cluster)plotReducedDim(scRNA3, 'TSNE', colour_by = 'hclust_10_cluster')

搞定！今天的聚类的内容就到这里，内容承接前面的scater对象。推文的难度也逐渐增加，从刚开始的数据下载到现在进行细胞聚类，一篇简单的5 单细胞分析的文章就已经完成了1/3了，你别不信，上次好几位学员发了自己感兴趣的5分左右的文章给我作为复现文章备选，按照那几个文章的工作量，今天的内容掌握之后，5分的文章就完成1/3啦，可怕！