目前主流的10x平台的单细胞转录组测序后得到的表达量矩阵里面的每个样品都是成百上千个细胞,因为技术本身就是牺牲了每个细胞的基因数量来换取更多的单细胞的产出,这样的话,每个细胞仅仅是测到一两千个基因就足够了,换句话说,就是表达量矩阵里面有非常多的0值,绝大部分基因在绝大部分细胞都是没有检测到的,也就是表达量为0,这个是单细胞技术的天然缺陷, 有一个专有名词:drop-out
单细胞表达量矩阵里面的0值非常多
我们以pbmc3k数据集为例:
library(Seurat)
library(SeuratData)
#InstallData("pbmc3k")
data("pbmc3k")
sce <- pbmc3k.final
sce@assays$RNA@counts[1:4,1:4]
可以看到,前面的4个细胞的4个基因都是0,在稀疏矩阵里面的0以小数点表示,如下所示是:
> sce@assays$RNA@counts[1:4,1:4]
4 x 4 sparse Matrix of class "dgCMatrix"
AAACATACAACCAC AAACATTGAGCTAC AAACATTGATCAGC AAACCGTGCTTCCG
AL627309.1 . . . .
AP006222.2 . . . .
RP11-206L10.2 . . . .
RP11-206L10.9 . . . .
>
也很容易计算到这个数据集到底是多少个0
> dim(sce@assays$RNA@counts)
[1] 13714 2638
> length(sce@assays$RNA@counts)
[1] 36177532
> sum(sce@assays$RNA@counts==0)
[1] 33938800
> sum(sce@assays$RNA@counts==0)/length(sce@assays$RNA@counts)
[1] 0.9381182
CD79A 和 CD79B都是B细胞的标记基因,所以相关性很好
既然具体到每个细胞来看,绝大部分基因都是0值,这样的话不同细胞之间很难计算相关性,比如FeatureScatter 函数专门是干这个事 :
library(patchwork)
FeatureScatter(sce,"nCount_RNA" ,"nFeature_RNA") +
FeatureScatter(sce,"CD14" ,"CD4") +
FeatureScatter(sce,"CD79A" ,"CD79B")
如下所示:
因为 nCount_RNA 和 nFeature_RNA是细胞的熟悉,所以没有0的干扰,这个相关性很好,而且是可靠的。
另外,因为 CD14 和 CD4 本来是髓系免疫细胞和cd4T细胞的标记基因,理论上就相关性应该是很差。
最后,CD79A 和 CD79B都是B细胞的标记基因,他们的相关性确实是应该是很好。
但是CD79A 和 CD79B在b细胞亚群里面是没有相关性的
看起来一切合情合理,但是如果我们具体到B细胞本身,就发现不对劲了。
Bsce = sce[,Idents(sce) == 'B']
FeatureScatter(Bsce,"nCount_RNA" ,"nFeature_RNA") +
FeatureScatter(Bsce,"CD79A" ,"CD79B")+
FeatureScatter(Bsce,"CD79A" ,"MS4A1")
他们的相关性居然完全没有了。。。。
这个时候有两个解释,首先是因为0值的存在,影响了相关性技术,其次是因为它们虽然都是B细胞的标记基因仅仅是说明它们都是应该在B细胞亚群里面高表达,并不能推理出来它们应该是正相关。
当然了,单细胞水平不同基因的表达量相关性本来就不应该是如此简单的计算,不过这个简单的探索,这两个简单的推理还是蛮有意思的的。
天色已晚,我不想写了,亲爱的读者们大家觉得应该是哪种可能性呢?