SeekSpace小鼠脑demo数据分析¶

背景信息¶

本⽂档是关于寻因SeekSpace技术所得到的数据及其seurat分析⽅法的说明。
寻因SeekSpace技术可以检测到单细胞精度的基因表达数据,同时也能定位出每个细胞在组织中的空间坐标。
SeekSpace的数据分析简便易学,能够很方便地兼容常见单细胞转录组的分析软件,例如seurat和scanpy。
SeekSpace技术中一个像素点的大小是约为 0.2653 微米,将像素点的坐标乘以 0.2653 即可转换计算细胞在真实空间上的距离。
本文档的输入数据为小鼠脑数据的分析结果。首先使用SeekSpace tools,在默认参数下进行分析;然后,使用seurat进行聚类;进一步,将具有高比例线粒体cluster的细胞过滤掉,过滤后的结果作为本文档输入数据。

⽂件说明¶

SeekSpace技术的配套的基础数据处理的软件是SeekSpace tools,该软件可以从测序文库中识别细胞的表达信息,同时可以识别出每个细胞在空间上的位置。
经过SeekSpace tools软件处理之后所得到的结果文件格式如下:

├── WTH1092_filtered_feature_bc_matrix # 表达矩阵⽬录,可以使⽤seurat的Read10X命令进⾏读取
│ ├── barcodes.tsv.gz
│ ├── features.tsv.gz
│ ├── matrix.mtx.gz
│ └── cell_locations.tsv.gz 为⼩⿏脑测序数据中细胞的空间坐标⽂件。第1列是barcode,顺序与matrix中的filtered_feature_bc_matrix/barcode中细胞的顺序⼀致;第2列和第3列分别为该barcode所代表的细胞的空间位置(即空间芯⽚上的像素坐标)。

├── WTH1092_aligned_DAPI.png # 为⼩⿏脑组织切⽚的DAPI染⾊图⽚
├── WTH1092_aligned_HE_TIMG.png # 为⼩⿏脑组织切⽚的DAPI染⾊图⽚
├── seekspace_of_seurat.ipynb # 为使⽤seurat分析该⼩⿏脑空间数据的jupyter示例⽂件
└── seekspace_of_scanpy.ipynb # 为使⽤scanpy分析该⼩⿏脑空间数据的jupyter示例⽂件
SeekSpace技术中一个像素点的大小是约为 0.2653 微米,将像素点的坐标乘以 0.2653 即可转换计算细胞在真实空间上的距离。

In [1]:
library(Seurat)
library(tidyverse)
library(ggplot2)
Attaching SeuratObject

Seurat v4 was just loaded with SeuratObject v5; disabling v5 assays and
validation routines, and ensuring assays work in strict v3/v4
compatibility mode

── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.0     ✔ stringr   1.5.1
✔ ggplot2   3.5.1     ✔ tibble    3.2.1
✔ lubridate 1.9.3     ✔ tidyr     1.3.1
✔ purrr     1.0.2     
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

创建seurat对象¶

读取SeekSpace的矩阵⽂件,这里直接使用seurat常用的读入和创建矩阵的方法即可。

In [2]:
mouse_brain.data <- Read10X('./Outs/WTH1092_filtered_feature_bc_matrix')
mouse_brain <- CreateSeuratObject(counts=mouse_brain.data,project='mouse_brain')

降维聚类¶

这里使用常见的默认参数进行降维聚类,实际分析的时候,可以按照样品的特征进行相应的修改。

In [3]:
mouse_brain <- NormalizeData(mouse_brain, normalization.method = "LogNormalize", scale.factor = 10000) %>%
FindVariableFeatures(selection.method = "vst", nfeatures = 2000) %>%
ScaleData() %>%
RunPCA() %>%
FindNeighbors(dims = 1:15) %>%
FindClusters(resolution = seq(0.2, 1.4, 0.3)) %>%
RunUMAP(dims = 1:15)
Idents(mouse_brain) <- 'RNA_snn_res.0.8'
Centering and scaling data matrix

