Chapter 3 scRNA-seq analysis

This chapter describes the workflows for individual sample analysis, sample integration, marker visualization, and reference-based cell annotation and mapping.

3.1 Single sample

Each sample is analyzed independently using Seurat. The filtered feature-barcode matrix generated by Cell Ranger is imported, and sample identifiers are added to the cell metadata

library(Seurat)
library(dplyr)
library(Matrix)
library(gplots)
library(matrixStats)
library(sva)
library(ggpubr)
library(openxlsx)
library(stringr)
library(scran)
library(ggthemes)
library(ggthemes)
library(grDevices)
library(reticulate)
library(Biobase)
library(scatterplot3d)
library(monocle)
library(pheatmap)
library(harmony)
library(SingleR)


out.path <- paste0("1.scRNA/20250801.", sample.name)
system(sprintf("mkdir %s", out.path))


sce.data <- Read10X(data.dir = paste0("/public/home/daiyt/projects/bmsc/analysis/nrctm/scell/cellranger_9/", sample.name.raw, "_RNA/outs/filtered_feature_bc_matrix") )
dim(sce.data)
colnames(sce.data) <- str_replace_all(colnames(sce.data), "1", sample.name)
message(sample.name, "  ", ncol(sce.data))



############ analysis
sce <- CreateSeuratObject(counts = sce.data, project = "sce", min.cells = 0, min.features = 0)
sce[["percent.mt"]] <- PercentageFeatureSet(sce, pattern = "^MT-")
sce
summary(sce@meta.data$percent.mt)

sce@meta.data$Sample = sample.name



sce <- subset(sce, subset = nFeature_RNA > 200 & percent.mt < 20)
sce


sce <- NormalizeData(sce, normalization.method = "LogNormalize", scale.factor = 10000)
sce <- FindVariableFeatures(sce, selection.method = "vst", nfeatures = 2000)

# Identify the 10 most highly variable genes
top10 <- head(VariableFeatures(sce), 10)

# plot variable features with and without labels
pdf(paste0(out.path, "/3.VariableFeaturePlot.pdf"), width = 12, height = 7)
plot1 <- VariableFeaturePlot(sce)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1 + plot2
dev.off()

all.genes <- rownames(sce)
sce <- ScaleData(sce, features = VariableFeatures(sce))
sce <- RunPCA(sce, features = VariableFeatures(object = sce))

p <- DimPlot(sce, reduction = "pca") + theme_few()
ggsave(paste0(out.path, "/4.PCA.pdf"), p, width = 9, height = 7)


sce <- FindNeighbors(sce, reduction = "pca", dims = 1:30)
sce <- FindClusters(sce, resolution = 0.8)

sce <- RunTSNE(sce, reduction = "pca", dims = 1:30, perplexity = 50)
sce <- RunUMAP(sce, reduction = "pca", dims = 1:30)


p <- DimPlot(sce, reduction = "umap", pt.size = 1, 
             cols = color.lib) + theme_few()
ggsave(paste0(out.path, "/5.UMAP.cluster.pdf"), p, width = 9, height = 7)

p <- DimPlot(sce, reduction = "umap", pt.size = 1, label = TRUE, label.size = 10,
             cols = color.lib) + theme_few()
ggsave(paste0(out.path, "/5.UMAP.cluster.label.pdf"), p, width = 9, height = 7)


p <- DimPlot(sce, reduction = "tsne", pt.size = 1, 
             cols = color.lib) + theme_few()
ggsave(paste0(out.path, "/6.tSNE.cluster.pdf"), p, width = 9, height = 7)

p <- DimPlot(sce, reduction = "tsne", pt.size = 1, label = TRUE, label.size = 10,
             cols = color.lib) + theme_few()
ggsave(paste0(out.path, "/6.tSNE.cluster.label.pdf"), p, width = 9, height = 7)



############### Cell cycle
cc.genes

sce <- CellCycleScoring(sce, s.features = cc.genes$s.genes, g2m.features = cc.genes$g2m.genes)

p <- DimPlot(sce, reduction = "umap", pt.size = 1, label = FALSE, 
             group.by = "Phase", cols = color.lib) + theme_few()
ggsave(paste0(out.path, "/7.CC.phase.pdf"), p, width = 9, height = 7)

