本次采用GSE181300数据集，包含了8个头颈部鳞癌患者的数据1、CSV（Comma-Separated Values）：
常用于存储基因表达矩阵的位置信息或细胞/位置的元数据。 例如，tissue_positions_list.csv 文件包含每个条形码的位置坐标、在组织切片上的位置、组织的像素坐标等信息，用于在空间图像上定位每个点。
2、H5（Hierarchical Data Format Version 5）：
通常包含压缩的基因表达矩阵。 filtered_feature_bc_matrix.h5 文件可以包含细胞条形码、基因和对应的表达值。这种格式高效地存储大量数据，便于在 R 和 Python 等程序中快速加载和读取。
3、HTML（HyperText Markup Language）：
用于存储交互式报告或分析结果。 一些空间转录组分析工具会生成 HTML 文件，用于可视化结果、查看细胞类型分布、表达图谱等。通常是报告的一部分，方便研究人员查看结果。
4、JSON（JavaScript Object Notation）：
常用于存储空间缩放因子等配置信息。 scalefactors_json.json 文件包含高分辨率和低分辨率图像的缩放比例，用于将空间坐标映射到图像中的像素坐标。例如，它会定义不同分辨率下组织切片的位置，以保证图像在多分辨率下的准确显示。
5、MTX（Matrix Market Format）： 一种稀疏矩阵格式，通常用于基因表达矩阵。 在空间转录组中，基因表达矩阵可以存储为 .mtx 文件，配合 .barcodes.tsv 和 .genes.tsv 文件，以节省存储空间并快速读取。MTX 文件存储稀疏矩阵的非零元素位置和值。
6、PNG（Portable Network Graphics）：
常用于存储组织切片的图像。 tissue_lowres_image.png 或 tissue_hires_image.png 文件分别包含低分辨率或高分辨率的组织切片图像，用于空间定位和可视化。通常这些图像是空间表达数据的基础，用于将基因表达结果映射到实际组织的结构上。
7、TSV（Tab-Separated Values）：
用于存储条形码和基因信息。 barcodes.tsv 和 genes.tsv 文件分别存储单个位置的条形码以及基因的信息，与 mtx 文件配合使用。这些文件提供了空间转录组中每个位置的标识和相应基因的列表，是空间基因表达数据的重要组成部分。
分析步骤
本次主要是涉及数据读取、部分流程探索及报错解决~
1.导入
rm(list = ls())
library(Seurat)
library(tidyverse)
library(sloop)
## 10x Visium
sce_s <- Load10X_Spatial(data.dir = "./RawData/",
+ filename = "GSM5494475_HNSCC201125T04_filtered_feature_bc_matrix.h5")
出现了如下的报错，笔者把文件下载之后直接丢到了RawData文件中之后进行了解压。这个报错的原因是路劲和命名出现了问题。因此重新整理数据，整理格式如下(供参考)：
把h5格式的文件放在每一个GSM编号的文件下面，然后把figure/json/position文件放在spatial文件夹下边。这里的spatial文件名称是根据报错的结果而命名的(一般都是这样)。
同时需要修改上述三个文件的名称，一般会改成下面图片中的样子。其中图片文件会有高像素和低像素的区分，需要额外注意，这个数据集中只有一种图片文件。
rm(list = ls())
library(Seurat)
library(tidyverse)
library(sloop)
# 10x Visium
sce_s <- Load10X_Spatial(data.dir = "./RawData/GSM5494475",
filename = "GSM5494475_HNSCC201125T04_filtered_feature_bc_matrix.h5")
# # 高像素
# image <- Read10X_Image(image.dir = "./RawData/GSM5494475/spatial/",
# image.name = "tissue_hires_image.png")
# sce_s <- Load10X_Spatial(data.dir = "./RawData/GSM5494475", image = image,
# filename = "GSM5494475_HNSCC201125T04_filtered_feature_bc_matrix.h5")
读进来之后可以check一下数据, 发现比一般的数据多一个images内容
2.check数据
class(sce_s[["slice1"]])
# [1] "VisiumV1" 表示是10X Genomics Visium 空间转录组数据格式
# attr(,"package")
# [1] "Seurat"
head(sce_s[["slice1"]]@coordinates)
# 表示Visium空间转录组数据中每个条形码(spot)的空间坐标信息
# @coordinates包含了每个条形码在组织切片图像上的位置
# tissue row col imagerow imagecol
# AAACACCAATAACTGC-1 1 59 19 1203 3953
# AAACAGCTTTCAGAAG-1 1 43 9 862 3030
# AAACAGGGTCTATATT-1 1 47 13 997 3260
# AAACATGGTGAGAGGA-1 1 62 0 572 4132
# AAACCGGGTAGGTACC-1 1 42 28 1494 2966
# AAACCTCATGAAGTTG-1 1 37 19 1192 2679
# row 和 col：该测序点在空间数据矩阵中的行和列位置。通常用于确定点在原始矩阵中的位置。
# imagerow 和 imagecol：该测序点在实际组织切片图像(即PNG图像)中的像素坐标，分别表示在图像中的行坐标和列坐标。这些坐标用于在组织切片图像上定位测序点的位置，从而将基因表达数据与空间位置关联起来。
class(sce_s[["slice1"]]@image)
# [1] "array"
dim(sce_s[["slice1"]]@image)
# [1] 1972 2000 3
# 1972 和 2000：表示图像的高度(1972像素)和宽度(2000像素)。
# 3：表示图像的颜色通道数，即RGB(红、绿、蓝)三通道。每个像素位置有三个值，分别对应于红色、绿色和蓝色的强度，用于合成彩色图像。
head(GetTissueCoordinates(sce_s))
# imagerow imagecol
# AAACACCAATAACTGC-1 134.38838 441.5937
# AAACAGCTTTCAGAAG-1 96.29492 338.4845
# AAACAGGGTCTATATT-1 111.37591 364.1780
# AAACATGGTGAGAGGA-1 63.89871 461.5900
# AAACCGGGTAGGTACC-1 166.89629 331.3349
# AAACCTCATGAAGTTG-1 133.15956 299.2739
head(GetTissueCoordinates(sce_s[["slice1"]]))
# imagerow imagecol
# AAACACCAATAACTGC-1 134.38838 441.5937
# AAACAGCTTTCAGAAG-1 96.29492 338.4845
# AAACAGGGTCTATATT-1 111.37591 364.1780
# AAACATGGTGAGAGGA-1 63.89871 461.5900
# AAACCGGGTAGGTACC-1 166.89629 331.3349
# AAACCTCATGAAGTTG-1 133.15956 299.2739
我们可以对比一下head(sce_s[["slice1"]]@coordinates)和head(GetTissueCoordinates(sce_s))得到的imagerow/imagecol的坐标不一样。这可能是因为这个数据集的缩放信息存在不匹配导致的。
3.数据预处理
plot1 <- VlnPlot(sce_s, features = "nCount_Spatial", pt.size = 0.2) +
NoLegend()
plot2 <- SpatialFeaturePlot(sce_s, features = "nCount_Spatial",
pt.size.factor = 2) +
theme(legend.position = "right")
wrap_plots(plot1, plot2)
# SCTransform归一化
sce_s <- SCTransform(sce_s, assay = "Spatial",
verbose = FALSE)
SpatialDimPlot(sce_s,alpha = 0)
Plot2中可以看到spot数据在HE切片上的映射，但这里问题来了，为什么看不到HE切片呢？
笔者认为是图片和缩放因子的问题，函数提示需要加载低像素的图片，但是由于没有这种低像素的数据，笔者把高像素的替换了进去。本次先不更换数据，先继续做尝试。
4.基因表达可视化
SpatialFeaturePlot(sce_s, features = c("TP53", "CD3E"),
pt.size.factor = 5)
library(ggplot2)
plot <- SpatialFeaturePlot(sce_s, features = c("TP53"), # alpha = c(0.1, 1)调整透明度
pt.size.factor = 5) + theme(legend.text = element_text(size = 0),
legend.title = element_text(size = 20),
legend.key.size = unit(1, "cm"));plot
ggsave(filename = "./output/images/spatial_tp53.pdf",
height = 9, width = 7)
5.降维聚类可视化
sce_s <- RunPCA(sce_s, assay = "SCT", verbose = FALSE)
sce_s <- FindNeighbors(sce_s, reduction = "pca", dims = 1:30)
sce_s <- FindClusters(sce_s, verbose = FALSE)
sce_s <- RunUMAP(sce_s, reduction = "pca", dims = 1:30)
p1 <- DimPlot(sce_s, reduction = "umap", label = TRUE)
p2 <- SpatialDimPlot(sce_s, label = TRUE, label.size = 3,
pt.size.factor = 5)
p1 + p2
6.提取感兴趣的cluster
SpatialDimPlot(sce_s,
cells.highlight = CellsByIdentities(object = sce_s,
idents = c(2,1,4,3,5,8)),
facet.highlight = TRUE, pt.size.factor = 5,ncol = 3)
7.交互式绘图
这个就不演示了，笔者目前觉得意义不大
# 设置interactive为True
SpatialDimPlot(sce_s, interactive = TRUE)
SpatialFeaturePlot(sce_s, features = "TP53", interactive = TRUE)
8.空间高变特征识别
一、根据组织内预先注释的解剖区域进行差异表达,比如上面通过FindClusters找到的不同cluster
# 指定比较clust5和6
de_markers <- FindMarkers(sce_s, ident.1 = 5, ident.2 = 6)
# 在空间位置上可视化差异表达基因
SpatialFeaturePlot(object = sce_s,
features = rownames(de_markers)[1:3],
alpha = c(0.1, 1), pt.size.factor = 5,
ncol = 3)
二、FindSpatiallyVariables中实现的另一种方法是在没有预注释的情况下搜索表现出空间模式的特征
# 默认method = 'markvariogram', 其他有moransi等
# 差异基因从高变基因中获取
sce_s <- FindSpatiallyVariableFeatures(sce_s,
assay = "SCT",
features = VariableFeatures(sce_s)[1:1000],
selection.method = "moransi")
# 提取6个差异基因看一看
top.features <- head(SpatiallyVariableFeatures(sce_s,selection.method = "moransi"),6)
# 异基因可视化/笔者这里自定义了
top.features <- c("CASQ1", "AEBP1", "MME")
SpatialFeaturePlot(sce_s, features = top.features,
ncol = 3, pt.size.factor = 5, alpha = c(0.1, 1))
今天暂时探索到这儿，后面的是需要对解剖区域(clusters)子集化，涉及到坐标信息的修整，但笔者这里的SPOT信息都没跟HE图片对应起来hhh，所以需要更换数据集之后再进行探索了~
本次内容主要是针对10XVisium数据的读取，此外还有其他类型数据的读取方式等之后再行演示。
参考资料：
10XGENOMICS：https://www.10xgenomics.com/cn Seurat空转：https://satijalab.org/seurat/articles/spatial_vignette.html 生信技能树： https://mp.weixin.qq.com/s/C57QHpg_uD_YkgrmnvRMPA https://mp.weixin.qq.com/s/xVYCwfJI4MUAuka1W99efQ 生信随笔：https://mp.weixin.qq.com/s/X8eP58iDEx4j7pSyGbze9Q
注：若对内容有疑惑或者有发现明确错误的朋友，请联系后台(欢迎交流)。更多内容可关注公众号：生信方舟
- END -