关键词
- FindTransferAnchors函数
- TransferData函数
- 细胞类型注释
- 映射数据集/Transfer 数据集
适用背景
细胞注释的标准一直备受争议，就目前来说仍没有一种统一的工具能实现让大多数人都信服的自动注释工具，主要还是靠人工根据不同细胞类型的经典marker基因来定义特定细胞类型。当然，这些软件注释出来的结果也可以当做一种辅助参考，可以同时结合软件和经验来更好注释细胞类型。
而Seurat作为一个专业的单细胞数据分析工具包，也推出了类似细胞类型注释的流程。与大多数细胞注释工具类似，其原理也是基于参考数据集的注释结果，通过找到待查数据集与参考数据集的锚点，然后给待查数据集的每一个细胞打上预测的细胞类型标签，然后可以根据最终结果进行人工校正。如果根据一下比较经典的marker基因也无法确定具体的细胞类型，则可以考虑通过transfer来进行鉴定。另外空间组数据与单细胞数据也可以直接通过transfer进行映射。
主函数
因为主函数有用到我用的比较多的调色板cor和快速可视化聚类的函数dimplot_2groups，需要先运行以下代码：
library(Seurat)
library(dplyr)
library(RColorBrewer)
cor<-c(brewer.pal(8, "Dark2"),brewer.pal(12, "Paired"),brewer.pal(8, "Set2"),brewer.pal(9, "Set1")）
dimplot_2groups <- function(obj,
group1="seurat_clusters", group2="orig.ident",
split.by=NULL,
wid=1800,hei=1000,
raster=NULL) {
if (is.null(raster)) {
png(paste0(group1,"_",group2,".png"),width = 1800, height = 1000)
print(DimPlot(obj, reduction = "umap", label = T, label.size = 5,repel =T)+scale_color_manual(values=cor) -
(DimPlot(obj, reduction = "umap", label = T, group.by = group2, label.size = 5,repel =T)+scale_color_manual(values=cor)))
dev.off()
if (!is.null(split.by)) {
png(paste0(split.by, "_split.png"),width = wid, height = hei)
print(DimPlot(obj, reduction = "umap", label = T, split.by = split.by, label.size = 5,repel =T)+scale_color_manual(values=cor))
dev.off()
}
}else{
png(paste0(group1,"_",group2,".png"),width = 1800, height = 1000)
print(DimPlot(obj, reduction = "umap", label = T, label.size = 5,repel =T,raster=FALSE)+scale_color_manual(values=cor) -
(DimPlot(obj, reduction = "umap", label = T, group.by = group2, label.size = 5,repel =T,raster=FALSE)+scale_color_manual(values=cor)))
dev.off()
if (!is.null(split.by)) {
png(paste0(split.by, "_split.png"),width = wid, height = hei)
print(DimPlot(obj, reduction = "umap", label = T, split.by = split.by, label.size = 5,repel =T,raster=FALSE)+scale_color_manual(values=cor))
dev.off()
}
}
}
我把流程封装在一个函数transfer_cell里，主要用到两个Seurat的函数FindTransferAnchors和TransferData，然后进行可视化。以下是参数介绍：
- ref 参考数据集，Seurat对象
- que 询问数据集，Seurat对象
- assay1 参考数据集的assay，默认为RNA
- assay2 询问数据集的assay，默认为RNA
- sct 是否对数据集进行SCTransform处理，默认为NULL，不作处理
- group1参考数据集进行映射的分组列名，默认为seurat_clusters
- group2参考数据集进行可视化的分组列名，默认为seurat_clusters，可不指定
- features FindTransferAnchors函数使用的基因集，默认为NULL，使用所有基因
- nor 是否对数据集使用NormalizeData，默认为NULL，会进行标准化处理
- nor_method 如果既没有SCTransform处理也没有NormalizeData，则需指定标准化方法，默认为LogNormalize
transfer_cell <- function(ref,que,assay1="RNA",assay2="RNA",sct=NULL,
group1="seurat_clusters",group2="seurat_clusters",
features=NULL, nor=NULL, nor_method= "LogNormalize"){
DefaultAssay(que) <- assay1
DefaultAssay(ref) <- assay2
#首先确保两者的基因名和总基因数一致，如果是跨物种则需要进行同源基因转换
c1 <- which(rownames(ref) %in% rownames(que))
ref <- subset(ref,features=rownames(ref)[c1])
c1 <- which(rownames(que) %in% rownames(ref))
que <- subset(que,features=rownames(que)[c1])
group1 <- group1
group2 <- group2
#查看参考数据集的细胞分组情况
ref$celltype.l1 <- ref[[group1]]
ref$celltype.l2 <- ref[[group2]]
p1 <- DimPlot(object = ref, reduction = "umap", group.by = "celltype.l2", label = TRUE, label.size = 3, repel = TRUE) + NoLegend() + ggtitle(group2)
+scale_color_manual(values=cor)
p2 <- DimPlot(object = ref, reduction = "umap", group.by = "celltype.l1", label = TRUE, label.size = 3, repel = TRUE) + NoLegend() + ggtitle(group1)
+scale_color_manual(values=cor)
p1 <- p1 - p2
pf <- paste0(group1,"_",group2,'_ref_umap.jpg')
jpeg(pf,1800,1000)
print(p1)
dev.off()
#对数据集进行标准化，建议提前做好标准化，因为这一步十分消耗时间和内存，容易报错
if (is.null(nor)) {
if (!is.null(sct)) {
que <- SCTransform(que, verbose = FALSE, vars.to.regress = "percent.mt")
ref <- SCTransform(ref, verbose = FALSE, vars.to.regress = "percent.mt")
nor_method <- "SCT"
}else{
que <- NormalizeData(que, normalization.method = "LogNormalize", scale.factor = 10000)
ref <- NormalizeData(ref, normalization.method = "LogNormalize", scale.factor = 10000)
nor_method <- "LogNormalize"
}
saveRDS(que,"query_nor.rds")
saveRDS(ref,"ref_nor.rds")
}else{
nor_method <- nor_method
}
#如果指定了基因集，则使用指定的基因集进行映射
if (is.null(features)) {
features <- rownames(que)
}
#以下进行映射的常规流程
anchors <- FindTransferAnchors(
reference = ref,
query = que,
features = features,
normalization.method = nor_method,
#normalization.method = "SCT",
#reference.reduction = "pca",
#reference.reduction = "spca",
dims = 1:30
)
#给询问数据集加上映射信息
predictions.assay <- TransferData(anchorset = anchors, refdata = ref@meta.data[,group1], prediction.assay = TRUE,
weight.reduction = que[["pca"]], dims = 1:30)
que[["predictions"]] <- predictions.assay
DefaultAssay(que) <- "predictions"
#保存映射信息
predictions <- TransferData(anchorset = anchors, refdata = ref@meta.data[,group1], prediction.assay = FALSE,
weight.reduction = que[["pca"]], dims = 1:30)
predic_mtrx <- as.data.frame(predictions)
write.table(predic_mtrx,"predic_mtrx.tsv",sep = '\t',row.names = T,col.names = T,quote = F)
#可视化映射结果1
gene <- unique(ref@meta.data[,group1])
png(paste0(group1,"_transfer.png"),width = 2100, height = 1500)
p1 <- FeaturePlot(que,features = gene,ncol=3)
print(p1)
dev.off()
#可视化映射结果2
que$predict_Celltype <- predic_mtrx$predicted.id
dimplot_2groups(que,group2="predict_Celltype")
#返回存有映射结果的Seurat对象
return(que)
}
使用示例
没有特殊要求的话，可以直接传入三个参数即可
ob1 <- transfer_cell(ref=ob1,que=ob2,group1="Celltype")
saveRDS(ob1,"final_result.rds")
参考数据集可视化图
ref
可视化结果1（有点类似于FeaturePlot）
结果1
可视化结果2
结果2
此外，还有一个predic_mtrx.tsv文件和final_result.rds文件。predic_mtrx.tsv文件第二列为预测的细胞类型，其它列存了预测分数；final_result.rds是包含预测信息的询问数据集Seurat对象，预测信息在predict_Celltype列。
小结与补充
根据实践结果来看，Seurat的映射结果还是比较可靠的，也可以根据结果进行人工修正。一般，我们会去掉线粒体基因进行映射，另外参考数据集如果太大，可以取子集，不然好很耗内存，而且运行很慢。