p <- VlnPlot(sce, features = "S.Score", pt.size = 0, 
             group.by = "seurat_clusters", cols = color.lib) + theme_few()
ggsave(paste0(out.path, "/7.CC.score.S.Score.pdf"), p, width = 10, height = 5)

p <- VlnPlot(sce, features = "G2M.Score", pt.size = 0, 
             group.by = "seurat_clusters", cols = color.lib) + theme_few()
ggsave(paste0(out.path, "/7.CC.score.G2M.Score.pdf"), p, width = 10, height = 5)



####################### SingleR annotation
ref <- readRDS(file = "hs.BlueprintEncodeData.RDS")
pred.BlueprintEncodeData <- SingleR(test = sce@assays$RNA$data, ref = ref, labels = ref$label.main)

ref <- readRDS(file = "hs.HumanPrimaryCellAtlasData.RDS")
pred.HumanPrimaryCellAtlasData <- SingleR(test = sce@assays$RNA$data, ref = ref, labels = ref$label.main)

ref <- readRDS(file = "NovershternHematopoieticData.RDS")
pred.NovershternHematopoieticData <- SingleR(test = sce@assays$RNA$data, ref = ref, labels = ref$label.main)


sce@meta.data$CellType.BlueprintEncodeData <- pred.BlueprintEncodeData$labels
sce@meta.data$CellType.HumanPrimaryCellAtlasData <- pred.HumanPrimaryCellAtlasData$labels
sce@meta.data$CellType.NovershternHematopoieticData <- pred.NovershternHematopoieticData$labels



p <- DimPlot(sce, reduction = "umap", pt.size = 1, label = TRUE, label.size = 4,
             group.by = "CellType.BlueprintEncodeData", cols = color.lib) + theme_few()
ggsave(paste0(out.path, "/8.SingleR.umap.BlueprintEncodeData.pdf"), p, width = 9, height = 7)

p <- DimPlot(sce, reduction = "umap", pt.size = 1, label = TRUE, label.size = 4,
             group.by = "CellType.HumanPrimaryCellAtlasData", cols = color.lib) + theme_few()
ggsave(paste0(out.path, "/8.SingleR.umap.HumanPrimaryCellAtlasData.pdf"), p, width = 9, height = 7)

p <- DimPlot(sce, reduction = "umap", pt.size = 1, label = TRUE, label.size = 4,
             group.by = "CellType.NovershternHematopoieticData", cols = color.lib) + theme_few()
ggsave(paste0(out.path, "/8.SingleR.umap.NovershternHematopoieticData.pdf"), p, width = 9, height = 7)


p <- DimPlot(sce, reduction = "tsne", pt.size = 1, label = TRUE, label.size = 4,
             group.by = "CellType.BlueprintEncodeData", cols = color.lib) + theme_few()
ggsave(paste0(out.path, "/8.SingleR.tsne.BlueprintEncodeData.pdf"), p, width = 9, height = 7)

p <- DimPlot(sce, reduction = "tsne", pt.size = 1, label = TRUE, label.size = 4,
             group.by = "CellType.HumanPrimaryCellAtlasData", cols = color.lib) + theme_few()
ggsave(paste0(out.path, "/8.SingleR.tsne.HumanPrimaryCellAtlasData.pdf"), p, width = 9, height = 7)

p <- DimPlot(sce, reduction = "tsne", pt.size = 1, label = TRUE, label.size = 4,
             group.by = "CellType.NovershternHematopoieticData", cols = color.lib) + theme_few()
ggsave(paste0(out.path, "/8.SingleR.tsne.NovershternHematopoieticData.pdf"), p, width = 9, height = 7)


p <- FeaturePlot(object = sce, features = c("CD34","MPO","SDC1","PTPRC","CD3G","CD79A","CA1","NGFR","AVP"), 
              cols = c("#CCCCCC", "red"), pt.size = 0.3, raster = FALSE, order = TRUE,
              reduction = "umap", ncol = 3)
ggsave(paste0(out.path, "/9.marker.umap.pdf"), p, width = 13, height = 12)
  

saveRDS(sce, file = paste0("obj/", sample.name ,"_RNA.rds"))

3.2 Merge sample

This section combines sample-level count matrices and metadata for joint analysis. Previously generated doublet labels are imported from sample-specific annotation files. The filtering criteria specify more than 300 detected genes, fewer than 50,000 RNA counts, a mitochondrial transcript percentage below 20%, and a non-doublet label.

