空间转录组|Load10X_Spatial函数修改适配多形式数据 + 空转标准流程

生信补给站

发布于 2023-08-25 10:20:11

1K3

发布于 2023-08-25 10:20:11

文章被收录于专栏：生信补给站

在空间转录组|数据读入，标准数据形式外，还有哪些"天残地缺"可以读取提到了多种形式的数据读取，在原函数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 + p2

3，降维聚类及可视化

和单细胞差异不大

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))