ixxmu / mp_duty

抓取网络文章到github issues保存
https://archives.duty-machine.now.sh/
116 stars 30 forks source link

用clusterProfiler一行代码对单细胞的marker genes做功能富集分析 #3001

Closed ixxmu closed 1 year ago

ixxmu commented 1 year ago

https://mp.weixin.qq.com/s/nvOXM2hpZq3UI_Npjhk8HQ

ixxmu commented 1 year ago

用clusterProfiler一行代码对单细胞的marker genes做功能富集分析 by YuLabSMU

曾老师发来他写的《单细胞各个亚群特异性高表达基因的数据库注释(包括GO,KEGG,ReactomePA)》,我看了一下,还能更简单一些。

我本来是想着,我是不是可以写点小代码衔接一下,结果啥也不用写,已然是衔接得很好了。

首先是数据

library(SeuratData)
data("pbmc3k")
sce <- pbmc3k.final

再者是marker基因的鉴定

library(Seurat)
sce.markers <- FindAllMarkers(object = sce, only.pos = TRUE
                            min.pct = 0.25,
                            thresh.use = 0.25)

总要过滤一下吧

library(dplyr)
markers <- sce.markers |> group_by(cluster) |>
    filter(p_val_adj < 0.001) |>
    ungroup()

clusterProfiler上场

首先假如我们需要换个ID,只需要用bitr函数。然后对着data.frame的结果,来个full_join即可,并不需要list倒腾过来,又倒腾过去。

library(clusterProfiler)
gid <- bitr(unique(markers$gene), 'SYMBOL''ENTREZID', OrgDb= 'org.Hs.eg.db')
markers <- full_join(markers, gid, by=c('gene' = 'SYMBOL'))

然后这样的data.frame就可以直接分析啦。如果你还不了解,那么你应该看我们clusterProfiler 4.0的文章。

x = compareCluster(ENTREZID ~ cluster, data = markers, fun='enrichKEGG')

然后就是出图了:

dotplot(x, label_format=40) + theme(axis.text.x = element_text(angle=45, hjust=1)) 


假如我们是用COSG

那么还是一样,先找marker:

library(COSG)
marker_cosg <- cosg(
      sce,
      groups='all',
      assay='RNA',
      slot='data',
      mu=1,
      n_genes_user=100)

输出的list的第一个元素是个这样子的data.frame:

> head(marker_cosg[[1]])
  Naive CD4 T Memory CD4 T CD14+ Mono         B CD8 T  FCGR3A+ Mono     NK
1        CCR7      TNFRSF4     S100A8     CD79A  GZMK        CDKN1C   GZMB
2        FHIT         AQP3     LGALS2     MS4A1  CCL5          HES4   GNLY
3        LEF1         IL7R     S100A9     TCL1A  CD8A         MS4A7  SPON2
4        LDHB          CD2       CD14     CD79B KLRG1           CKB FGFBP2
5   PRKCQ-AS1        CRIP2     MS4A6A LINC00926  LAG3 RP11-290F20.3 AKR1C3
6     PIK3IP1        TRADD       FCN1    VPREB3  NKG7        MS4A4A   PRF1
        DC   Platelet
1   FCER1A        GP9
2 SERPINF1        PF4
3  CLEC10A AP001189.4
4     ENHO     ITGA2B
5    CLIC2      GNG11
6   CLEC4C     TMEM40

每一个column是一个cluster,有多少个clusters就有多少个columns,然后因为大家都取前100个基因,所以等长,形成一个data.frame。大家都知道compareCluster的默认输入是一个named list,很多人在学R的时候,也了解到data.frame其实是个等长的list,这时候能够变通的人就很少了,它是个list,而且是个named list,名字就是column names啊。所以这个结果,直接就可以是clusterProfiler的输入了。

于是一行代码搞定。

y = compareCluster(marker_cosg[[1]], fun='enrichGO'
          OrgDb = 'org.Hs.eg.db', keyType = 'SYMBOL', ont="MF")

翠花上酸菜图片:

dotplot(y, label_format=60) + theme(axis.text.x = element_text(angle=45, hjust=1))

clusterProfiler以其易用性著称,许多人说这是他们学生信所学的第一个R包。但显而易见,这个包,比大家想象中的还容易用的。

ixxmu commented 11 months ago

每一个column是一个cluster,有多少个clusters就有多少个columns,然后因为大家都取前100个基因,所以等长,形成一个data.frame。大家都知道compareCluster的默认输入是一个named list,很多人在学R的时候,也了解到data.frame其实是个等长的list,这时候能够变通的人就很少了,它是个list,而且是个named list,名字就是column names啊。所以这个结果,直接就可以是clusterProfiler的输入了。

data.frame 就是一个等长的list,column names就是 column names, 大佬的理解就是不一样啊!