Chapter 4 Spatial

This chapter covers spatial transcriptomic data preprocessing, clustering, visualization, and cell-type signature scoring.

4.1 Preprocessing

Visium HD outputs are imported into Seurat, and spatial coordinates and tissue images are prepared for quality assessment and visualization. The workflow includes expression normalization, variable gene selection, dimensionality reduction, clustering, and visualization of clusters and selected gene expression patterns in tissue sections.

sample.name <- ""
v3d <- Load10X_Spatial(data.dir = paste0("spaceranger_4/", sample.name,"_STR/outs/"), slice = "slice1", bin.size = c(8, 16, "polygons") )
v3d$SampleID <- sample.name

v3d@assays

v3d.img <- image_read( paste0("", sample.name, "/", sample.name, "_PIC.jpg" ) ) 
message(paste("Pic size:", image_info(v3d.img)$width, "x", image_info(v3d.img)$height))


out.path <- paste0("spatial_", sample.name)
system(sprintf("mkdir %s", out.path))


#############################################################
####   Normalize coord
#############################################################
DefaultAssay(v3d) <- "Spatial.008um"
coord.dott <- cbind(v3d@meta.data[!is.na(v3d$nCount_Spatial.008um), ], GetTissueCoordinates(v3d) )
coord.dott$width_nor <- (coord.dott$y - min(coord.dott$y)) / max(coord.dott$y - min(coord.dott$y)) * 100
coord.dott$height_nor <- abs( max(coord.dott$x) - coord.dott$x ) / max( abs(max(coord.dott$x) - coord.dott$x) ) * 100 

DefaultAssay(v3d) <- "Spatial.016um"
coord.rect <- cbind(v3d@meta.data[!is.na(v3d$nCount_Spatial.016um), ], GetTissueCoordinates(v3d) )
coord.rect$width_nor <- (coord.rect$y - min(coord.rect$y)) / max(coord.rect$y - min(coord.rect$y)) * 100
coord.rect$height_nor <- abs( max(coord.rect$x) - coord.rect$x ) / max( abs(max(coord.rect$x) - coord.rect$x) ) * 100 

DefaultAssay(v3d) <- "Spatial.Polygons"
coord.core <- cbind(v3d@meta.data[!is.na(v3d$nCount_Spatial.Polygons), ], GetTissueCoordinates(v3d) )
coord.core$width_nor <- (coord.core$x - min(coord.core$x)) / max(coord.core$x - min(coord.core$x)) * 100
coord.core$height_nor <- (coord.core$y - min(coord.core$y)) / max(coord.core$y - min(coord.core$y)) * 100


#############################################################
####   Crop pic
#############################################################
summary(coord.rect$y)
summary(coord.rect$x)

image_info(v3d.img)$width
image_info(v3d.img)$height

img.width <- max(coord.rect$y) - min(coord.rect$y)  
img.height <-  max(coord.rect$x) - min(coord.rect$x)
geometry <- paste0(img.width, "x", img.height, "+", min(coord.rect$y),"+", min(coord.rect$x))

img.cropped <- v3d.img %>% image_crop(geometry = geometry)

image_info(img.cropped)$width
image_info(img.cropped)$height


image_write(img.cropped, paste0(out.path, "/0.Pic_cropped.png") )
img.raster <- as.raster(img.cropped)



#############################################################
####  basic analysis
#############################################################
coord.core$QC <- "YES"
coord.core$QC[which(coord.core$nCount_Spatial.Polygons < 5) ] <- "NO"
table(coord.core$QC)
p <- ggplot() + annotation_raster(img.raster, xmin = 0, xmax = 100, ymin = 0, ymax = 100)
p <- p + geom_point(data = coord.core, aes(x = width_nor, y = height_nor, color = QC, fill = QC), size = 0.5, shape = 16)
p <- p + ggtitle(paste0("QC YES: ",  table(coord.core$QC)["YES"], "; NO: ",  table(coord.core$QC)["NO"] ) )
p <- p + scale_color_manual(values = c("#CE0000", "black"))
p <- p + theme_blank() + coord_fixed(ratio = 1) + guides(color = guide_legend(override.aes = list(size = 5))) 
ggsave(paste0(out.path, "/1.nCount_Spatial.Polygons_QC.png"), p, width = 9.2, height = 8, dpi = 300)
saveRDS(coord.core, paste0(out.path, "/1.nCount_Spatial.Polygons_QC.rds"))