PC_ 1 
Positive:  Rbfox1, Kcnip4, Grin2b, Snap25, Tenm2, Grin2a, Nrxn3, Atp1a3, Opcml, Dlgap2 
	   Ppp3r1, Erc2, Atp1b1, Lrrc7, Grm5, Olfm1, Ncdn, Nrgn, Aldoa, Cacna1e 
	   Cntnap2, Atp2b1, Cnksr2, Enc1, Ppp3cb, Ppp3ca, Kctd16, Nell2, Frmpd4, Slit3 
Negative:  Plp1, Trf, Epas1, Glul, Mag, Atp1a2, Cldn11, Mbp, Cnp, Selenop 
	   Mal, Mobp, Ptgds, Flt1, Prr5l, Slc1a3, Gab1, St18, Car2, Fa2h 
	   Apod, Myrf, Cryab, Plpp3, Slco1c1, Mog, Npas3, Kcnj10, Ptprb, Ugt8a 
PC_ 2 
Positive:  Nrgn, 2010300C02Rik, Enc1, Ptk2b, Cnksr2, Wipf3, Cacna1e, Slc8a2, Neurod2, Hpca 
	   Ppp3ca, Tafa1, Slit3, Nptx1, Camk2a, Cpne6, Chrm1, Phactr1, Zbtb20, Egr3 
	   Iqgap2, Nptxr, Slc17a7, Ano3, Rprml, C1ql3, Neurod6, Pde1a, Slit1, Dlgap2 
Negative:  Rbms3, Tcf7l2, Ntng1, Synpo2, Kcnc1, Zfhx3, Kcnc2, Btbd11, Nxph1, Prkcd 
	   Grik1, Grip1, Robo1, Gad1, Pcp4l1, Nap1l5, Gad2, Cit, Inpp4b, Ank1 
	   Slc32a1, Slc17a6, Amotl1, 6330411D24Rik, Alk, Rgs16, Adarb2, Adarb1, Rcan2, Gria4 
PC_ 3 
Positive:  Plp1, Tmeff2, Mag, Pde4b, Trf, Cldn11, St18, Cnp, Fa2h, Mal 
	   Myrf, Enpp2, Mbp, Prr5l, Slc24a2, Mog, Dock10, Cryab, Mobp, Tspan2 
	   Ugt8a, Gpr37, Rnf220, Plekhh1, Pex5l, Ptprd, D7Ertd443e, Phldb1, Gjc3, Gm16168 
Negative:  Flt1, Epas1, Ptprb, Rgs5, Slc1a2, Atp1a2, Utrn, Slco1a4, Fn1, Pltp 
	   Slco1c1, Eng, Ebf1, Tns1, Vtn, Ahnak, Myl9, Igfbp7, Podxl, Adgrl4 
	   Cxcl12, Itga1, Cdh5, Igf2, Adgrf5, Ltbp4, Pitpnc1, Mecom, Pdgfrb, Col4a1 
PC_ 4 
Positive:  Dpp10, Egr1, Ngef, Mef2c, Nrg1, Satb2, Pdzrn3, Kcnh7, Bc1, Car10 
	   Stx1a, Tbr1, Sgcz, Atp1a1, Camk2n1, Homer1, Tmem132d, Lamp5, Efna5, Unc5d 
	   Rgs6, Ipcef1, Slc1a2, Kcnq5, Kcnh5, Atp2b4, Etl4, Garnl3, Snap25, Nrgn 
Negative:  Shisa6, Zbtb20, Nr3c2, Adarb2, Gria1, Tanc1, Galnt17, Prox1, Epha6, Gm20754 
	   St6galnac5, Nkain3, Cacnb2, Slc44a5, Wipf3, Epha7, Sema5a, Lsamp, Iqgap2, Sipa1l3 
	   Grm1, C1ql2, Dscaml1, Synpr, Slit1, Glis3, Ryr3, Pip5k1b, Dpp6, Dgkg 