length(sample.list)
length(sample.list %>% unique() )

sel.name <- c("orig.ident","nCount_RNA","nFeature_RNA","percent.mt","Sample","S.Score","G2M.Score","Phase","CellType","CellType.BlueprintEncodeData","CellType.HumanPrimaryCellAtlasData","CellType.NovershternHematopoieticData","Doublet")
for (kkk in 1:length(sample.list)) {
message(kkk, "   ", sample.list[kkk])

obj.raw <- readRDS(paste0("obj/", sample.list[kkk], "_RNA.rds"))
demux <- read.xlsx(paste0("obj/", sample.list[kkk], "_demux.xlsx"))

obj.raw@meta.data$CellType <- NA
obj.raw@meta.data$Doublet <- demux$Doublet[match(rownames(obj.raw@meta.data), demux$Cell)]

exp.data <- cbind(exp.data, obj.raw@assays$RNA$counts )
meta.data <- rbind(meta.data, obj.raw@meta.data[, sel.name] )

}

table(meta.data$Sample)
length(table(meta.data$Sample))

##############################
# Analysis and add meta data
sce <- CreateSeuratObject(counts = exp.data, project = "sce", min.cells = 0, min.features = 0)
sce[["percent.mt"]] <- PercentageFeatureSet(sce, pattern = "^MT-")
sce
summary(sce@meta.data$percent.mt)

sce@meta.data$SampleID <- meta.data$Sample[match(rownames(sce@meta.data), rownames(meta.data) )]
length(unique(sce@meta.data$SampleID))
table(sce@meta.data$Sample)

data.frame(table(sce$SampleID))
sce <- subset(sce, subset = nFeature_RNA > 300 & percent.mt < 20 & nCount_RNA < 50000 & Doublet == FALSE)
data.frame(table(sce$SampleID))



############ analysis
pdf(paste0(out.path, "/2.filter.vlnplot.pdf"), width = 16, height = 7)
VlnPlot(sce, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), pt.size = 0, group.by = "SampleID", ncol = 3)
dev.off()