coord.rect$QC <- "YES"
coord.rect$QC[which(coord.rect$nCount_Spatial.016um < 3) ] <- "NO"
table(coord.rect$QC)
p <- ggplot() + annotation_raster(img.raster, xmin = 0, xmax = 100, ymin = 0, ymax = 100)
p <- p + geom_point(data = coord.rect, aes(x = width_nor, y = height_nor, color = QC, fill = QC), size = 0.05, shape = 15)
p <- p + ggtitle(paste0("QC YES: ",  table(coord.rect$QC)["YES"], "; NO: ",  table(coord.rect$QC)["NO"] ) )
p <- p + scale_color_manual(values = c("#CE0000", "black"))
p <- p + theme_blank() + coord_fixed(ratio = 1) + guides(color = guide_legend(override.aes = list(size = 5))) 
ggsave(paste0(out.path, "/1.nCount_Spatial.016um_QC.png"), p, width = 9.2, height = 8, dpi = 300)
saveRDS(coord.rect, paste0(out.path, "/1.nCount_Spatial.016um_QC.rds"))


coord.dott$QC <- "YES"
coord.dott$QC[which(coord.dott$nCount_Spatial.008um < 3) ] <- "NO"
table(coord.dott$QC)
p <- ggplot() + annotation_raster(img.raster, xmin = 0, xmax = 100, ymin = 0, ymax = 100)
p <- p + geom_point(data = coord.dott, aes(x = width_nor, y = height_nor, color = QC, fill = QC), size = 0.05, shape = 15)
p <- p + ggtitle(paste0("QC YES: ",  table(coord.dott$QC)["YES"], "; NO: ",  table(coord.dott$QC)["NO"] ) )
p <- p + scale_color_manual(values = c("#CE0000", "black"))
p <- p + theme_blank() + coord_fixed(ratio = 1) + guides(color = guide_legend(override.aes = list(size = 5))) 
ggsave(paste0(out.path, "/1.nCount_Spatial.008um_QC.png"), p, width = 9.2, height = 8, dpi = 300)
saveRDS(coord.dott, paste0(out.path, "/1.nCount_Spatial.008um_QC.rds"))


