Cell-ID发表在NBT上面：Genesignature extraction and cell identity recognition at the single-cell levelwith Cell-ID。

最近我看了很多细胞类型注释的思路。在谈到自动注释时，又基于marker gene的方法，例如scCATCH，SCCA等前面都介绍过了。还有基于相关性也就是elation分析的，例如SingleR等。还有基于机器学习的，例如Garnett，ACTINN，srClassify等。类似的软件还有很多，原理上可能五花八门。很多软件在做注释时会基于某个细胞和参考细胞表达谱的相似性，今天看到的这个软件使用的是另一种思路：将某个细胞的特征基因集与表征细胞类型的参考基因集做富集分析，当在某个细胞类型的基因集上显著富集时，就定义为这个细胞类型（和GO富集很像）。所以这种方法的关键点是：找单个细胞的特征基因集，以及找细胞类型的参考基因集。

==安装====

devtools::install_github("RausellLab/CelliD",ref = "legacy")

==例子测试==

library(CellID)

library(Seurat)

library(tidyverse)

#获取原始表达矩阵

BaronMatrix   <-readRDS(url("https://storage.googleapis.com/cellid-cbl/BaronMatrix.rds"))

//也可以下载其它几组数据

BaronMetaData <- readRDS("BaronMetaData.rds")

SegerMatrix   <- readRDS("SegerstolpeMatrix.rds")

SegerMetaData <- readRDS("SegerstolpeMetaData2.rds")

#仅考虑编码蛋白质的基因

data("HgProteinCodingGenes")

BaronMatrixProt <- BaronMatrix[rownames(BaronMatrix) %in% HgProteinCodingGenes,]

SegerMatrixProt <- SegerMatrix[rownames(SegerMatrix) %in% HgProteinCodingGenes,]

#这几步类似Seurat的标准流程

Baron <- CreateSeuratObject(counts = BaronMatrixProt, project = "Baron", min.cells = 5, meta.data = BaronMetaData)

Seger <- CreateSeuratObject(counts = SegerMatrixProt, project = "Segerstolpe", min.cells = 5, meta.data = SegerMetaData)

Baron <- NormalizeData(Baron)

Baron <- ScaleData(Baron, features = rownames(Baron))

Baron <- RunMCA(Baron) //#该软件的主要分析函数，将细胞和基因同时降维到一个空间，离细胞近的基因被定义为细胞的特征基因

DimPlotMC(Baron, reduction = "mca", group.by = "cell.type", features = c("CTRL", "INS", "MYZAP", "CDH11"), as.text = TRUE) + ggtitle("MCA with some key gene markers") 

//还可以将MCA的结果和PCA，UMAP以及TSNE做对比

Baron <- RunPCA(Baron, features = rownames(Baron))

Baron <- RunUMAP(Baron, dims = 1:30)

Baron <- RunTSNE(Baron, dims = 1:30)

PCA  <- DimPlot(Baron, reduction = "pca", group.by = "cell.type")  + ggtitle("PCA") + theme(legend.text = element_text(size =10), aspect.ratio = 1)

tSNE <- DimPlot(Baron, reduction = "tsne", group.by = "cell.type")+ ggtitle("tSNE") + theme(legend.text = element_text(size =10), aspect.ratio = 1)

UMAP <- DimPlot(Baron, reduction = "umap", group.by = "cell.type") + ggtitle("UMAP") + theme(legend.text = element_text(size =10), aspect.ratio = 1)

MCA <- DimPlot(Baron, reduction = "mca", group.by = "cell.type")  + ggtitle("MCA") + theme(legend.text = element_text(size =10), aspect.ratio = 1)

ggarrange(PCA, MCA, common.legend = TRUE, legend = "top")

ggarrange(tSNE, UMAP, common.legend = TRUE, legend = "top")

同样的，cell-ID也可以利用预先知道的marker集合。

下载参考基因集

panglao <-read_tsv("https://panglaodb.se/markers/PanglaoDB_markers_27_Mar_2020.tsv.gz")

#根据器官过滤，示例数据就是胰腺的

panglao_pancreas <- panglao %>%filter(organ == "Pancreas")

#物种包含人

panglao_pancreas <- panglao_pancreas%>%  filter(str_detect(species,"Hs"))

#下面两步会得到一个列表，列表的每一个元素是由基因集构成的字符向量，且每个元素被命名为对应的细胞类型的名字