PC_ 5 
Positive:  Lsamp, Gpc5, Npas3, Gm3764, Gad1, Trpm3, Sox6, Luzp2, Slc6a1, Gad2 
	   Nxph1, Ptprt, Ptprz1, Slc32a1, Slc4a4, Rgs20, Gli2, Erbb4, Plpp3, Grip1 
	   Mertk, Maml2, Slc1a3, Gli3, Egfr, Slc6a11, Nwd1, Reln, Htra1, Kcnmb2 
Negative:  Lef1, Flt1, Ptprb, Rgs5, Synpo2, Epas1, Rnf220, Heg1, Slco1a4, Fn1 
	   Plekhg1, Prkcd, Plp1, Adgrf5, Igfbp7, Kitl, Eng, Pakap.1, Vtn, Ebf1 
	   Ltbp4, Pltp, Myl9, Cxcl12, Podxl, Jcad, Tns1, Adgrl4, Cdh5, Amotl1 

Computing nearest neighbor graph

Computing SNN

Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 22870
Number of edges: 834752

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9677
Number of communities: 15
Elapsed time: 3 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 22870
Number of edges: 834752

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9429
Number of communities: 17
Elapsed time: 3 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 22870
Number of edges: 834752

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9234
Number of communities: 22
Elapsed time: 3 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 22870
Number of edges: 834752

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.9076
Number of communities: 26
Elapsed time: 3 seconds
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 22870
Number of edges: 834752

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.8952
Number of communities: 32
Elapsed time: 3 seconds
Warning message:
“The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
This message will be shown once per session”
10:46:03 UMAP embedding parameters a = 0.9922 b = 1.112

10:46:03 Read 22870 rows and found 15 numeric columns

10:46:03 Using Annoy for neighbor search, n_neighbors = 30

10:46:03 Building Annoy index with metric = cosine, n_trees = 50

0%   10   20   30   40   50   60   70   80   90   100%

[----|----|----|----|----|----|----|----|----|----|

*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
*
|

10:46:06 Writing NN index file to temp file /tmp/RtmpCReRpJ/file10268882ade

10:46:06 Searching Annoy index using 1 thread, search_k = 3000

10:46:15 Annoy recall = 100%

10:46:15 Commencing smooth kNN distance calibration using 1 thread
 with target n_neighbors = 30

10:46:16 Initializing from normalized Laplacian + noise (using irlba)

10:46:19 Commencing optimization for 200 epochs, with 1000636 positive edges

10:46:33 Optimization finished

添加空间坐标¶

接下来我们使用seurat的CreateDimReducObject函数,给seurat对象中的每个细胞添加一个空间的坐标。
注意!!
CreateDimReducObject函数读入空间坐标矩阵的细胞的顺序必须与mouse_brain里面的meta.data里面细胞的顺序相一致,否则在多样品整合分析的时候会出现细胞顺序错乱的情况。
这里我们使用了一个简单的一行代码对空间坐标矩阵进行了排序。

In [4]:
spatial_df <- read.table('./Outs/WTH1092_filtered_feature_bc_matrix/cell_locations.tsv.gz', row.names = 1, sep = '\t',header = T)
colnames(spatial_df) <- c("spatial_1","spatial_2")
spatial_matrix <- as.matrix(spatial_df)
spatial_matrix_sorted <- spatial_matrix[match(row.names(mouse_brain@meta.data),row.names(spatial_matrix)), ]
mouse_brain@reductions$spatial <- CreateDimReducObject(embeddings = spatial_matrix_sorted, key='spatial_', assay='RNA')
In [5]:
head(mouse_brain@reductions$spatial@cell.embeddings)
A matrix: 6 × 2 of type int
spatial_1spatial_2
AAGTCGTCGACTCTAGAGCTCTGCTTC3448017330
ACAGGGATACTTCCAGTCCTAATAACG3163412314
ACCACAAATGGCGTTTCGCTCTGCTTC1340716626
ACTATGGTACAACTCCTCTAAGTACGA29719 4091
AGACTCAATGTTCTGGCTAGTCATACG32670 4484
AGCCAATCGAACGGGATCCTAATAACG1387415676

经过这步处理之后,我们可以在seurat对象中看到一个新的坐标系”spatial“,这个坐标系即为每个细胞在空间位置上的坐标。

In [6]:
str(mouse_brain@reductions$spatial)
Formal class 'DimReduc' [package "SeuratObject"] with 9 slots
  ..@ cell.embeddings           : int [1:22870, 1:2] 34480 31634 13407 29719 32670 13874 19186 18630 22301 37179 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:22870] "AAGTCGTCGACTCTAGAGCTCTGCTTC" "ACAGGGATACTTCCAGTCCTAATAACG" "ACCACAAATGGCGTTTCGCTCTGCTTC" "ACTATGGTACAACTCCTCTAAGTACGA" ...
  .. .. ..$ : chr [1:2] "spatial_1" "spatial_2"
  ..@ feature.loadings          : num[0 , 0 ] 
  ..@ feature.loadings.projected: num[0 , 0 ] 
  ..@ assay.used                : chr "RNA"
  ..@ global                    : logi FALSE
  ..@ stdev                     : num(0) 
  ..@ jackstraw                 :Formal class 'JackStrawData' [package "SeuratObject"] with 4 slots
  .. .. ..@ empirical.p.values     : num[0 , 0 ] 
  .. .. ..@ fake.reduction.scores  : num[0 , 0 ] 
  .. .. ..@ empirical.p.values.full: num[0 , 0 ] 
  .. .. ..@ overall.p.values       : num[0 , 0 ] 
  ..@ misc                      : list()
  ..@ key                       : chr "spatial_"

