日常瞎掰
寻找细胞亚群是单细胞的一个重要的基本任务，后续的分析基本都是基于这些分好的细胞亚群进行的。可见细胞分群是多么的重要，这就有点成也萧何，败也萧何的感觉了！能找到一个比较合理的分群方式才会有后续的分析内容。
当然了，今天我们并不是要讨论如何分群，今天我们想说一下如果分群时有些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，乏了，我要去躺平了，今天就说这么多了！
网友评论