日常瞎掰
寻找细胞亚群是单细胞的一个重要的基本任务,后续的分析基本都是基于这些分好的细胞亚群进行的。可见细胞分群是多么的重要,这就有点成也萧何,败也萧何的感觉了!能找到一个比较合理的分群方式才会有后续的分析内容。
当然了,今天我们并不是要讨论如何分群,今天我们想说一下如果分群时有些cluster找不到marker时该怎么办?能不能分好群其实主要受制于以下两个因素:
1、研究方面的marker多不多,这个是直接影响分群的因素,没有marker基因何谈将cluster分为某个细胞类型;
2、数据质量好不好,这里说的数据质量包括两层含义:一是数据本身的质量,例如测序质量好不好,这个直接关乎数据能不能使用,二是数据生物学质量,不是所有测序质量好的数据都能得出很有意义的结果,这么说也许大家都可以理解了,想要好的且有意义的分析结果关键还是在于实验做的好不好!
第一个因素算是技术上的限制,也许随着后续研究的发展,更多的marker基因会被发现,这个没什么好说,只能基于现有的marker基因分析。第二因素就不是那么容易界定了,毕竟测序数据都有了,测序质量也没有问题,这时如果说数据有问题,也许最难接受的恐怕是做实验的人了。
说到这里让我想起来一个不是很恰当的但具有讽刺意义的笑话:证明蜘蛛的听力在脚上。为了证明蜘蛛的听力在脚上,有一个研究人员先抓来一只蜘蛛,然后朝着蜘蛛吼了一声,蜘蛛被吓跑了。接着把那只蜘蛛再抓回来,然后把蜘蛛的脚都拔了,再朝着他吼一声,发现蜘蛛不跑了。于是研究人员得出结果,蜘蛛的听力在脚上!
好吧,说了这么一大堆大家都知道的废话,下面我回到今天的主题,当找不到cluster的marker基因时,我们应该怎么办呢?当然是查找一下可能的原因,这时不妨再细看一下这两个QC的图,也许可以找到一些线索。
示例展示
这里用一个Seurat对象做为演示,仅是为了做一个展示,代码不是很完整:
>library(Seurat)
>library(dplyr)
>library(scater)
>sceagg
An object of class Seurat
36615 features across 69000 samples within 1 assay
Active assay: RNA (36615 features, 0 variable features)
>sce <- as.SingleCellExperiment(sceagg)
>mt <- grepl("^MT-", rownames(sce))
>sce <- addPerCellQC(sce, subsets=list(Mito=mt))
>sce <- addPerFeatureQC(sce, subsets=list(Mito=mt))
>scemeta <- as.data.frame(colData(sce))
>options(repr.plot.width = 8, repr.plot.height=7)
>ggplot(scemeta, aes(ident,sum,colour=ident)) + geom_violin()
>ggplot(scemeta, aes(ident,detected,colour=ident)) + geom_violin()
每个细胞nCount分布图:
每个细胞nCount检测到基因数的分布图:
上面这两张图,本质上跟Seurat用nFeature、nCount画的点图类似,只不过小提琴可以看出数量分布特征,而当我们分好cluster以后再来画一个这样的图看看,就会得到一些额外的信息。例如,这里共分成了15个cluster,基本上所有的cluster中绝大部分细胞的count数都小于5000,其中0、1、3、5、9、12这些cluster中大部分细胞的count数在1250左右。再来看看下图张图,基本上所有的cluster中绝大部分细胞表达的基因数都小于2000,其中0、1、3、5、9、12这些cluster中大部分细胞表达的基因数在750左右。从这两张图,我们可以想到一些什么呢?很想听听大家的看法。
结束语
Seurat基本上已经是分析细胞数据的明星软件了,友好到没话说,但这里还是想说一下scater这个包,这个包用来做QC还是值得一用的,比如addPerCellQC,addPerFeatureQC这两个函数就可以统计一下细胞层面的信息,还可以用来统计ERCC的信息,还有可以一次性统计所有这些信息的函数,功能太多,这里就不详细说明了,大家可以自行摸索。emm,乏了,我要去躺平了,今天就说这么多了!