空间转录组|Load10X_Spatial函数修改适配多形式数据 + 空转标准流程
空间转录组|Load10X_Spatial函数修改适配多形式数据 + 空转标准流程
在空间转录组|数据读入，标准数据形式外，还有哪些"天残地缺"可以读取提到了多种形式的数据读取，在原函数Load10X_Spatial的基础上进行了 “简陋”的修改，添加一些判断使之可以读取上文提到的多种数据形式。
后台回复 “函数” 既可以获取Load10X_Spatial_change.R 和 测试数据文件。第二部分为空转的Seurat标准分析流程，重点介绍了和单细胞转录组的差异。
一 Load10X_Spatial_change函数测试
加载R包和函数，然后使用不用场景进行测试。
library(Seurat)
library(jsonlite)
library(png)
library(tidyverse)
library(ggpubr)
library(patchwork)
source("Load10X_Spatial_change.R")1, 标准数据读取
含有h5文件 + spatial文件夹
T0 <- Load10X_Spatial_change(data.dir = ".\\outs\\",filename = "filtered_feature_bc_matrix.h5")
for (i in colnames((T0@images$slice1@coordinates))) {
T0@images$slice1@coordinates[[i]] <- as.integer(T0@images$slice1@coordinates[[i]])
}
SpatialDimPlot(T0,alpha = 0)注：循环修改integer 的操作，也许你的数据不需要
2，filtered_feature_bc_matrix文件夹 + spatial文件夹（无 h5文件）
data.dir = ".\\outs_noh5\\"
T1 <- Load10X_Spatial_change(data.dir = data.dir,
filename = "filtered_feature_bc_matrix")
for (i in colnames((T1@images$slice1@coordinates))) {
T1@images$slice1@coordinates[[i]] <- as.integer(T1@images$slice1@coordinates[[i]])
}
p1 <- SpatialDimPlot(T1,alpha = 0)
p2 <- SpatialFeaturePlot(T1, features = "nFeature_Spatial",
pt.size.factor = 3)
p1 + p2
3, h5 + spatial文件夹(只有tissue_hires_image.png )
通过image.name指定使用的图片
data.dir = ".\\outs_onlyhires_h5\\"
T2 <- Load10X_Spatial_change(data.dir = data.dir,
filename = "filtered_feature_bc_matrix.h5",
#slice = "sli",
image.name = "tissue_hires_image.png")
for (i in colnames((T2@images[[1]]@coordinates))) {
T2@images[[1]]@coordinates[[i]] <- as.integer(T2@images[[1]]@coordinates[[i]])
}
SpatialDimPlot(T2,alpha = 0)
SpatialFeaturePlot(T2, features = "nFeature_Spatial",
pt.size.factor = 8)4, filtered_feature_bc_matrix + spatial文件夹(只有tissue_hires_image.png )
data.dir = ".\\outs_onlyhires\\"
T3 <- Load10X_Spatial_change(data.dir = data.dir,
filename = "filtered_feature_bc_matrix",
image.name = "tissue_hires_image.png")
for (i in colnames((T3@images$slice1@coordinates))) {
T3@images$slice1@coordinates[[i]] <- as.integer(T3@images$slice1@coordinates[[i]])
}
SpatialDimPlot(T3,alpha = 0)
SpatialFeaturePlot(T3, features = "nFeature_Spatial",
pt.size.factor = 8)二 空转标准分析
空转流程和单细胞转录组的流程有一些细微的区别，重点介绍这部分
1，数据预处理（过滤，NormalizeData）
关于空转数据是否进行质控过滤，以及根据哪些条件进行过滤差异“很大”，且都有较多文章支撑，这里仅列出过滤基因 以及 spot 的方式，是否过滤根据自己的数据情况确定。
注：Seurat的官网并未进行质控，直接SCT分析。
1）计算线粒体/核糖体/特定基因集的比例
mt.genes <- grep(pattern = "^MT-", x = rownames(T0), value = TRUE)
T0$percent.mito <- (Matrix::colSums(T0@assays$Spatial@counts[mt.genes, ])/Matrix::colSums(T0@assays$Spatial@counts))*100如果想针对核糖体或者特定基因集合，只需进行对应的改即可。
2）过滤/保留基因
保留在所有spot中的表达量之和大于5的基因
genes_to_keep <- setdiff(names(which(Matrix::rowSums(T0@assays$Spatial@counts )>5)),mt.genes)也可以是自定义的基因集合。
3）subset过滤spot
T0 <- subset(T0,
features =genes_to_keep, #针对基因
subset = nFeature_Spatial > 300 & percent.mito < 30 #针对spots
)
plot1 <- VlnPlot(T0, features = "nCount_Spatial", pt.size = 0.1) + NoLegend()
plot2 <- SpatialFeaturePlot(T0, features = "nCount_Spatial",pt.size.factor = 3) +
theme(legend.position = "right")
wrap_plots(plot1, plot2)
NormalizeData
相较于单细胞转录组LogNormalize，空转数据推荐使用SCTransform。下图是来源于Seurat官网（https://satijalab.org/seurat/articles/seurat5_spatial_vignette.html#gene-expression-visualization）的比较结果。
T0 <- SCTransform(T0, assay = "Spatial", verbose = FALSE)
2，基因可视化
空转的可视化函数扩充了SpatialFeaturePlot 和 SpatialDimPlot （以空转切片为背景） ，相较于单细胞的FeaturePlot，DimPlot函数。
通过pt.size.factor设置点的大小，仍然可以使用ggplot2进行一些图形的修饰（theme，legend等）。
SpatialFeaturePlot(T0, features = c("COL1A1", "EPCAM"),
pt.size.factor = 4)
library(ggplot2)
SpatialFeaturePlot(T0, features = c("COL1A1", "EPCAM"),pt.size.factor = 3) +
theme(legend.text = element_text(size = 0),
legend.title = element_text(size = 20),
legend.key.size = unit(1, "cm"))
p1 <- SpatialFeaturePlot(T0, features = "COL1A1", pt.size.factor = 1)
p2 <- SpatialFeaturePlot(T0, features = "COL1A1", alpha = c(0.1, 1))
p1 + p23，降维聚类及可视化
和单细胞差异不大
T0 <- RunPCA(T0, assay = "SCT", verbose = FALSE)
T0 <- FindNeighbors(T0, reduction = "pca", dims = 1:30)
T0 <- FindClusters(T0, verbose = FALSE)
T0 <- RunUMAP(T0, reduction = "pca", dims = 1:30)
p1 <- DimPlot(T0, reduction = "umap", label = TRUE)
p2 <- SpatialDimPlot(T0, label = TRUE, label.size = 3)
p1 + p2
可以使用SpatialDimPlot 的cells.highlight参数， 只展示感兴趣的簇。
SpatialDimPlot(T0,
cells.highlight = CellsByIdentities(object = T0,
idents = c(2, 0, 3)),
facet.highlight = TRUE, ncol = 3 ,pt.size.factor = 3 )
4，空间表达变量基因
由于空转的特殊性，Seurat提供了两种计算的方式 。
1）FindMarkers
根据定义好的cluster或者celltype进行计算（和scRNA类似）
#指定比较的group
de_markers <- FindMarkers(T0, ident.1 = 0, ident.2 = 3)
SpatialFeaturePlot(object = T0, features = rownames(de_markers)[1:3],
alpha = c(0.1, 1), pt.size.factor = 3,
ncol = 3)
#全部比较，注意activate的ident
de_markers2 <- FindAllMarkers(T0)
注：FindAllMarkers函数注意当前activate的ident
2）FindSpatiallyVariableFeatures
在没有cluster 或者 celltype等预先注释的情况下使用，倾向于返回在某些切片区域高表达的基因。
T0 <- FindSpatiallyVariableFeatures(T0, assay = "SCT", features = VariableFeatures(T0)[1:1000],
selection.method = "moransi")
top.features <- head(SpatiallyVariableFeatures(T0, selection.method = "moransi"), 6)
SpatialFeaturePlot(T0, features = top.features,
pt.size.factor = 3,
ncol = 3, alpha = c(0.1, 1))
两种之间会有一些overlap，但是也会有一些差异。