我从Science中,“偷”到了这个实用的聚类分析技能!真舍不得分享!(附代码)
一键学习Scater流程
大家好,我是风。欢迎来到风风的从零开始单细胞系列。前面我们已经学习了数据下载、构建分析对象和数据质控。如果你的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...
$k_estimation# 查看最佳聚类的数目,发现为6类 sc3
# [1] 6
结果发现SC3聚类识别到的聚类数目和原始数据的细胞类型大相径庭,SC3认为6类就是最佳聚类,这个数量比我们前面看到的10类要少得多,那到底哪个更加准确呢?口说无凭,我们进行PCA分析看一下结果:
plotPCA(scRNA1, colour_by = "cell_type1")
# Warning: call 'runPCA' explicitly to compute results
scRNA2 <- sc3(scRNA,
ks = 10,
biology = TRUE,
gene_filter = TRUE,
pct_dropout_min = 10,
pct_dropout_max = 90,
n_cores = 4
)
sc3_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_clusters
sc3_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_type2
rawtype
# [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啦,可怕!
好啦,今天的内容就到这里,后台回复“feng43”获取数据和代码,我们下期见啦!
单细胞分析专栏传送门
碎碎念专栏传送门(完结)
风之美图系列传送门(完结)
—
END—
撰文丨风 风
排版丨四金兄
值班 | 风 风
主编丨小雪球
欢迎大家关注解螺旋生信频道-挑圈联靠公号~
社群的小伙伴们看过来!
-首次酸谈社群日-来了!
---酸谈社群日---
聊聊酸谈社群的成长历程
聊聊大家学习中遇到的问题
还有超多福利好礼回馈大家
---直播---
直播时间:
3月27日(周六)晚18点-20点直播内容:
《酸谈社群日》专场直播
带着你的小伙伴们
一起相约解螺旋直播间吧!
领悟科研 优人一步
最新评论
推荐文章
作者最新文章
你可能感兴趣的文章
Copyright Disclaimer: The copyright of contents (including texts, images, videos and audios) posted above belong to the User who shared or the third-party website which the User shared from. If you found your copyright have been infringed, please send a DMCA takedown notice to [email protected]. For more detail of the source, please click on the button "Read Original Post" below. For other communications, please send to [email protected].
版权声明:以上内容为用户推荐收藏至CareerEngine平台,其内容(含文字、图片、视频、音频等)及知识版权均属用户或用户转发自的第三方网站,如涉嫌侵权,请通知[email protected]进行信息删除。如需查看信息来源,请点击“查看原文”。如需洽谈其它事宜,请联系[email protected]。
版权声明:以上内容为用户推荐收藏至CareerEngine平台,其内容(含文字、图片、视频、音频等)及知识版权均属用户或用户转发自的第三方网站,如涉嫌侵权,请通知[email protected]进行信息删除。如需查看信息来源,请点击“查看原文”。如需洽谈其它事宜,请联系[email protected]。