pdf(paste0(out.path, "/2.filter.geneplot.pdf"), width = 14, height = 7)
plot1 <- FeatureScatter(sce, feature1 = "nCount_RNA", feature2 = "percent.mt", group.by = "SampleID", shuffle = T)
plot2 <- FeatureScatter(sce, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", group.by = "SampleID", shuffle = T)
plot1 + plot2
dev.off()


sce <- NormalizeData(sce, normalization.method = "LogNormalize", scale.factor = 100000 )


sce <- FindVariableFeatures(sce, selection.method = "vst", nfeatures = 2500)
# Identify the 10 most highly variable genes
top10 <- head(VariableFeatures(sce), 10)
top10
VariableFeatures(sce) <- VariableFeatures(sce)[! VariableFeatures(sce) %in% gene.id.ig$V6]
VariableFeatures(sce) <- VariableFeatures(sce)[! VariableFeatures(sce) %in% gene.id.tr$V6]
VariableFeatures(sce) <- VariableFeatures(sce)[! VariableFeatures(sce) %in% gene.id.hb$V6]
VariableFeatures(sce) <- VariableFeatures(sce)[! VariableFeatures(sce) %in% gene.id.rp$V6]
VariableFeatures(sce) <- VariableFeatures(sce)[! VariableFeatures(sce) %in% cc.genes$s.genes]
VariableFeatures(sce) <- VariableFeatures(sce)[! VariableFeatures(sce) %in% cc.genes$g2m.genes]
length(VariableFeatures(sce))
top10 <- head(VariableFeatures(sce), 10)
top10
VariableFeatures(sce) <- head(VariableFeatures(sce), 2000)


all.genes <- rownames(sce)
sce <- ScaleData(sce, features = VariableFeatures(sce))
sce <- RunPCA(sce, features = VariableFeatures(object = sce))


# plot variable features with and without labels
pdf(paste0(out.path, "/3.VariableFeaturePlot.pdf"), width = 12, height = 7)
plot1 <- VariableFeaturePlot(sce)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
plot1 + plot2
dev.off()

all.genes <- rownames(sce)
sce <- ScaleData(sce, features = VariableFeatures(sce))
sce <- RunPCA(sce, features = VariableFeatures(object = sce))

p <- DimPlot(sce, reduction = "pca", group.by = "SampleID", cols = color.lib, shuffle = T) + theme_few()
ggsave(paste0(out.path, "/4.PCA.pdf"), p, width = 11, height = 7)


length(VariableFeatures(sce))
sce <- RunHarmony(sce, "SampleID", max_iter = 3, sigma = 0.050)
sce <- FindNeighbors(sce, reduction = "harmony", dims = 1:20)
sce <- FindClusters(sce, resolution = 1)
set.seed(1)
sce <- RunUMAP(sce, reduction = "harmony", dims = 1:20)
set.seed(1)
kmeans.cluster <- stats::kmeans(sce@reductions$harmony@cell.embeddings[, 1:20], centers = 100, iter.max = 100)
summary(as.numeric(table(kmeans.cluster$cluster)))
sce$Kmeans <- kmeans.cluster$cluster
sce$Cluster <- paste0("C", str_pad(sce$seurat_clusters, 2, "left", "0") )



pt.size = 0.005
set.raster = FALSE

p <- DimPlot(sce, reduction = "umap", pt.size = pt.size, raster = set.raster, label = F, label.size = 5,
             group.by = "SampleID", cols = c(color.lib,color.lib), shuffle = T) + theme_few()
ggsave(paste0(out.path, "/5.UMAP.SampleID.png"), p, width = 15, height = 7, dpi = 300)

p <- DimPlot(sce, reduction = "umap", pt.size = pt.size, raster = set.raster, label = T, label.size = 5,
             group.by = "Cluster", cols = c(color.lib,color.lib) ) + theme_few()
ggsave(paste0(out.path, "/5.UMAP.cluster.png"), p, width = 10, height = 7, dpi = 300)

3.3 Marker plot

For each gene, the code calculates the mean expression and the percentage of cells with expression values greater than zero within each cell type. Gene order and cell-type order follow the supplied marker list and annotation palette.

out.path.plot <- paste0(out.path, "/DotMarker")
cmd <- sprintf("mkdir %s", out.path.plot)
system(cmd)

gene.list <- c("AVP","SPINK2","CD34","GATA2","ITGA2B","GATA1","PF4","MPO","AZU1","ELANE","RNASE2","LYZ","DEFA4","CAMP","CRISP3","MMP8","MMP9","CMTM2","FCGR3B","CD14","FCN1","S100A9","C1QA","C1QB","CD1C","CLEC10A","LILRA4","IL3RA","IRF8","HMBS","KLF1","ALAS2","IFIT1B","TMCC2","DNTT","IGLL1","MME","MKI67","CD79A","DTX1","ACSM3","MS4A1","BANK1","FCER2","SDC1","TNFRSF13B","CD3E","CD8A","GZMH","VWF","CDH5","FABP4","RGS5","ACTA2","COMP","DCN","FN1","COL3A1","COL1A1","CXCL12","CHL1","ALPL","IBSP","SPP1","BGLAP")
length(gene.list)


plot.data <- cbind(sce@meta.data, FetchData(sce, vars = c("umap_1","umap_2", gene.list) ))
plot.data$Group <- plot.data$CellType
plot.mean <- aggregate(plot.data[, gene.list], list(Type = plot.data$Group), mean )
plot.perc <- aggregate(plot.data[, gene.list] > 0, list(Type = plot.data$Group), mean )

plot.info <- NULL
for (i in 1:length(gene.list)) {
  sub <- data.frame(Type = plot.mean$Type, 
                    MeanExp = plot.mean[, gene.list[i]], 
                    MeanScale = scale(as.numeric(plot.mean[, gene.list[i]])), 
                    Percent = plot.perc[, gene.list[i]] * 100, Gene = gene.list[i] ) 
  plot.info <- rbind(plot.info, sub)
}
plot.info$Type <- factor(as.character(plot.info$Type), levels = rev(names(color.bmsc)) )
plot.info$Gene <- factor(as.character(plot.info$Gene), levels = c(gene.list) )
p <- ggscatter(plot.info, x = "Gene", y = "Type", 
               color = "MeanScale", size = "Percent", shape = 16,
               #order = order.all.sample.2, 
               xlab = "", ylab = "")
p <- p + theme_blank() + theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5, size = 8))
p <- p + theme(axis.text.y = element_text(size = 12))
p <- p + scale_size(range = c(0, 6))
p <- p + gradient_color(c("#93e0ff","#c1f4ff","#EEEEEE","#ffb459","#e8613c","#b70909","#b70909","#b70909"))
ggsave(paste0(out.path.plot, "/20260701.marker_s42_scale.pdf"), p, width = 14, height = 7)

