在單細(xì)胞數(shù)據(jù)分析中，比較棘手的一件事莫過于細(xì)胞類型的鑒定了：雖然細(xì)胞可以用聚類算法聚成好幾類，每個(gè)類型的細(xì)胞是什么呢?這就考驗(yàn)技術(shù)了。其實(shí)我們介紹過好幾種細(xì)胞類型鑒定的方法了：
- 單細(xì)胞數(shù)據(jù)挖掘||Garnett :基于監(jiān)督學(xué)習(xí)的細(xì)胞注釋算法
- 10X 單細(xì)胞細(xì)胞類型手動(dòng)定義||Identifying Cell Types with Loupe Cell Browser
- CHETAH: sc-RNA-seq 細(xì)胞類型定義工具
- SingleR || 單細(xì)胞細(xì)胞類型定義工具
- sc-RNA-seq Key issues guide ||重新理解細(xì)胞類型
- celaref ||單細(xì)胞細(xì)胞類型定義工具
- single cell marker 基因數(shù)據(jù)庫
- SuperCT: sc-RNA-seq 單細(xì)胞細(xì)胞類型定義在線工具
- cellassign：基于Tensorfiow 框架的sc-RNA細(xì)胞類型鑒定方法
- sc-Review||Best Practices:聚類與細(xì)胞鑒定
今天在讀clusterProfiler 的文檔的時(shí)候，突然發(fā)現(xiàn)，這個(gè)也可以用來做細(xì)胞類型的鑒定啊。我們拿每個(gè)cluster的差異基因（或者marker基因）來與cellMark數(shù)據(jù)庫里的基因做富集分析不就可以看出某一亞群富集到那種細(xì)胞類型中了嗎？！
首先請(qǐng)出我們熟悉的數(shù)據(jù)集pbmc，并計(jì)算Cluster 1 的差異基因：
library(Seurat)
pbmc<-readRDS("G:\\Desktop\\Desktop\\Novo周運(yùn)來\\SingleCell\\scrna_tools\\pbmc3k_final.rds") # 已經(jīng)分過群了的
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函數(shù)轉(zhuǎn)化一下：
library("clusterProfiler")
library(org.Hs.eg.db)
gene.df <- bitr(rownames(marker), fromType = "SYMBOL", #fromType是指你的數(shù)據(jù)ID類型是屬于哪一類的 ENTREZID
toType = c("ENSEMBL", "ENTREZID"), #toType是指你要轉(zhuǎn)換成哪種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的數(shù)據(jù)：
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
最關(guān)鍵的一步來了，做富集：
y <- enricher(gene.df$ENTREZID, TERM2GENE=cell_markers, minGSSize=1)
DT::datatable(as.data.frame(y))
看到了嗎？Cluster 1的差異基因與cellmarker的富集結(jié)果出來啦。我們可以根據(jù)我們的組織樣本類型以及富集結(jié)果的指標(biāo)來進(jìn)一步判斷該群的細(xì)胞類型。
我們看到展示的是Gene ID，如果我想在展示gene SYMBOL 該如何批量替換呢？
id <- gene.df$ENTREZID
names(id) <- gene.df$SYMBOL # 用R語言的names函數(shù)構(gòu)建類似字典的數(shù)據(jù)結(jié)構(gòu)
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 # 以防數(shù)據(jù)被搞壞，新建一個(gè)
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，有轉(zhuǎn)換多的情況，取一個(gè)吧。
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)
當(dāng)然，細(xì)胞類型鑒定需要結(jié)合具體的生物學(xué)意義以及大量的背景知識(shí)，clusterProfiler只是一個(gè)工具，其實(shí)關(guān)于富集我們可以做很多工作，想象力可以改變世界。
