clusterProfiler也能做细胞类型鉴定？

在单细胞数据分析中，比较棘手的一件事莫过于细胞类型的鉴定了：虽然细胞可以用聚类算法聚成好几类，每个类型的细胞是什么呢?这就考验技术了。其实我们介绍过好几种细胞类型鉴定的方法了：

今天在读clusterProfiler 的文档的时候，突然发现，这个也可以用来做细胞类型的鉴定啊。我们拿每个cluster的差异基因（或者marker基因）来与cellMark数据库里的基因做富集分析不就可以看出某一亚群富集到那种细胞类型中了吗？！

首先请出我们熟悉的数据集pbmc，并计算Cluster 1 的差异基因：

library(Seurat)
pbmc<-readRDS("G:\\Desktop\\Desktop\\Novo周运来\\SingleCell\\scrna_tools\\pbmc3k_final.rds") # 已经分过群了的
marker<-FindMarkers(pbmc,group.by ="seurat_clusters",ident.1 = 1,logfc.threshold = 0.5, min.pct = 0.5)
head(marker)

               p_val avg_logFC pct.1 pct.2     p_val_adj
S100A9  0.000000e+00  3.860873 0.996 0.215  0.000000e+00
S100A8  0.000000e+00  3.796640 0.975 0.121  0.000000e+00
LGALS2  0.000000e+00  2.634294 0.908 0.059  0.000000e+00
FCN1    0.000000e+00  2.352693 0.952 0.151  0.000000e+00
CD14   2.856582e-294  1.951644 0.667 0.028 3.917516e-290
TYROBP 3.190467e-284  2.111879 0.994 0.265 4.375406e-280

大家看到这些基因不是ENTREZID ，所以我用clusterProfiler的bitr函数转化一下：

library("clusterProfiler")
library(org.Hs.eg.db)
gene.df <- bitr(rownames(marker), fromType = "SYMBOL", #fromType是指你的数据ID类型是属于哪一类的  ENTREZID
                toType = c("ENSEMBL", "ENTREZID"), #toType是指你要转换成哪种ID类型，可以写多种，也可以只写一种
                OrgDb =org.Hs.eg.db   )

head(gene.df)

  SYMBOL         ENSEMBL ENTREZID
1 S100A9 ENSG00000163220     6280
2 S100A8 ENSG00000143546     6279
3 LGALS2 ENSG00000100079     3957
4   FCN1 ENSG00000085265     2219
5   CD14 ENSG00000170458      929
6 TYROBP ENSG00000011600     7305

接下来我加载cellmarker的数据：

cell_markers <- vroom::vroom('http://bio-bigdata.hrbmu.edu.cn/CellMarker/download/Human_cell_markers.txt') %>%
  tidyr::unite("cellMarker", tissueType, cancerType, cellName, sep=", ") %>% 
  dplyr::select(cellMarker, geneID) %>%
  dplyr::mutate(geneID = strsplit(geneID, ', '))
cell_markers

   cellMarker                                             geneID   
   <chr>                                                  <list>   
 1 Kidney, Normal, Proximal tubular cell                  <chr [1]>
 2 Liver, Normal, Ito cell (hepatic stellate cell)        <chr [1]>
 3 Endometrium, Normal, Trophoblast cell                  <chr [1]>
 4 Germ, Normal, Primordial germ cell                     <chr [1]>
 5 Corneal epithelium, Normal, Epithelial cell            <chr [1]>
 6 Placenta, Normal, Cytotrophoblast                      <chr [1]>
 7 Periosteum, Normal, Periosteum-derived progenitor cell <chr [4]>
 8 Amniotic membrane, Normal, Amnion epithelial cell      <chr [2]>
 9 Primitive streak, Normal, Primitive streak cell        <chr [2]>
10 Adipose tissue, Normal, Stromal vascular fraction cell <chr [1]>
# ... with 2,858 more rows

最关键的一步来了，做富集：


y <- enricher(gene.df$ENTREZID, TERM2GENE=cell_markers, minGSSize=1)
DT::datatable(as.data.frame(y))

看到了吗？Cluster 1的差异基因与cellmarker的富集结果出来啦。我们可以根据我们的组织样本类型以及富集结果的指标来进一步判断该群的细胞类型。

我们看到展示的是Gene ID，如果我想在展示gene SYMBOL 该如何批量替换呢？

id <- gene.df$ENTREZID
names(id) <- gene.df$SYMBOL  # 用R语言的names函数构建类似字典的数据结构
id

  IL32      LTB      LTB      LTB      LTB      LTB      LTB      LTB      LTB     CD3D     IL7R     LDHB      CD2     CD3E  HLA-DRA  HLA-DRA  HLA-DRA  HLA-DRA 
  "9235"   "4050"   "4050"   "4050"   "4050"   "4050"   "4050"   "4050"   "4050"    "915"   "3575"   "3945"    "914"    "916"   "3122"   "3122"   "3122"   "3122" 
 HLA-DRA  HLA-DRA  HLA-DRA HLA-DRB1 HLA-DRB1 HLA-DRB1 HLA-DRB1 HLA-DRB1 HLA-DRB1     CD74 HLA-DPA1 HLA-DPA1 HLA-DPA1 HLA-DPA1 HLA-DPA1 HLA-DPA1 HLA-DPA1 HLA-DPA1 
  "3122"   "3122"   "3122"   "3123"   "3123"   "3123"   "3123"   "3123"   "3123"    "972"   "3113"   "3113"   "3113"   "3113"   "3113"   "3113"   "3113"   "3113" 
HLA-DPB1 HLA-DPB1 HLA-DPB1 HLA-DPB1 HLA-DPB1 HLA-DPB1 HLA-DPB1 HLA-DPB1     CYBA      FTL     CTSS     OAZ1      LYZ  GABARAP     PSAP     SAT1     FTH1     SRGN 
  "3115"   "3115"   "3115"   "3115"   "3115"   "3115"   "3115"   "3115"   "1535"   "2512"   "1520"   "4946"   "4069"  "11337"   "5660"   "6303"   "2495"   "5552" 
  S100A6    COTL1 
  "6277"  "23406" 

y@result-> res  # 以防数据被搞坏，新建一个

res$genesym <- unlist(lapply(y@result$geneID,  FUN =function(x){paste( unlist(lapply(unlist(str_split(x,"/")),FUN=function(x){names(id[which(id ==x)])})) , collapse = "/")} ))
DT::datatable(res)

似乎没有取unique，有转换多的情况，取一个吧。

id <- unique(gene.df$ENTREZID)
names(id) <- unique(gene.df$SYMBOL)

id
y@result-> res

res$genesym <- unlist(lapply(res$geneID,  FUN =function(x){paste( unlist(lapply(unlist(str_split(x,"/")),FUN=function(x){names(id[which(id ==x)])})) , collapse = "/")} ))
DT::datatable(res)

当然，细胞类型鉴定需要结合具体的生物学意义以及大量的背景知识，clusterProfiler只是一个工具，其实关于富集我们可以做很多工作，想象力可以改变世界。

clusterProfiler-book