exp.data <- v3d.sub@assays$Spatial.Polygons$counts
colnames(exp.data) <- str_replace_all(colnames(exp.data), "-1", paste0("-", sample.name) )
message(sample.name, "  ", ncol(exp.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$SampleID <- sample.name
sce$width_nor <- enroll.cell$width_nor
sce$height_nor <- enroll.cell$height_nor




#############################################################
####  basic analysis
#############################################################
sce <- NormalizeData(sce, normalization.method = "LogNormalize", scale.factor = 10000)



length(VariableFeatures(sce))
sce <- FindVariableFeatures(sce, selection.method = "vst", nfeatures = 2000)
length(VariableFeatures(sce))

coord.core <- sce@meta.data

sce$width_group <- cut(sce$width_nor, breaks = 10)
sce$height_group <- cut(sce$height_nor, breaks = 10)
sce$width_group <- paste0("W", str_pad(as.numeric(sce$width_group), 2, "left", "0") )
sce$height_group <- paste0("H", str_pad(as.numeric(sce$height_group), 2, "left", "0") )
sce$group <- paste0(sce$width_group, "_", sce$height_group)
length(unique(sce$group))

scell.subc.exp <- NULL
subc.name <- names(table(sce$group))
for (i in 1:length(subc.name)) {
  if (i %% 50 == 1) message(i)
  if (length(which(sce$group == subc.name[i] ) ) > 1) {
    sub <- rowMeans(sce@assays$RNA$data[, rownames(sce@meta.data)[which(sce$group == subc.name[i] )] ])
  } else {
    sub <- sce@assays$RNA$data[, rownames(sce@meta.data)[which(sce$group == subc.name[i] )] ]
  }
  scell.subc.exp <- cbind(scell.subc.exp, sub)
}
colnames(scell.subc.exp) <- subc.name
scell.subc.exp <- limma::normalizeBetweenArrays(scell.subc.exp)


sce.gene.vars <- data.frame(Gene = rownames(scell.subc.exp), Vars = rowVars(scell.subc.exp), 
                            MeanExp = rowMeans(scell.subc.exp))
sce.gene.vars <- sce.gene.vars[order(sce.gene.vars$Vars, decreasing = T), ]
sce.gene.vars$Index <- 1:nrow(sce.gene.vars)
sce.gene.vars$isVar <- 0
sce.gene.vars$isVar[sce.gene.vars$Gene %in% VariableFeatures(sce) ] <- 1
sce.gene.vars$isCoding <- 0
sce.gene.vars$isCoding[sce.gene.vars$Gene %in% gene.id.coding$V6 ] <- 1

sub.1 <- sce.gene.vars[which(sce.gene.vars$Index <= 2000 & sce.gene.vars$MeanExp > 0.01 & sce.gene.vars$isCoding == 1), ]
sub.2 <- sce.gene.vars[which(sce.gene.vars$Index <= 5000 & sce.gene.vars$MeanExp > 0.01 & sce.gene.vars$isCoding == 1 & sce.gene.vars$isVar == 1), ]

length(union(sub.1$Gene, sub.2$Gene))


VariableFeatures(sce) <- union(sub.1$Gene, sub.2$Gene)
length(VariableFeatures(sce))


saveRDS(sce.gene.vars, paste0(out.path, "/2.sce.gene.vars.rds")   )



#############################################################
####  cluster analysis
#############################################################
sce <- ScaleData(sce, features = VariableFeatures(sce))
sce <- RunPCA(sce, reduction.name = "pca")
sce <- FindNeighbors(sce, reduction = "pca", dims = 1:20)
sce <- RunUMAP(sce, reduction = "pca", reduction.name = "umap", dims = 1:20, n.neighbors = 30, min.dist = 0.1)
sce <- FindClusters(sce, resolution = 0.8)
sce <- FindClusters(sce, resolution = 1.0)
sce <- FindClusters(sce, resolution = 1.5)


sce$Cluster <- paste0("C", str_pad(sce$RNA_snn_res.0.8, 2, "left", "0") )
table(sce$Cluster)


set.seed(1)
kmeans.cluster <- stats::kmeans(sce@reductions$pca@cell.embeddings[, 1:20], centers = 50, iter.max = 100)
summary(as.numeric(table(kmeans.cluster$cluster)))
sce$Kmeans <- paste0("K", str_pad(kmeans.cluster$cluster, 2, "left", "0") ) 


#############################################################
####   Visualization
#############################################################
coord.core <- cbind(sce@meta.data, FetchData(sce, vars = c("umap_1","umap_2")) )

p1 <- DimPlot(sce, group.by = "Cluster", reduction = "umap", pt.size = 0.1, raster = FALSE,
              label = TRUE, label.size = 6, cols = color.lib) + theme_blank() 
p2 <- ggplot() + annotation_raster(img.raster, xmin = 0, xmax = 100, ymin = 0, ymax = 100) + 
      geom_point(data = coord.core, aes(x = width_nor, y = height_nor, color = Cluster, fill = Cluster), size = 0.5, shape = 16) + 
      ggtitle(paste0("Cluster" ) ) + scale_color_manual(values = color.lib ) + theme_blank() + 
      coord_fixed(ratio = 1) + guides(color = guide_legend(override.aes = list(size = 5))) 
p <- ggarrange(p1, p2, ncol = 2, nrow = 1, widths = c(1, 1) )
ggsave(paste0(out.path, "/3.UMAP.Cluster.png"), p, width = 15, height = 7, dpi = 300)



gexp <- cbind(sce@meta.data, FetchData(sce, vars = c("umap_1","umap_2", "CXCL12", "CXCR4" ) ) )
gexp$PlotTag <- "other"
gexp$PlotTag[which(gexp$CXCL12 > 0 )] <- "CXCL12"
gexp$PlotTag[which(gexp$CXCR4 > 0 )] <- "CXCR4"
table(gexp$PlotTag)
gexp$PlotTag <- factor(as.character(gexp$PlotTag), levels = c("other","CXCL12","CXCR4"))
gexp$PlotSize = 0.5
gexp$PlotSize[which(gexp$PlotTag != "other")] = 1
plot.data <- gexp
set.seed(1)
sub <- plot.data[sample(1:nrow(plot.data), nrow(plot.data) * 0.6), ]
sub <- sub[order(sub$PlotTag), ]
plot.data[rownames(plot.data) %in% rownames(sub), ] <- sub
p1 <- ggplot() + #annotation_raster(img.raster, xmin = 0, xmax = 100, ymin = 0, ymax = 100) + 
      geom_point(data = plot.data, aes(x = umap_1, y = umap_2, color = PlotTag, size = PlotSize) , shape = 16, alpha = 1) + 
      ggtitle(paste0( "CXCL12-CXCR4" ) ) + theme_blank() + scale_size(range = c(0.5, 1)) + 
      scale_color_manual(values = c("#DDDDDD","#7cb0ff","#cc0846"))
p2 <- ggplot() + #annotation_raster(img.raster, xmin = 0, xmax = 100, ymin = 0, ymax = 100) + 
      geom_point(data = plot.data, aes(x = width_nor, y = height_nor, color = PlotTag, size = PlotSize), shape = 16, alpha = 1) + 
      ggtitle(paste0( "CXCL12-CXCR4" ) ) + theme_blank() + scale_size(range = c(0.5, 1)) + 
      coord_fixed(ratio = 1) + scale_color_manual(values = c("#DDDDDD","#7cb0ff","#cc0846"))
p <- ggarrange(p1, p2, ncol = 2, nrow = 1, widths = c(1, 1) )
ggsave(paste0(out.path.plot, "/CXCL12-CXCR4.png"), p, width = 15, height = 7, dpi = 300)

4.2 AUCell

AUCell is used to score predefined cell-type marker gene sets in the spatial expression data. The resulting scores are saved for downstream annotation and analysis of spatial cell-type distributions.

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

library(AUCell)

marker.list.all <- list(
  HSC_LMPP = c("SPINK2", "PRSS57", "SMIM24", "CD34", "PROM1", "HLF", "HOXA9", "AVP", "BAALC", "CRHBP", "MMRN1", "MSI2", "ETV6", "RUNX1", "CDK6", "SATB1", "SELL", "TSC22D1", "ANKRD28", "SERPINB6"),
  MEP_MKP = c("ITGA2B", "GP9", "GP1BB", "MPIG6B", "PTGS1", "TREML1", "GATA1", "NFE2", "CMTM5", "RAB27B", "FERMT1", "STAT5A", "PDLIM1", "PRKAR2B", "SLC39A3", "PARVB", "TPM1", "NME1", "HADH", "SPN", "PF4"),
  GMP_MDP = c("MPO", "PRTN3", "CTSG", "MS4A3", "CEBPA", "AZU1", "LYZL6", "RNASE2", "RNASE3", "RAB32"),
  Early_myeloid = c("ELANE", "BPI", "DEFA4", "SLPI", "PGLYRP1", "LCN2", "LTF", "AZU1", "DEFA3", "CD24", "RETN", "CEBPE"),
  Late_myeloid = c("CAMP", "CRISP3", "MMP8", "OLFM4", "TCN1", "CEACAM8", "PADI4", "ARG1", "RETN", "ANXA3", "CYBB", "GCA", "MNDA", "CDA", "LTA4H", "NCF4"),
  Neutrophil = c("FCGR3B", "FPR1", "C5AR1", "CMTM2", "AQP9", "VNN2", "MMP25", "FCAR", "TREM1", "IL1R2", "CSF3R", "S100A12", "CXCL8", "G0S2", "BCL2A1", "NCF1", "MMP9"),
  Monocyte = c("CD14", "SLC11A1", "CD300E", "MS4A6A", "ANPEP", "NLRP3", "PLIN2", "LILRB2", "C5AR1", "NCF2"),
  Macrophage = c("C1QA", "C1QB", "C1QC", "CD163", "FOLR2", "VSIG4", "MSR1", "GPNMB", "STAB1", "MRC1", "LGMN", "CD68", "FCGRT"),
  cDC = c("CD1C", "FCER1A", "CPVL", "ZNF385A", "IL4I1", "CD1E", "FCGR2B", "CCL22", "CD83", "SAMHD1", "LGALS9", "SYNGR2", "CLEC10A", "IFI30", "CTSH", "AIF1", "CD300C", "FGR", "KLF4", "RAB31"),
  pDC = c("LILRA4", "CLEC4C", "IL3RA", "TCF4", "IRF7", "PLD4", "SCT", "LRRC26", "PACSIN1", "SLC15A4", "PTCRA", "ITM2C", "PLAC8", "GNA15", "PPP1R14B", "RGS1", "TSPAN13", "CYB561A3", "LIME1", "SH2B3", "AREG", "SLC7A5", "BHLHE40", "GPR183", "IRF8"),
  Proerythroblast = c("KLF1", "GFI1B", "HEMGN", "TMEM14C", "SLC25A37", "SLC25A39", "HMBS", "UROD", "FECH", "CD36", "TFRC", "CA1", "CA2", "HBM", "HBA2", "HBD", "EIF2AK1", "BLVRB", "PRDX2", "TESC", "NME4", "GLRX5", "BOLA3", "IDH2"),
  Erythroblast = c("GYPA", "GYPC", "SLC4A1", "RHAG", "EPB42", "SPTA1", "SPTB", "ANK1", "TRIM58", "AHSP", "BPGM", "ALAS2", "EPB41", "BNIP3L", "FECH", "KLF1", "HMBS", "HEMGN"),
  Pro_B = c("DNTT", "VPREB1", "VPREB3", "IGLL1", "RAG1", "RAG2", "LCN6", "MME", "EBF1", "TCF3", "CD79B", "SOX4", "MEF2C", "FOXO1", "SMIM3", "CD9", "S1PR4"),
  Immature_B = c("TCL1A", "FCRLA", "DTX1", "BACH2", "NIBAN3", "KLHL14", "ADAM23", "LAMP5", "CD24", "SYK", "CMTM7", "RCSD1", "MTSS1", "TLE1", "CD79B", "AFF3", "PAX5"),
  Mature_B = c("MS4A1", "CD22", "BANK1", "FCER2", "TNFRSF13C", "FCRL1", "FCRL2", "CD40", "CXCR5", "BLK", "CD69", "LTB", "FCMR", "BIRC3", "CD79A"),
  Plasma = c("PRDM1", "MZB1", "SDC1", "TNFRSF17", "TXNDC5", "FKBP11", "SEC11C", "PRDX4", "PDIA4", "SDF2L1", "MANF", "DNAJB9", "MYDGF", "LMAN1", "ERN1", "TXNDC11", "PDIA6", "PIM2", "ELL2", "RPN2", "POU2AF1", "GPRC5D"),
  T_NK_cell = c("CD3D", "CD3E", "CD3G", "CD247", "TRAC", "LCK", "ZAP70", "KLRD1", "KLRK1", "PRF1", "GZMA", "GZMH", "CCL5", "NKG7", "FYN", "ETS1", "CD48", "RHOH", "GNG2", "FYB1", "SOCS1"),
  Endothelial = c("PECAM1", "CDH5", "ECSCR", "VWF", "CLDN5", "MYCT1", "ESAM", "CCL14", "SHANK3", "PLVAP"),
  Vascular = c("TINAGL1", "COX4I2", "RGS5", "ACTA2", "SEPTIN4", "HIGD1B", "ADIRF", "GJA4", "MYH11"),
  Fibroblast = c("LUM", "COL6A2", "FN1", "DCN", "BGN", "COL1A2", "COL6A1", "COL3A1", "PCOLCE", "HTRA1", "NNMT", "SPARC", "COL5A2", "MGP", "LRP1", "CCDC80", "COL6A3", "FSTL1", "COL1A1", "SERPING1"),
  Mesenchymal = c("IGF2", "COL1A2", "COL6A2", "MGP", "NNMT", "CXCL12", "DCN", "CYP1B1", "CALD1", "COL6A1", "ANGPTL4", "SERPING1", "IGFBP5", "PCOLCE", "IGFBP4", "C1S", "GGT5", "LUM"),
  Osteoblast = c("IBSP", "BGLAP", "RUNX2", "SATB2", "DLX5", "FAM20C", "CPE", "OMD", "BAMBI", "SERPINF1", "COL8A1", "FAM20A", "TMEM119", "CDH11", "TNC", "CCN2", "CLEC11A", "COL1A1", "COL1A2", "COL12A1"),
  Macro_C1Q = c("C1QA", "C1QB", "C1QC", "TIMD4", "FOLR2", "VSIG4", "CD163", "MS4A7", "MRC1", "SELENOP", "SLC40A1", "MSR1", "LGMN", "NCOA4", "SLC48A1", "LIPA", "IGSF6", "FGL2", "CTSF", "TIMP3", "GPX3", "ADI1"),
  Macro_LYVE1 = c("C1QA", "C1QB", "C1QC", "LYVE1", "FOLR2", "MRC1", "STAB1", "CD163", "VSIG4", "MSR1", "MERTK", "TIMD4", "SELENOP", "F13A1", "RNASE1", "IGF1", "MAFB", "SLC40A1", "LIPA", "PLTP", "TREM2", "LGMN", "FGL2", "ABCA1", "KLF4"),
  Macro_CX3CR1 = c("C1QA", "C1QB", "C1QC", "CX3CR1", "FCGR3A", "MS4A7", "AIF1", "LILRB2", "C5AR1", "SLC11A1", "RHOC", "ITGAX", "CXCL8", "CCL3", "CCL4", "CTSS"),
  Macro_CLEC10A = c("C1QA", "C1QB", "C1QC", "CLEC10A", "ITGAX", "LILRB2", "PLAUR", "PILRA", "CTSH", "SOCS3", "ADGRE5", "FCGR2B", "CCL22", "CD1C", "CPVL", "GPR183", "CSTA"),
  Macro_GPNMB = c("C1QA", "C1QB", "C1QC", "GPNMB", "TREM2", "CD9", "OLR1", "NR1H3", "SCARB2", "ITGA5", "CAPG", "MGLL", "SLC6A6", "FPR3", "ITGAX", "SLC11A1", "CD163", "MSR1", "LGMN"),
  Osteoclast = c("ACP5", "MMP9", "TNFRSF11A", "TCIRG1", "CLCN7", "CTSK", "ITGAV", "SPI1", "ATP6V1A", "ATP6V1B2", "ATP6V1E1", "ATP6V1H", "MATK", "CD109", "SLC37A2", "CLEC4A", "CAPN2", "GNPTAB", "MAP4K4", "OXR1", "SPHK1", "ORAI1", "FOSL2", "SPARC"),
  BMSC_Osteo = c("COL1A1", "COL1A2", "COL16A1", "ITGA10", "ITGA11", "ADAMTS10", "SHOX2", "SULF1", "PTPRS", "MRC2", "UNC5B", "PTPRD", "CDK14", "PLEKHA4", "GXYLT2", "C12orf75", "RUNX2", "SATB2", "DLX5", "FAM20C"),
  BMSC_LEPR = c("LEPR", "CXCL12", "LPL", "APOE", "CFD", "RARRES2", "GGT5", "PTGDS", "PAPPA", "THBS1", "MDK", "APOC1", "ABCA8", "ID4", "ANGPT1", "VCAM1", "FGF7", "CHL1", "CHRDL1"),
  BMSC_THY1 = c("THY1", "ANGPT1", "VCAM1", "FGF7", "TMEM176A", "TMEM176B", "CHRDL1", "CHL1", "GPM6B", "GAS6", "FRZB", "FBN1", "COL14A1", "NTRK2", "S1PR3", "CYP1B1", "LEPR", "CXCL12", "LPL", "APOE"),
  Chondrocyte = c("APOD", "COL2A1", "COL11A1", "SOX9", "SOX5", "FMOD", "OGN", "DCN", "SOD3", "ECRG4", "SERPINA3", "KCNMA1", "PDPN", "BCAT1", "SMOC2", "COL27A1", "PAPSS2", "GPC6", "CRTAC1", "CDO1", "SERPINE2"),
  Fibrochondrocyte = c("PRG4", "CRTAC1", "DPT", "CFB", "GFPT2", "CDO1", "HTRA1", "PROCR", "CLU", "SERPINE2", "CHI3L2", "SLC2A1", "SLC39A14", "INHBA", "FN1", "PDGFRL", "LUM", "ACKR3", "PRRX1"),
  VSMC = c("ACTA2", "TAGLN", "MYL9", "MYLK", "MCAM", "NOTCH3", "EDNRA", "RGS16", "PPP1R14A", "ADIRF", "TINAGL1", "NDUFA4L2", "GPRC5C", "NR2F2", "CSPG4", "HSPB6", "HSPB2", "RBPMS", "TPM2", "PPP1R12A", "SORBS2", "SLIT3", "KANK2", "FILIP1L", "SPARCL1"),
  AEC = c("PECAM1", "CDH5", "CLDN5", "VWF", "SOX18", "CD93", "ESAM", "FLT1", "EFNB2", "EMCN", "ECSCR", "ACVRL1", "ICAM2", "ENG", "JAM2", "S1PR1", "COL15A1", "NES", "BCAM", "ADGRL4"),
  SEC = c("DNASE1L3", "RNASE1", "THBD", "PLAT", "CALCRL", "ADGRL4", "TM4SF1", "TGM2", "GNG11", "BST2", "ITGA6", "CCL14", "CD36", "EGFL7", "CD93", "PECAM1", "CDH5", "ESAM", "CLDN5", "EFNB2"),
  SEC_CD34n = c("CCL14", "STAB1", "FABP4", "RBP7", "TFF3", "EGFL7", "SHANK3", "FAM167B", "CAVIN2", "KDR", "ROBO4", "DNASE1L3", "CD36", "CLDN5", "RNASE1", "PECAM1", "CDH5", "VWF", "FLT1")
)




for (i in 1:length(marker.list.all)) {
  message(i, "  ", marker.list.all[[i]][! marker.list.all[[i]] %in% rownames(sce) ] )
}

cells_rankings <- AUCell_buildRankings(sce@assays$RNA$data, nCores=4)
module_scores <- AUCell_calcAUC(marker.list.all, cells_rankings)


saveRDS(module_scores, paste0(out.path.plot, "/AUC_feature.score.rds") )