添加组织图像信息¶

下面我们展示一下在空间数据上添加背景图片。

In [7]:
samplename = 'WTH1092'
size_x = 55128
size_y = 19906
#如果芯片编号的倒数第二个字母是A,则size_x=55128,size_y=19906;
#如果芯片编号的倒数第二个字母是B,则size_x=55050,size_y=19906;
#如果芯片编号的倒数第二个字母是C,则size_x=55050,size_y=19906.
#(可以在samplename_report.html质控报告的Summary Table中查看芯片编号)
mouse_brain@misc$info[[`samplename`]]$size_x = as.integer(size_x)
mouse_brain@misc$info[[`samplename`]]$size_y = as.integer(size_y)

添加DAPI图片¶

In [8]:
img = './Outs/WTH1092_aligned_DAPI.png'  

#base64格式
img_64 = base64enc::dataURI(file = img)
mouse_brain@misc$info[[`samplename`]]$img = img_64
#png格式
img_gg <- png::readPNG(img)
img_grob <- grid::rasterGrob(img_gg, interpolate = FALSE, width = grid::unit(1,"npc"), height = grid::unit(1, "npc"))
mouse_brain@misc$info[[`samplename`]]$img_gg = img_grob

添加HE图片¶

In [9]:
img = './Outs/WTH1092_aligned_HE_TIMG.png'  

