我从Science中，偷学到这个聚类分析技能！真舍不得分享
大家好，我是风。欢迎来到风风的从零开始单细胞系列。前面我们已经学习了数据下载、构建分析对象和数据质控。如果你的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啦，可怕！
赞 (0)