3.4 KNN

This section compares selected query populations with an annotated hematopoietic reference. Reference and query cells are combined for variable feature selection, scaling, and PCA. A separate nearest-neighbor procedure generates coordinates for displaying query cells on the reference UMAP. This step uses the existing CellTypePred annotations. Within each annotated cell type, up to 1,000 reference cells are retained, and up to 100 reference neighbors are identified for each query cell in PCA space. One neighbor is randomly selected, and its reference UMAP coordinates are assigned to the query cell.

Gaussian jitter with a standard deviation of 0.1 is added to reduce point overlap. The resulting coordinates are used for visualization, and the query-cell plot is exported with cells colored by CellTypePred.

out.path.knn <- paste0(out.path, "/KNN_hspc")
cmd <- sprintf("mkdir %s", out.path.knn)
system(cmd)


sce.ref <- readRDS("obj/sce.hemato.rds")
color.this <- c(
"HSC_LMPP" = "#CE0000",
"MEP_MKP" = "#ffb82b",
"GMP_MDP" = "#ed4e80",
"Early_myeloid" = "#ff9397",
"Late_myeloid" = "#ffbfc3",
"Neutrophil" = "#ffdbde",
"Monocyte" = "#aa3a79",
"Macrophage" = "#7c1540",
"cDC" = "#ba784c",
"pDC" = "#aa4400",
"Proerythroblast" = "#ffd57a",
"Erythroblast" = "#ffef99",
"Pro_B" = "#c4b2ff",
"Immature_B" = "#8860ff",
"Mature_B" = "#5c2ce0")

p <- DimPlot(sce.ref, reduction = "umap", pt.size = 0.2, raster = FALSE, label = T, label.size = 3,
             group.by = "CellType", cols = color.this, shuffle = T) + theme_blank()
ggsave(paste0(out.path.knn, "/1.ref.CellType.png"), p, width = 8.5, height = 7, dpi = 300)


sce.query <- subset(sce, CellType %in% c("Prog_mye","Macrophage","Erythroblast","B_cell") )
p <- DimPlot(sce.query, reduction = "umap", pt.size = 0.3, raster = FALSE, label = T, label.size = 3,
             group.by = "CellType", cols = color.spin, shuffle = T) + theme_blank()
ggsave(paste0(out.path.knn, "/1.query.CellType.png"), p, width = 8.5, height = 7, dpi = 300)



common.col <- intersect( colnames(sce.ref@meta.data), colnames(sce.query@meta.data) )
meta.data <- rbind( sce.ref@meta.data[, common.col], sce.query@meta.data[, common.col] )
count.data <- cbind( sce.ref@assays$RNA$count, sce.query@assays$RNA$count )


sce.merge <- CreateSeuratObject(counts = count.data, project = "sce", min.cells = 0, min.features = 0)
sce.merge@meta.data <- cbind(sce.merge@meta.data, meta.data[, 4:ncol(meta.data)] )
sce.merge <- NormalizeData(sce.merge, normalization.method = "LogNormalize", scale.factor = ncol(sce) )
sce.merge@assays$RNA$data <- cbind( sce.ref@assays$RNA$data, sce.query@assays$RNA$data )
table(sce.merge$Tissue)


sce.merge <- FindVariableFeatures(sce.merge, selection.method = "vst", nfeatures = 2500)
VariableFeatures(sce.merge) <- VariableFeatures(sce.merge)[! VariableFeatures(sce.merge) %in% gene.id.ig$V6]
VariableFeatures(sce.merge) <- VariableFeatures(sce.merge)[! VariableFeatures(sce.merge) %in% gene.id.tr$V6]
VariableFeatures(sce.merge) <- VariableFeatures(sce.merge)[! VariableFeatures(sce.merge) %in% gene.id.hb$V6]
VariableFeatures(sce.merge) <- VariableFeatures(sce.merge)[! VariableFeatures(sce.merge) %in% gene.id.rp$V6]
VariableFeatures(sce.merge) <- VariableFeatures(sce.merge)[! VariableFeatures(sce.merge) %in% cc.genes$s.genes]
VariableFeatures(sce.merge) <- VariableFeatures(sce.merge)[! VariableFeatures(sce.merge) %in% cc.genes$g2m.genes]
VariableFeatures(sce.merge) <- head(VariableFeatures(sce.merge), 2000)
length(VariableFeatures(sce.merge))