#base64格式
img_64 = base64enc::dataURI(file = img)
mouse_brain@misc$info[[`samplename`]]$img_he = img_64
#png格式
img_gg <- png::readPNG(img)
img_grob <- grid::rasterGrob(img_gg, interpolate = FALSE, width = grid::unit(1,"npc"), height = grid::unit(1, "npc"))
mouse_brain@misc$info[[`samplename`]]$img_he_gg = img_grob
In [10]:
str(mouse_brain@misc$info)
List of 1
 $ WTH1092:List of 6
  ..$ size_x   : int 55128
  ..$ size_y   : int 19906
  ..$ img      : chr "data:;base64,iVBORw0KGgoAAAANSUhEUgAAD6AAAAWnCAAAAAByYdMoAAAgAElEQVR4AezBaa9l6WGe5/t537X22vOZ6lSdququrq5ukt2kRI"| __truncated__
  ..$ img_gg   :List of 12
  .. ..$ raster     : 'raster' chr [1:1447, 1:4000] "#000000" "#000000" "#000000" "#000000" ...
  .. ..$ x          : 'simpleUnit' num 0.5npc
  .. .. ..- attr(*, "unit")= int 0
  .. ..$ y          : 'simpleUnit' num 0.5npc
  .. .. ..- attr(*, "unit")= int 0
  .. ..$ width      : 'simpleUnit' num 1npc
  .. .. ..- attr(*, "unit")= int 0
  .. ..$ height     : 'simpleUnit' num 1npc
  .. .. ..- attr(*, "unit")= int 0
  .. ..$ just       : chr "centre"
  .. ..$ hjust      : NULL
  .. ..$ vjust      : NULL
  .. ..$ interpolate: logi FALSE
  .. ..$ name       : chr "GRID.rastergrob.11"
  .. ..$ gp         : list()
  .. .. ..- attr(*, "class")= chr "gpar"
  .. ..$ vp         : NULL
  .. ..- attr(*, "class")= chr [1:3] "rastergrob" "grob" "gDesc"
  ..$ img_he   : chr "data:;base64,iVBORw0KGgoAAAANSUhEUgAAD6AAAAWnCAIAAADYaBujAAAgAElEQVR4AZThUZYlx4JlR8rW4Jy5+NvzJBcHkq6HZtcdQACJV5"| __truncated__
  ..$ img_he_gg:List of 12
  .. ..$ raster     : 'raster' chr [1:1447, 1:4000] "#DEDDE0" "#777982" "#74757F" "#727179" ...
  .. ..$ x          : 'simpleUnit' num 0.5npc
  .. .. ..- attr(*, "unit")= int 0
  .. ..$ y          : 'simpleUnit' num 0.5npc
  .. .. ..- attr(*, "unit")= int 0
  .. ..$ width      : 'simpleUnit' num 1npc
  .. .. ..- attr(*, "unit")= int 0
  .. ..$ height     : 'simpleUnit' num 1npc
  .. .. ..- attr(*, "unit")= int 0
  .. ..$ just       : chr "centre"
  .. ..$ hjust      : NULL
  .. ..$ vjust      : NULL
  .. ..$ interpolate: logi FALSE
  .. ..$ name       : chr "GRID.rastergrob.12"
  .. ..$ gp         : list()
  .. .. ..- attr(*, "class")= chr "gpar"
  .. ..$ vp         : NULL
  .. ..- attr(*, "class")= chr [1:3] "rastergrob" "grob" "gDesc"

细胞注释结果¶

接下来我们再把细胞的注释结果添加到每个细胞上,这样我们就可以看不同的细胞类型在空间上的分布情况了。
SeekSpace细胞注释的过程,跟普通的单细胞注释的过程完全一致。
我们之前已经对demo数据进行了大群和亚群的注释,注释的结果文件放在了anntation.csv文件中。
接下来我们简单的读入处理一下。

In [11]:
anno <- read.csv("./annotation.csv",row.names = 1,header = TRUE)
In [12]:
mouse_brain <- AddMetaData(mouse_brain, anno)
In [13]:
head(mouse_brain)
A data.frame: 10 × 12
orig.identnCount_RNAnFeature_RNARNA_snn_res.0.2RNA_snn_res.0.5RNA_snn_res.0.8RNA_snn_res.1.1RNA_snn_res.1.4seurat_clusterssc_barcodeMain_CellTypeSub_CellType
<fct><dbl><int><fct><fct><fct><fct><fct><fct><chr><chr><chr>
AAGTCGTCGACTCTAGAGCTCTGCTTCmouse_brain200156111215151717AAGTCGTCGACTCTAGAGCTCTGCTTCOPC OPC
ACAGGGATACTTCCAGTCCTAATAACGmouse_brain2001463 3 1 0 2 2 ACAGGGATACTTCCAGTCCTAATAACGAstroAstro
ACCACAAATGGCGTTTCGCTCTGCTTCmouse_brain200145141621242929ACCACAAATGGCGTTTCGCTCTGCTTCInh Ependymal
ACTATGGTACAACTCCTCTAAGTACGAmouse_brain2001571 1 117 8 8 ACTATGGTACAACTCCTCTAAGTACGAOligoOligo
AGACTCAATGTTCTGGCTAGTCATACGmouse_brain2001204 1518162323AGACTCAATGTTCTGGCTAGTCATACGExt nb
AGCCAATCGAACGGGATCCTAATAACGmouse_brain2001351 1 3 7 5 5 AGCCAATCGAACGGGATCCTAATAACGOligoOligo
AGGTTGTCGAGCGTATTGCTCTGCTTCmouse_brain2001503 3 1 0 2 2 AGGTTGTCGAGCGTATTGCTCTGCTTCAstroAstro
AGTAGCTCGACATCGCGGGACTGTGGAmouse_brain2001404 5 2 2 7 7 AGTAGCTCGACATCGCGGGACTGTGGAEndo Endo
ATGGATCTACTACGGTAAGCGACGGATmouse_brain2001261 1 3 7 8 8 ATGGATCTACTACGGTAAGCGACGGATOligoOligo
CCTAAGAGCTACCGACCAAGGCACTATmouse_brain2001204 5 2 2 7 7 CCTAAGAGCTACCGACCAAGGCACTATEndo Endo
In [14]:
saveRDS(mouse_brain,"WTH1092_demo_mouse_brain.rds")

