Seurat_StepVI_CellType
此文章是通过学习瑞典国家生物信息学基础设施(NBIS) 所开放的单细胞分析教程加上网上所查找的资料，自身的理解所形成的，可能会有不足之处。该部分是对不同集群的细胞类型进行识别。
参考来源：https://github.com/NBISweden/workshop-scRNAseq/blob/master/labs/seurat/seurat_06_celltype.Rmd
感兴趣的话可以阅读原文。
细胞类型预测
可以预测每个细胞的细胞的细胞类型，也可以对整个细胞集群上进行预测，所有的预测都是基于与其他的数据集，单细胞数据集或分类好的RNA-seq数据集之间的相似性，或是自身所储备的有关细胞类型的marker基因相关知识进行预测。
在此教程中，从covid数据集中选择一个ctrl_13样本进行细胞类型的预测。
有一些方法会根据与已知的细胞类型之间的相似性进行细胞类型的预测，即使这个细胞类型没有参考文献的支撑（问题来了，既然没有参考文献的支撑，那怎么判断自己所分析的这个细胞类型是对的呢？），而其他包含令人混淆的基因类型的细胞集群就没有办法进行预测，因为他们与已知细胞类型marker之间的相似性较低。
目前，有很多方法可以对细胞类型进行预测，在这只介绍一些方法，剩下的感兴趣可以自己查阅资料。
加载数据
# devtools::install_github("powellgenomicslab/scPred")
suppressPackageStartupMessages({
library(Seurat)
library(dplyr)
library(cowplot)
library(ggplot2)
library(pheatmap)
library(rafalib)
library(scPred)
})
# 加载数据集并选择ctrl_13样本
alldata <- readRDS("data/results/covid_qc_dr_int_cl.rds")
ctrl = alldata[, alldata$orig.ident == "ctrl_13"]
# 将默认基因表达矩阵设置为RNA，并删除CCA（批次整合矩阵）
ctrl@active.assay = "RNA"
ctrl[["CCA"]] = NULL
ctrl
# An object of class Seurat
# 18121 features across 1135 samples within 1 assay
# Active assay: RNA (18121 features, 0 variable features)
# 4 dimensional reductions calculated: umap, tsne, harmony, umap_harmony
参考数据集注释
预测前处理
- 参考数据集处理
# 参考数据集，可以使用标准化、可变基因选择、缩放和降维运行正常分析的所有步骤
reference <- scPred::pbmc_1
reference
# An object of class Seurat
# 32838 features across 3500 samples within 1 assay
# Active assay: RNA (32838 features, 0 variable features)
reference <- reference %>%
NormalizeData() %>%
FindVariableFeatures() %>%
ScaleData() %>%
RunPCA(verbose = F) %>%
RunUMAP(dims = 1:30)
# 绘图查看细胞类型
DimPlot(reference, group.by = "cell_type", label = TRUE, repel = TRUE) + NoAxes()
- ctrl_13样本处理
# 选择分辨率为0.3
ctrl <- SetIdent(ctrl, value = "CCA_snn_res.0.5")
ctrl <- ctrl %>%
NormalizeData() %>%
FindVariableFeatures() %>%
ScaleData() %>%
RunPCA(verbose = F) %>%
RunUMAP(dims = 1:30)
DimPlot(ctrl, label = TRUE, repel = TRUE) + NoAxes()
Seurat预测
首先，使用整合不同批次数据中类似的方法，使用自带函数FindTransferAnchors找到两个数据集不同细胞类型合并的锚点，默认方法是基于PCA空间pcaproject，将预测数据集投影到参考数据集的PCA空间中，然后将参考数据集对应细胞类型的标签就传递给了预测数据集，实现了细胞类型的注释。
transfer.anchors <- FindTransferAnchors(reference = reference, query = ctrl, dims = 1:30)
predictions <- TransferData(anchorset = transfer.anchors, refdata = reference$cell_type,
dims = 1:30)
ctrl <- AddMetaData(object = ctrl, metadata = predictions)
DimPlot(ctrl, group.by = "predicted.id", label = T, repel = T) + NoAxes()
也可以绘制出不同细胞集群中，细胞类型的分布。
ggplot(ctrl@meta.data, aes(x = CCA_snn_res.0.5, fill = predicted.id)) + geom_bar() +
theme_classic()
scPred预测
scPred 将基于所有主成分训练分类器,是一个机器学习的方法。首先，getFeatureSpace 将创建一个 scPred 对象，存储在 @misc 插槽中，从中提取最能区分不同细胞类型的 PC。然后 trainModel 将对每种细胞类型进行训练。
# 获取每个细胞类型的代表性PC
reference <- getFeatureSpace(reference, "cell_type")
# 针对这些代表性基因进行模型训练
reference <- trainModel(reference)
# 使用get_scpred可以查看每种细胞类型所使用的PC数，ROC值与灵敏度/特异性
get_scpred(reference)
## 'scPred' object
## ✔ Prediction variable = cell_type
## ✔ Discriminant features per cell type
## ✔ Training model(s)
## Summary
##
## |Cell type | n| Features|Method | ROC| Sens| Spec|
## |:-----------|----:|--------:|:---------|-----:|-----:|-----:|
## |B cell | 280| 50|svmRadial | 1.000| 0.964| 1.000|
## |CD4 T cell | 1620| 50|svmRadial | 0.997| 0.971| 0.975|
## |CD8 T cell | 945| 50|svmRadial | 0.985| 0.902| 0.978|
## |cDC | 26| 50|svmRadial | 0.995| 0.547| 1.000|
## |cMono | 212| 50|svmRadial | 0.994| 0.958| 0.970|
## |ncMono | 79| 50|svmRadial | 0.998| 0.582| 1.000|
## |NK cell | 312| 50|svmRadial | 0.999| 0.936| 0.996|
## |pDC | 20| 50|svmRadial | 1.000| 0.700| 1.000|
## |Plasma cell | 6| 50|svmRadial | 1.000| 0.800| 1.000|
可以通过更改参数和测试不同类型的模型来优化每个数据集的参数，更多信息请参见：https://powellgenomicslab.github.io/scPred/articles/introduction.html。但就目前而言，继续使用此模型。
使用scPredict函数预测样本细胞的类型。
ctrl <- scPredict(ctrl, reference)
DimPlot(ctrl, group.by = "scpred_prediction", label = T, repel = T) + NoAxes()
ggplot(ctrl@meta.data, aes(x = CCA_snn_res.0.5, fill = scpred_prediction)) + geom_bar() +
theme_classic()
结果比较
crossTab(ctrl, "predicted.id", "scpred_prediction")
# B cell CD4 T cell CD8 T cell cDC cMono ncMono NK cell pDC Plasma cell
# B cell 103 1 2 0 1 0 0 0 0
# CD4 T cell 0 200 0 0 0 0 0 0 0
# CD8 T cell 0 7 228 0 1 0 3 0 0
# cDC 0 0 0 13 8 0 0 0 0
# cMono 0 6 2 0 197 5 0 0 0
# ncMono 0 0 0 0 8 101 0 0 0
# NK cell 0 0 15 0 0 0 151 0 0
# pDC 0 0 0 0 0 0 0 1 0
# Plasma cell 0 1 0 0 0 0 0 0 2
# unassigned 0 10 10 0 57 1 1 0 0
可以看出两种方法对于细胞类型的注释，结果是相似的。
CellMarker注释
除了上述的基于参考数据集来进行细胞注释之外，还可以使用在不同细胞类型之间分布具有差异的marker基因进行注释，类似于我们在差异基因分析时所使用基因富集分析类似，现在有一些类似于通路信息的细胞类型的基因集信息，如MSigDB中的CellMarker，PangiaoDN或celltype基因集，除了这些现有的基因集之外，还可以使用参考数据集中分析得到的marker基因对目标数据进行细胞类型的注释。
参考数据集差异基因集注释
- 目标数据集差异基因分析
# 首先对目标基因集进行差异基因的分析
alldata <- SetIdent(alldata, value = "CCA_snn_res.0.5")
DGE_table <- FindAllMarkers(alldata, logfc.threshold = 0, test.use = "wilcox", min.pct = 0.1,
min.diff.pct = 0, only.pos = TRUE, max.cells.per.ident = 20, return.thresh = 1,
assay = "RNA")
# 将每个集群中的差异基因单独成表，最后合并成一个列表
DGE_list <- split(DGE_table, DGE_table$cluster)
unlist(lapply(DGE_list, nrow))
# 0 1 2 3 4 5 6 7 8 9
# 3167 3318 2522 4100 2091 3799 2473 3308 2233 2502
- 参考数据集差异基因分析
# 计算参考数据集的差异基因分布
reference <- SetIdent(reference, value = "cell_type")
reference_markers <- FindAllMarkers(reference, min.pct = 0.1, min.diff.pct = 0.2,
only.pos = T, max.cells.per.ident = 20, return.thresh = 1)
# 在参考数据集不同细胞类型的前100个marker基因中选择50个差异倍数变化最大的
reference_markers <- reference_markers[order(reference_markers$avg_log2FC, decreasing = T),
]
reference_markers %>%
group_by(cluster) %>%
top_n(-100, p_val) %>%
top_n(50, avg_log2FC) -> top50_cell_selection
# 将这些marker基因转化为列表
ref_list = split(top50_cell_selection$gene, top50_cell_selection$cluster)
unlist(lapply(ref_list, length))
# CD8 T cell CD4 T cell cMono B cell NK cell pDC ncMono
# 30 15 50 50 50 50 50
# cDC Plasma cell
# 50 50
现在就可以对目标数据集中的差异基因集进行GESA分析，并检查目标数据集的差异基因在参考数据集的细胞类型代表基因中的富集程度。
suppressPackageStartupMessages(library(fgsea))
# 在目标数据集的每个集群上运行GESA
res <- lapply(DGE_list, function(x) {
gene_rank <- setNames(x$avg_log2FC, x$gene)
fgseaRes <- fgsea(pathways = ref_list, stats = gene_rank, nperm = 10000)
return(fgseaRes)
})
names(res) <- names(DGE_list)
# 可以基于 ES 富集分数, NES 标准富集分数 or pvalue p值过滤某些信息
res <- lapply(res, function(x) {
x[x$pval < 0.1, ]
})
res <- lapply(res, function(x) {
x[x$size > 2, ]
})
res <- lapply(res, function(x) {
x[order(x$NES, decreasing = T), ]
})
head(res,10)
# $`0`
# pathway pval padj ES NES nMoreExtreme size
# 1: cMono 0.0000999900 0.0002999700 0.9578612 2.087484 0 48
# 2: ncMono 0.0000999900 0.0002999700 0.8406934 1.826483 0 45
# 3: cDC 0.0000999900 0.0002999700 0.8246721 1.782681 0 41
# 4: pDC 0.0004011231 0.0009025271 0.7760109 1.598105 3 20
# 5: NK cell 0.0154702970 0.0278465347 0.7605442 1.477318 149 11
# 6: B cell 0.0316654826 0.0474982239 0.7043057 1.403353 311 14
# leadingEdge
# 1: S100A8,S100A9,LYZ,S100A12,VCAN,FCN1,...
# 2: CTSS,TYMP,CST3,S100A11,AIF1,SERPINA1,...
# 3: LYZ,GRN,TYMP,CST3,AIF1,SPI1,...
# 4: GRN,MS4A6A,CST3,MPEG1,CTSB,TGFBI,...
# 5: TYROBP,FCER1G,SRGN,CCL3,CD63,MYO1F,...
# 6: NCF1,LY86,MARCH1,PHACTR1,HLA-DRB5,POU2F2,...
接下来就可以根据预测的标签来命名细胞类型了，但是，如果某些集群的p值都不好，那这些集群的细胞类型可靠度大大降低，有可能是因为所使用的参考数据集没有办法涵盖在目标数据集中的细胞类型，所以只能预测最相似的细胞，此外，也有一些细胞类型的marker基因之间会比较相似，那这样所得到的p值就会比较低。
new.cluster.ids <- unlist(lapply(res, function(x) {
as.data.frame(x)[1, 1]
}))
alldata$ref_gsea <- new.cluster.ids[as.character(alldata@active.ident)]
cowplot::plot_grid(ncol = 2, DimPlot(alldata, label = T, group.by = "CCA_snn_res.0.5") +
NoAxes(), DimPlot(alldata, label = T, group.by = "ref_gsea") + NoAxes())
已知细胞类型基因集
就如上所述，除了使用参考数据集进行差异分析从而得到marker基因外，现在已经有了对特定细胞类型marker基因整合好的数据库，如Cellmarker，接下来就可以使用cellmarker数据库进行后续注释。
使用数据库中的marker基因比使用参考数据集进行差异分析得到的marker基因来说，基因的种类更全，细胞类型更全，但是会收到数据库的局限性，不能及时更新，所以需要我们自己时刻关注最新的单细胞数据，及时更新自身的单细胞参考基因集。
# 下载marker基因集
if (!dir.exists("data/CellMarker_list/")) {
dir.create("data/CellMarker_list")
download.file(url = "http://yikedaxue.slwshop.cn/Cell_marker_Human.xlsx",
destfile = "./data/CellMarker_list/Human_cell_markers.xlsx")
download.file(url = "http://yikedaxue.slwshop.cn/Cell_marker_Mouse.xlsx",
destfile = "./data/CellMarker_list/Mouse_cell_markers.xlsx")
}
# 加载cellmarker基因集
markers <- readxl::read_xlsx("data/CellMarker_list/Human_cell_markers.xlsx",sheet = 1)
markers <- markers[markers$Species == "Human", ]
markers <- markers[markers$`Cancer type` == "Normal", ]
# 也可以按组织过滤（能够减少计算时间并且能够进行组织特异性的分类）
# sort(unique(markers$`Tissue type`))
# grep('blood',unique(markers$`Tissue type`),value = T)
# markers <- markers[markers$`Tissue type` %in% c('Blood','Venous blood', 'Serum','Plasma','Spleen','Bone marrow','Lymph node'), ]
# 移除掉陌生的字符串
celltype_list <- lapply(unique(markers$`Cell name`), function(x) {
x <- paste(markers$Marker[markers$`Cell name` == x], sep = ",")
x <- gsub("[[]|[]]| |-", ",", x)
x <- unlist(strsplit(x, split = ","))
x <- unique(x[!x %in% c("", "NA", "family")])
x <- casefold(x, upper = T)
})
names(celltype_list) <- unique(markers$`Cell name`)
# celltype_list <- lapply(celltype_list , function(x) {x[1:min(length(x),50)]}
# )
celltype_list <- celltype_list[unlist(lapply(celltype_list, length)) < 100]
celltype_list <- celltype_list[unlist(lapply(celltype_list, length)) > 5]
- GESA分析-细胞类型判别
# 对每个集群按照该基因集进行GESA分析
res <- lapply(DGE_list, function(x) {
gene_rank <- setNames(x$avg_log2FC, x$gene)
fgseaRes <- fgsea(pathways = celltype_list, stats = gene_rank, nperm = 10000)
return(fgseaRes)
})
names(res) <- names(DGE_list)
# 可以基于 ES 富集分数, NES 标准富集分数 or pvalue p值过滤某些信息
res <- lapply(res, function(x) {
x[x$pval < 0.01, ]
})
res <- lapply(res, function(x) {
x[x$size > 5, ]
})
res <- lapply(res, function(x) {
x[order(x$NES, decreasing = T), ]
})
# 展示每个集群前三个p值最大的细胞类型.
lapply(res, head, 3)
# $`0`
# pathway pval padj ES NES nMoreExtreme size
# 1: Eosinophil 0.0001039285 0.00795053 0.9039267 1.740575 0 10
# 2: Classical monocyte 0.0001000500 0.00795053 0.8228162 1.737254 0 28
# 3: CD1C+_B dendritic cell 0.0000999900 0.00795053 0.7911922 1.721561 0 48
# leadingEdge
# 1: S100A8,S100A9,LYZ,CD14,RETN,CD33,...
# 2: S100A8,S100A9,LYZ,S100A12,VCAN,FCN1,...
# 3: S100A8,S100A9,LYZ,S100A12,VCAN,FCN1,...
# 细胞类型绘图
new.cluster.ids <- unlist(lapply(res, function(x) {
as.data.frame(x)[1, 1]
}))
alldata$cellmarker_gsea <- new.cluster.ids[as.character(alldata@active.ident)]
# 结果比较
cowplot::plot_grid(ncol = 2, DimPlot(alldata, label = T, group.by = "ref_gsea") +
NoAxes(), DimPlot(alldata, label = T, group.by = "cellmarker_gsea") + NoAxes())
保存数据
saveRDS(ctrl, "data/results/ctrl13_qc_dr_int_cl_celltype.rds")
1595
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