sce.merge <- ScaleData(sce.merge, features = VariableFeatures(sce.merge))
sce.merge <- RunPCA(sce.merge, features = VariableFeatures(object = sce.merge))

table(sce.merge$Tissue)


sce.ref.sub <- subset(sce.merge, cells = colnames(sce.ref) )
sce.query.sub <- subset(sce.merge, cells = colnames(sce.query) )

sce.anchors <- FindTransferAnchors(reference = sce.ref.sub, query = sce.query.sub, dims = 1:20, reference.reduction = "pca")

pred <- TransferData(anchorset = sce.anchors, refdata = sce.ref.sub$CellType, dims = 1:20)
sce.query <- AddMetaData(sce.query, pred)

table(pred$predicted.id)
table(sce.query$predicted.id)
table(sce.query$CellTypePred)



sum(sce.query$CellTypePred != sce.query$predicted.id)

#sce.query$CellTypePred <- sce.query$predicted.id
table(sce.query@meta.data[, c("CellTypePred","AgeRange")])




library(FNN)
library(BiocNeighbors)

k.search <- 100
max_ref_cells_per_type <- 1000
pca_dims <- 1:20
jitter_size <- 0.1


ref.pca <- Embeddings(sce.ref.sub, "pca")[ , pca_dims]
query.pca <- Embeddings(sce.query.sub, "pca")[ , pca_dims]
ref.umap <- cbind(UMAP_1 = sce.ref$umap_anno_1, UMAP_2 = sce.ref$umap_anno_2)
rownames(ref.umap) <- colnames(sce.ref)


cell_types <- unique(sce.query$CellTypePred)
final_query_umap <- vector("list", length(cell_types))
names(final_query_umap) <- cell_types


for (ct in cell_types) {
  ref_idx <- which(sce.ref$CellType == ct)
  query_idx <- which(sce.query$CellTypePred == ct)
  query_cells <- colnames(sce.query)[query_idx]

  if (length(ref_idx) > max_ref_cells_per_type) {
    ref_idx <- sample(ref_idx, max_ref_cells_per_type)
  }

  k_eff <- min(k.search, length(ref_idx))
  knn <- queryKNN(
    query = query.pca[query_idx, , drop = FALSE],
    X = ref.pca[ref_idx, , drop = FALSE],
    k = k_eff,
    BNPARAM = AnnoyParam()
  )

  selected_idx <- apply(knn$index, 1, \(x) sample(x[!is.na(x)], 1))
  query_umap <- ref.umap[ref_idx[selected_idx], , drop = FALSE]
  
  rownames(query_umap) <- query_cells
  final_query_umap[[ct]] <- query_umap
}


query.umap <- do.call(rbind, final_query_umap)
query.umap <- query.umap[colnames(sce.query), ]


set.seed(1)
sce.query$umap_pred_1 <- query.umap[,1] + rnorm(nrow(query.umap), sd = jitter_size)
sce.query$umap_pred_2 <- query.umap[,2] + rnorm(nrow(query.umap), sd = jitter_size)




color.this <- color.bmsc[1:17]
color.this

sce.query$CellTypePred <- factor(as.character(sce.query$CellTypePred), levels = names(color.this))
sce.query@reductions$umap@cell.embeddings[, "umap_1"] <- sce.query$umap_pred_1 
sce.query@reductions$umap@cell.embeddings[, "umap_2"] <- sce.query$umap_pred_2 

p <- DimPlot(sce.query, reduction = "umap", pt.size = 0.5, raster = FALSE, label = F, label.size = 5,
             group.by = "CellTypePred", shuffle = T, cols = color.this) + theme_blank()
ggsave(paste0(out.path.knn, "/2.CellTypePred.png"), p, width = 8.8, height = 7, dpi = 300)