结果可视化¶

经过seurat处理后的数据保存为rds,我们可以继续使用seurat包默认的一系列函数进行绘图

空间坐标绘图¶

In [15]:
#dimplot umap
options(repr.plot.height=7, repr.plot.width=7)
DimPlot(mouse_brain, reduction = 'umap')
No description has been provided for this image
In [16]:
options(repr.plot.height=10, repr.plot.width=10)
FeaturePlot(mouse_brain, reduction = 'umap', features=c('Mbp','Mobp', 'Olig1','Plp1'), pt.size = 1)
No description has been provided for this image

当我们把 DimPlot 函数中的 reduction 参数设置为 "spatial" 的时候,我们即可将细胞绘制在空间水平上。

In [17]:
options(repr.plot.height=7, repr.plot.width=15)
DimPlot(mouse_brain, reduction = 'spatial',pt.size = 1)
No description has been provided for this image
In [18]:
options(repr.plot.height=7, repr.plot.width=15)
DimPlot(mouse_brain, reduction = 'spatial',pt.size = 1,group.by = "Sub_CellType")
No description has been provided for this image

同样,我们可以在使用FeaturePlot 函数时,将reduction参数设置为"spatial",来查看基因在空间位置上的表达情况。

In [19]:
options(repr.plot.height=10, repr.plot.width=20)
FeaturePlot(mouse_brain, reduction = 'spatial', features = c('Mbp','Mobp', 'Olig1','Plp1'),pt.size = 0.7)
No description has been provided for this image

SeekSpace数据对Seurat的兼容性比较好,取子集等等函数均可运行。

In [20]:
options(repr.plot.height=7, repr.plot.width=15)
sub_mouse_brain <- subset(mouse_brain, Main_CellType == 'Ext')
DimPlot(sub_mouse_brain, reduction = 'spatial',pt.size = 1.5, group.by = "Sub_CellType")
No description has been provided for this image

带有DAPI/HE的图片绘制¶

自定义了ImageSpacePlot和FeatureSpacePlot两个函数绘制细胞在空间上的分群信息和细胞连续变量的指标

