这算是不一样吗

在《单细胞天地》的交流群，看到有粉丝提问关于 FindMarkers与AverageExpression 两个函数的差异 ：

请问下怎么得到FindMarkers计算时ident.1和ident.2指定的cluster的平均表达量，因为输出结果中只会给avg_log2FC，我知道AverageExpression()函数可以计算平均表达量，但对应计算出来的logFC对应不上！都弄了2天了，一直找不到原因和解决办法，求大神解答！

而且还很贴心的给出来了一个示意图，如下所示：

 

这样的提问太棒了，所以我就让他把这个数据save成为rdata文件，然后把代码复制粘贴到交流群里面。我马上开始debug模式！

我的代码如下所示：

library(Seurat)
load('pbmc_Platelet_DC.Rdata')
sce=pbmc_Platelet_DC
table(sce$seurat_clusters)
table(Idents(sce))
#计算平均表达量
av <-AverageExpression(sce, assays = "RNA") 
av=as.data.frame(av$RNA)
av$diff = log(av[,2]+0.1) - log(av[,1]+0.1)

#计算差异基因
diff  <- FindMarkers(sce,
                     ident.1 = "Platelet", 
                     ident.2 = "DC",assay = "RNA")

head(diff)
head(av)

很快就能重复出来 FindMarkers与AverageExpression 两个函数的结果 ：

> head(diff)
                 p_val avg_logFC pct.1 pct.2    p_val_adj
PPBP      1.511444e-10  6.028905 1.000 0.031 2.072795e-06
PF4       3.476540e-10  4.885316 1.000 0.062 4.767728e-06
GNG11     4.483201e-10  4.581688 0.929 0.000 6.148262e-06
HIST1H2AC 4.483201e-10  4.060724 0.929 0.000 6.148262e-06
SPARC     4.483201e-10  3.554182 0.929 0.000 6.148262e-06
GP9       4.483201e-10  3.553585 0.929 0.000 6.148262e-06
> head(av)
                      DC Platelet       diff
AL627309.1    0.00000000        0  0.0000000
AP006222.2    0.00000000        0  0.0000000
RP11-206L10.2 0.08462847        0 -0.6131753
RP11-206L10.9 0.00000000        0  0.0000000
LINC00115     0.00000000        0  0.0000000
NOC2L         0.32861136        0 -1.4553804

如果这样的肉眼去查看，就是粉丝提问的截图，当然是看不出来啊！会错误的以为，两个函数得到的结果差异很大！

我这里做一个小小的可视化，代码如下：

ids=intersect(rownames(diff),
              rownames(av))
comp=data.frame(FindMarkers=diff[ids,'avg_log2FC'],
           AverageExpression=av[ids,'diff'])
head(comp)
plot(comp)
head(comp)
library(ggpubr)
ggscatter(comp, x = "FindMarkers", y = "AverageExpression",
          color = "black", shape = 21, size = 3, # Points color, shape and size
          add = "reg.line",  # Add regressin line
          add.params = list(color = "blue", fill = "lightgray"), # Customize reg. line
          conf.int = TRUE, # Add confidence interval
          cor.coef = TRUE, # Add correlation coefficient. see ?stat_cor
          cor.coeff.args = list(method = "pearson",  label.sep = "\n")
)

结果如下：

 

可以看到，两个函数，FindMarkers与AverageExpression，其实计算得到 差异变化几乎是没有差异，这样的相关性已经是非常好了。（可能是粉丝对数据的一致性这个判断没有相应的统计学背景知识）

那么，有没有可能让两个函数 FindMarkers与AverageExpression得到的结果完全一致呢？

那就期待明天的讲解哈！

另外，如果你也想拿到前面的 pbmc_Platelet_DC.Rdata 文件，走这个代码，其实超级容易，就是pbmc3k数据集哦。

往期回顾

---

---

---

如果你对单细胞转录组研究感兴趣，但又不知道如何入门，也许你可以关注一下下面的课程

看完记得顺手点个“在看”哦！

生物 | 单细胞 | 转录组丨资料每天都精彩

长按扫码可关注