panglao_pancreas <- panglao_pancreas%>% 

 group_by(`cell type`) %>% 

 summarise(geneset = list(`official gene symbol`))

pancreas_gs <-setNames(panglao_pancreas$geneset, panglao_pancreas$`cell type`)

得到的相当于每个cell中的marker gene list

#富集分析，用的超几何检验

HGT_pancreas_gs <- RunCellHGT(Baron,pathways = pancreas_gs, dims = 1:50, n.features = 200) #每个细胞选200个特征基因

富集分析之后会返回一个矩阵，矩阵中的每一个值表示p值在校正之后，求-log10，即-log10 corrected p-value，当然这个值越大越好。

pancreas_gs_prediction <-rownames(HGT_pancreas_gs)[apply(HGT_pancreas_gs, 2, which.max)]

#矩阵的每一列依次：判断是否是最大值，得到一列布尔值，结合矩阵的行名会返回行名中的一个元素（也就是最有可能的细胞类型）。

#所有列运行完了之后会得到所有细胞最可能的注释

pancreas_gs_prediction_signif <-ifelse(apply(HGT_pancreas_gs, 2, max)>2, yes = pancreas_gs_prediction,"unassigned")

#如果`-log10 correctedp-value`的值小于等于2，则认为不显著，注释更正为"unassigned"

Baron$pancreas_gs_prediction <-pancreas_gs_prediction_signif

所以到此，细胞类型的预测就结束了。

color <- c("#F8766D", "#E18A00", "#BE9C00", "#8CAB00", "#24B700", "#00BE70", "#00C1AB", "#00BBDA", "#00ACFC", "#8B93FF", "#D575FE", "#F962DD", "#FF65AC", "grey")

ggcolor <- setNames(color,c(sort(unique(Baron$cell.type)), "unassigned"))

OriginalPlot <- DimPlot(Baron, reduction = "tsne", group.by = "cell.type") + 

  scale_color_manual(values = ggcolor) + 

  theme(legend.text = element_text(size =10), aspect.ratio = 1)

Predplot1 <- DimPlot(Baron, reduction = "tsne", group.by = "pancreas_gs_prediction") + 

  scale_color_manual(values = ggcolor) + 

  theme(legend.text = element_text(size =10), aspect.ratio = 1)

ggarrange(OriginalPlot, Predplot1, legend = "top",common.legend = TRUE)

//可以把原先注释的和现在注释的做个对比，可以看出准确率还是可以的

//还可以试下用库中所有的marker，而不是挑选的Pancreas子集。

HGT_all_gs <- RunCellHGT(Baron, pathways = all_gs, dims = 1:50)

all_gs_prediction <- rownames(HGT_all_gs)[apply(HGT_all_gs, 2, which.max)]

all_gs_prediction_signif <- ifelse(apply(HGT_all_gs, 2, max)>2, yes = all_gs_prediction, "unassigned")

Baron$all_gs_prediction <- ifelse(all_gs_prediction_signif %in% c(names(pancreas_gs), "Schwann cells", "Endothelial cells", "Macrophages", "Mast cells", "T cells","Fibroblasts", "unassigned"), all_gs_prediction_signif,"other")

color <- c("#F8766D", "#E18A00", "#BE9C00", "#8CAB00", "#24B700", "#00BE70", "#00C1AB", "#00BBDA", "#00ACFC", "#8B93FF", "#D575FE", "#F962DD", "#FF65AC", "#D575FE", "#F962DD", "grey", "black")

ggcolor <- setNames(color,c(sort(unique(Baron$cell.type)), "Fibroblasts", "Schwann cells", "unassigned", "other"))

Baron$pancreas_gs_prediction <- factor(Baron$pancreas_gs_prediction,c(sort(unique(Baron$cell.type)), "Fibroblasts", "Schwann cells", "unassigned", "other"))

PredPlot2 <- DimPlot(Baron, group.by = "all_gs_prediction", reduction = "tsne") + 

             scale_color_manual(values = ggcolor, drop = FALSE) + 

             theme(legend.text = element_text(size = 10), aspect.ratio = 1)

ggarrange(OriginalPlot, PredPlot2, legend = "top", common.legend = TRUE)

//效果还是不如指定organ的好

做过GO或者pathway的人能体会的，参考GO库和pathway库的重要性。因此，可以看出这个方法非常依赖于PanglaoDB的准确性。但是对于不知道类型的cluster就比较方便，可以像GO一样推测可能的细胞类型。

本文使用 文章同步助手 同步