In [21]:
######################绘图细胞分群函数
ImageSpacePlot = function(obj, group_by, type="DAPI", sample=names(obj@misc$info)[1], size=1, alpha=1,color=MYCOLOR){
    MYCOLOR=c(
          "#6394ce", "#2a4c87", "#eed500", "#ed5858",
          "#f6cbc2", "#f5a2a2", "#3ca676", "#6cc9d8",
          "#ef4db0", "#992269", "#bcb34a", "#74acf3",
          "#3e275b", "#fbec7e", "#ec4d3d", "#ee807e",
          "#f7bdb5", "#dbdde6", "#f591e1", "#51678c",
          "#2fbcd3", "#80cfc3", "#fbefd1", "#edb8b5",
          "#5678a8", "#2fb290", "#a6b5cd", "#90d1c1",
          "#a4e0ea", "#837fd3", "#5dce8b", "#c5cdd9",
          "#f9e2d6", "#c64ea4", "#b2dfd6", "#dbdfe7",
          "#dff2ec", "#cce8f3", "#e74d51", "#f7c9c4",
          "#f29c81", "#c9e6e0", "#c1c5de", "#750000"
          )
        
        raster_type <- switch(type,
                          HE = "img_he_gg",
                          DAPI = "img_gg",
                          stop("Invalid type. Must be 'HE' or 'DAPI'.")
                             )

        spatial_coord1 <- as.data.frame(obj[[group_by]])
        colnames(spatial_coord1) <- group_by
        spatial_coord2 <- as.data.frame(obj@reductions$spatial@cell.embeddings)
        spatial_coord <-cbind(spatial_coord2,spatial_coord1)
 
    
        ImageSpacePlot <- ggplot2::ggplot() + ggplot2::annotation_custom(grob = obj@misc$info[[sample]][[raster_type]],
        xmin = 0, xmax = obj@misc$info[[sample]]$size_x, 
        ymin = 0, ymax = obj@misc$info[[sample]]$size_y) +
        ggplot2::geom_point(data = spatial_coord, ggplot2::aes(x = spatial_1,y = spatial_2, color = !!sym(group_by), 
                            fill = !!sym(group_by)), size=size, alpha=alpha)+
        labs(size = group_by) + guides(alpha = "none")+ 
        ggplot2::theme_classic()+
        scale_color_manual(values = color)+ coord_fixed()
    return(ImageSpacePlot)
}

###################绘图基因表达函数
FeatureSpacePlot = function(obj, feature, type="DAPI", sample=names(obj@misc$info)[1], size=1, alpha=c(1,1),color=c("lightgrey","blue")){
    raster_type <- switch(type,
                          HE = "img_he_gg",
                          DAPI = "img_gg",
                          stop("Invalid type. Must be 'HE' or 'DAPI'.")
                             )

    spatial_coord1 <- as.data.frame(obj@reductions$spatial@cell.embeddings)
    spatial_coord2 <- FetchData(obj,feature)
    colnames(spatial_coord2) <- feature
    spatial_coord <-cbind(spatial_coord1,spatial_coord2)

    FeatureSpacePlot <-ggplot2::ggplot() + ggplot2::annotation_custom(grob = obj@misc$info[[sample]][[raster_type]],
        xmin = 0, xmax = obj@misc$info[[sample]]$size_x, ymin = 0, ymax = obj@misc$info[[sample]]$size_y) +
        ggplot2::geom_point(data = spatial_coord, ggplot2::aes(x = spatial_1, y = spatial_2,color = !!sym(feature),alpha = !!sym(feature)), size=size)+
        labs(color = feature)+
        guides(alpha = "none")+
        ggplot2::theme_classic()+
        ggplot2::scale_alpha_continuous(range=alpha)+
        scale_color_gradient(low=color[1],high = color[2])+ coord_fixed()
    return(FeatureSpacePlot)
    }
In [22]:
# DAPI背景的细胞分群图
options(repr.plot.height=7, repr.plot.width=15)
ImageSpacePlot(obj=mouse_brain, group_by = "Sub_CellType",type="DAPI",size=0.7)
No description has been provided for this image
In [23]:
# HE背景的细胞分群图
options(repr.plot.height=7, repr.plot.width=15)
ImageSpacePlot(obj=mouse_brain, group_by = "Sub_CellType",type="HE")
No description has been provided for this image
In [24]:
# DAPI背景的基因表达图
options(repr.plot.height=7, repr.plot.width=15)
FeatureSpacePlot(obj=mouse_brain, feature="Hpca",type="DAPI")
No description has been provided for this image
In [25]:
# HE背景的基因表达图
options(repr.plot.height=7, repr.plot.width=15)
FeatureSpacePlot(obj=mouse_brain, feature="Hpca",type="HE")
No description has been provided for this image
In [ ]: