scmap运行过程

1.
导入singleR的参考数据
library(celldex)
ref <- celldex::DatabaseImmuneCellExpressionData()

refdata <-  get(load("/home/wy/.cache/R/ExperimentHub/ref_Human_all.RData"))

2.
将参考数据结构进行修改
细胞类型标签位于列名为“cell_type1”的列中
colData(ref)$cell_type1 <- colData(ref)$label.fine
基因名称位于列名为“feature_symbol”的列中
rowData(ref)$feature_symbol <- rownames(ref)
数据格式转化为SingleCellExperiment的对象
ref_sce <- SingleCellExperiment::SingleCellExperiment(assays=list(logcounts=Matrix::Matrix(assays(ref)$logcounts)),colData=colData(ref),rowData=rowData(ref))

3.（参考数据集为自定义的数据集）从summarizedexperiment转化为singlecellexperiment对象

library(SingleCellExperiment)
 library(SummarizedExperiment)

loaded_seurat_obj <- readRDS(file = "/media/public1/wy/singleR/Segerstolpe/Segerstolpe.ref2.rds")

loaded_seurat_obj2 <- as(loaded_seurat_obj, "SingleCellExperiment")

colData(loaded_seurat_obj2)$cell_type1 <- colData(loaded_seurat_obj2)$ref_label

rowData(loaded_seurat_obj2)$feature_symbol <- rownames(loaded_seurat_obj2)

ref_sce <- loaded_seurat_obj2

3.
library(scmap)
scmap挑出信息量最大的基因
ref_sce <- scmap::selectFeatures(ref_sce,suppress_plot=FALSE)
查看挑选的前50个基因
rownames(ref_sce)[which(rowData(ref_sce)$scmap_features)][1:50]
查看信息量最大features有多少
length(rownames(ref_sce)[which(rowData(ref_sce)$scmap_features)])

4.
可以自己加标记基因
创建自己感兴趣的标记基因
my_key_markers = c("TRAC","TRBC1","TRBC2","IGKC")
将创建的标记基因加入到scmap挑选的基因中
rowData(ref_sce)$scmap_features[rownames(ref_sce) %in% my_key_markers] <- TRUE
查看scmap 挑选的基因是否发生了变化
length(rownames(ref_sce)[which(rowData(ref_sce)$scmap_features)])

同样可以删除基因
mt_genes <- rownames(ref_sce)[grep("^MT-",rownames(ref_sce))]
rowData(ref_sce)$scmap_features[rownames(ref_sce) %in% mt_genes] <- FALSE
length(rownames(ref_sce)[which(rowData(ref_sce)$scmap_features)])

5.
为scmap feature后指定一个新变量
scmap_feature_genes <- rownames(ref_sce)[which(rowData(ref_sce)$scmap_features)]
length(scmap_feature_genes)

6.
构建scmap-cluster 使用的参考文件，用于细胞类型注释
ref_sce <- scmap::indexCluster(ref_sce)
格式化数据，让它可以作为ggplot2的输入，方便做图可视化特征
cormat <- reshape2::melt(as.matrix(metadata(ref_sce)$scmap_cluster_index))
画图可视化特征
ggplot2::ggplot(cormat,ggplot2::aes(x=Var2,y=Var1,fill=value))+ggplot2::geom_tile()+ggplot2::scale_fill_gradient2(low="blue",high="darkred",name="Expression value")+ggplot2::theme_minimal()+ggplot2::theme(axis.text.x=ggplot2::element_text(angle=90,vjust=1,size=18,hjust=1),axis.text.y=ggplot2::element_text(size=15),axis.title.x=ggplot2::element_blank(),axis.title.y=ggplot2::element_blank())

使用以上得到的ref_sce 对象来构建scmap-cell索引
ref_sce <- scmap::indexCell(ref_sce)

7.
在注释之前要完成正常的预处理，包括 QC 过滤、标准化和对数转换
对象为seurat，转化为适合scmap使用的对象
query_sce <- SingleCellExperiment::SingleCellExperiment(assays=list(counts=data))

构建seurat对象
query_seur <- Seurat::CreateSeuratObject(data)

query_seur <- readRDS(file = "/home/wy/xiangmu/mymodel_data/final_data/pancreas/Baron(GSE84133)/Baron_R.rds")
数据格式转换
query_sce <- Seurat::as.SingleCellExperiment(query_seur)

library(scater)
使用scater包进行标准化
query_sce <- scater::logNormCounts(query_sce)
增加名为features_symbol的列
rowData(query_sce)$feature_symbol <-rownames(query_sce)

8.
用scmap-cluster注释查询数据
scmap_cluster_res <- scmap::scmapCluster(projection=query_sce,index_list=list(ref_pancreas = metadata(ref_sce)$scmap_cluster_index),threshold=0.1)

scmap_cluster_labs <- scmap_cluster_res[[1]]

write.csv(scmap_cluster_labs, file = "/media/public1/wy/scmap_cluster/pancreas/baron_scmap_cluster_labs.csv", row.names = FALSE)

统计注释结果
par(mar=c(13,4,1,0))
barplot(table(scmap_cluster_res$combined_labs),las=2)

存储注释结果，可视化分群注释结果
colData(query_sce)$scmap_cluster <-scmap_cluster_res$combined_labs
query_sce <-scater::runUMAP(query_sce)
scater::plotReducedDim(query_sce,dimred="UMAP",colour_by="scmap_cluster")

9.
scmap-cell

nearest_neighbours <- scmap::scmapCell(projection=query_sce,index_list=list(ref_pancreas_cell = metadata(ref_sce)$scmap_cell_index),w=10)

model_label <- function(neighbours,metadata=ref_sce$cell_type1){
+ freq <- table(metadata[neighbours])
+ label <- names(freq)[which(freq ==max(freq))]
+ if (length(label) >1) {return("ambiguous")}
+ return(label)
+ }

scmap_cell_labs <- apply(nearest_neighbours$ref_pancreas_cell$cells,2,model_label)

write.csv(scmap_cluster_labs, file = "/media/public1/wy/scmap_cell/pancreas/baron_scmap_cell_labs.csv", row.names = FALSE)