我们正常单细胞表达量矩阵,是有2万个基因。但是通常情况下我们走单细胞流程,会仅仅是挑选2000个高变基因做后续分析,所以FindMarkers函数默认就是针对这2000个基因在做分析。
所以如果我们使用FindMarkers函数对两个分组单细胞进行差异分析,它本身默认参数就不可能返回全部的基因,大家可以使用help函数看这个FindMarkers函数的帮助文档,尤其是注意下面的3个参数:
logfc.threshold = 0.25,
test.use = "wilcox",
min.pct = 0.1,
默认用法(针对两个单细胞亚群进行差异分析)如下所示:
markers_2 <- FindMarkers(sce.sub.Endo, ident.1="HE", ident.2="CI")
head(x = markers_2) # 就143个基因
> head(markers_2)
p_val avg_log2FC pct.1 pct.2 p_val_adj
CXCL8 1.920167e-57 0.4110767 0.239 0.485 3.840335e-54
ITPKC 2.792237e-43 0.2522640 0.248 0.414 5.584474e-40
THY1 5.018897e-42 -0.3035701 0.991 0.557 1.003779e-38
FSTL3 6.718088e-37 0.3469674 0.284 0.408 1.343618e-33
IL13RA2 2.964407e-36 0.3333894 0.342 0.490 5.928815e-33
IGFBP5 6.739079e-31 -0.8307001 0.872 0.527 1.347816e-27
也就是说,我们对两个单细胞亚群进行差异分析, 其实就拿到了 143个基因,当然了,本身两个细胞亚群差异可能并不大,因为我们针对的内皮细胞的两个细分亚群,差异本来就不大。但是它其实是可以返回全部的基因,因为仅仅是143个基因我们只能说进行超几何分布检验去注释到go和kegg等数据库,没办法使用gsea分析。
如果是需要gsea分析,一般来说,得拿到全部的基因在两个单细胞亚群的差异情况,而不仅仅是返回统计学显著的差异基因列表。所以全部的使用方法如下所示;
markers <- FindMarkers(sce.sub.Endo , ident.1="HE", ident.2="CI"
assay = 'RNA',slot = 'counts',
logfc.threshold =0,min.pct = 0 )
dim(markers)
> dim(markers)
[1] 21592 5
> head(markers)
p_val avg_log2FC pct.1 pct.2 p_val_adj
RPS4Y1 1.081453e-207 1.478420 0.741 0.004 2.335073e-203
IGHA1 8.766553e-153 1.723276 0.964 0.555 1.892874e-148
MT-CYB 2.117111e-139 -1.508095 0.995 0.997 4.571265e-135
XIST 2.315355e-139 -1.588707 0.000 0.780 4.999314e-135
RPS4X 2.801748e-131 -1.866850 0.957 0.997 6.049533e-127
JCHAIN 3.343925e-124 1.155754 0.745 0.156 7.220203e-120
可以看到,主要是修改了3个参数。
但是,我为什么知道需要修改这3个参数呢?熟能生巧啊,如果你也想学单细胞,参考前面的例子:人人都能学会的单细胞聚类分群注释 ,我们演示了第一层次的分群。
如果你对单细胞数据分析还没有基础认知,可以看基础10讲: