STUtility || 空间转录组多样本分析框架（一）
男，
一个长大了才会遇到的帅哥，
稳健，潇洒，大方，靠谱。
一段生信缘，一棵技能树。
生信技能树核心成员，单细胞天地特约撰稿人，简书创作者，单细胞数据科学家。
空间转录组学是一种通过结合基因表达数据和显微图像数据来可视化和定量分析组织切片中转录组定量的方法。在前几期的文章中，我们主要讲述的是单个空间转录组样本的分析，今天要讲的是如何分析多张切片。处理过单细胞转录组的同学不会忘记，多样本分析和单样本是多么的不同。在空间这里关键的一点是多个图像的处理（对齐）。STUtility的开发者Ludvig Larsson和Joseph Bergenstrahle是Joakim Lundebergs教授团队的博士生，该团队是空间转录组技术(ST)的最初发明者，后来被10X Genomics收购。团队工作室位于瑞典斯德哥尔摩的生命科学实验室(SciLifeLab)。为了给大家一个宏观的视角，来看看人家的实验室是怎样的：
开发STUtility的目标是分析多个空转切片。与所有生物数据一样，使用多个样本增加了分析的力量，同时减少不确定性的影响。这里，我们将展示如何输入多个切片，以期查看诸如批次效应和缺失数据等情况以及应用已有的方法纠正这些情况，获得数据的整体印象，并以各种方式进行更深入的分析。
我们已经广泛尝试了处理ST（空转）数据的不同方法和工作流程。虽然条条大路通罗马，但截止到撰写本文时，我们发现Seurat方法最适合这种类型的数据。Seurat是一个专为单细胞RNAseq数据设计的R包。显然，这偏离了ST技术目前产生的数据，因为阵列上的分辨率意味着每个捕获点由源自多个细胞的转录本组成。然而，ST数据的特征在很大程度上与scRNAseq相似（我们会在下一篇文章中用数据说明之）。注意，STUtility包依赖Seurat v3.0或更高版本。
本节主要内容有：
多切片数据读取 空转数据质控 图像对齐（旋转，切割） 绘制感兴趣的区域 空间数据3D可视化
我们载入R包和10X的空转数据：
library(STutility)
library(SeuratData)
library(Seurat)
library(tidyverse)
df <- data.frame(samples =c("E:\\learnscanpy\\data\\V1_Breast_Cancer_Block_A_Section_1_spatial\\filtered_feature_bc_matrix.h5",
'E:\\learnscanpy\\data\\V1_Breast_Cancer_Block_A_Section_2\\filtered_feature_bc_matrix.h5'),
spotfiles = c("E:\\learnscanpy\\data\\V1_Breast_Cancer_Block_A_Section_1_spatial\\spatial\\tissue_positions_list.csv",
'E:\\learnscanpy\\data\\V1_Breast_Cancer_Block_A_Section_2\\spatial\\tissue_positions_list.csv'),
imgs = c("E:\\learnscanpy\\data\\V1_Breast_Cancer_Block_A_Section_1_spatial\\spatial\\tissue_hires_image.png",
'E:\\learnscanpy\\data\\V1_Breast_Cancer_Block_A_Section_2\\spatial\\tissue_hires_image.png'),
json = c("E:\\learnscanpy\\data\\V1_Breast_Cancer_Block_A_Section_1_spatial\\spatial\\scalefactors_json.json",
'E:\\learnscanpy\\data\\V1_Breast_Cancer_Block_A_Section_2\\spatial\\scalefactors_json.json'))
df
和Seurat的思想很像，其实很多表达量分析的函数也是直接借用Seurat的。
?InputFromTable # 依然会过滤一些
se <- InputFromTable(infotable = df,
min.gene.count = 100,
min.gene.spots = 5,
min.spot.count = 500,
platform = "Visium")
查看数据对象结构，俨然一个Seurat对象：
se
An object of class Seurat
16273 features across 7798 samples within 1 assay
Active assay: RNA (16273 features, 0 variable features)
> head(se@meta.data)
orig.ident nCount_RNA nFeature_RNA id labels
CAGGATCCGCCCGACC-1_1 SeuratProject 5406 2555 1 Default
CACGATTGGTCGTTAA-1_1 SeuratProject 5740 2622 2 Default
GGTTGTATCGTGAAAT-1_1 SeuratProject 7272 3165 3 Default
TCTTATGGGTAGTACC-1_1 SeuratProject 4208 2081 4 Default
TACAAGCTGTTCACTG-1_1 SeuratProject 8111 3276 5 Default
GTATCTTGTTGCTCAC-1_1 SeuratProject 7407 3118 6 Default
head(se@tools$Staffli@meta.data)
x y adj_x adj_y pixel_x pixel_y sample original_x original_y
CAGGATCCGCCCGACC-1_1 49 1 49 1 842.6568 353.3828 1 842.6568 353.3828
CACGATTGGTCGTTAA-1_1 50 0 50 0 853.9604 333.8284 1 853.9604 333.8284
GGTTGTATCGTGAAAT-1_1 51 1 51 1 865.1815 353.4653 1 865.1815 353.4653
TCTTATGGGTAGTACC-1_1 52 0 52 0 876.4851 333.9109 1 876.4851 333.9109
TACAAGCTGTTCACTG-1_1 53 1 53 1 887.7063 353.4653 1 887.7063 353.4653
GTATCTTGTTGCTCAC-1_1 54 0 54 0 899.0099 333.9109 1 899.0099 333.9109
查看数据分布：
library(foreach)
library(parallel)
plt <- function(i){
print(ggplot() +
geom_histogram(data = se[[]], aes(get(i)), fill = "red", alpha = 0.7, color = "gray", bins = 50) +
Seurat::DarkTheme() +
ggtitle(paste0("Total", i , " per spots")))
}
pl = list()
pl1 <- foreach::foreach(i = c('nFeature_RNA','nCount_RNA'),.packages = c("Seurat","ggplot2")) %dopar% plt(i)
gene_attr <- data.frame(nUMI = Matrix::rowSums(se@assays$RNA@counts),
nSpots = Matrix::rowSums(se@assays$RNA@counts > 0))
plt2 <- function(i){
print(ggplot() +
geom_histogram(data = gene_attr, aes(get(i)), fill = "red", alpha = 0.7, color = "gray", bins = 50) +
Seurat::DarkTheme() +
ggtitle(paste0("Total", i , " per gene")))
}
pl2 <- foreach::foreach(i = c('nUMI','nSpots'),.packages = c("Seurat","ggplot2")) %dopar% plt2(i)
cowplot::plot_grid(plotlist = c(pl1,pl2))
线粒体核糖体基因等质量指标：
# Collect all genes coded on the mitochondrial genome
mt.genes <- grep(pattern = "^MT-", x = rownames(se), value = TRUE)
se$percent.mito <- (Matrix::colSums(se@assays$RNA@counts[mt.genes, ])/Matrix::colSums(se@assays$RNA@counts))*100
# Collect all genes coding for ribosomal proteins
rp.genes <- grep(pattern = "^RPL|^RPS", x = rownames(se), value = TRUE)
se$percent.ribo <- (Matrix::colSums(se@assays$RNA@counts[rp.genes, ])/Matrix::colSums(se@assays$RNA@counts))*100
head(se@meta.data)
orig.ident nCount_RNA nFeature_RNA id labels percent.mito percent.ribo
CAGGATCCGCCCGACC-1_1 SeuratProject 5406 2555 1 Default 9.507954 10.229375
CACGATTGGTCGTTAA-1_1 SeuratProject 5740 2622 2 Default 8.832753 11.672474
GGTTGTATCGTGAAAT-1_1 SeuratProject 7272 3165 3 Default 8.044554 10.588559
TCTTATGGGTAGTACC-1_1 SeuratProject 4208 2081 4 Default 8.911597 11.644487
TACAAGCTGTTCACTG-1_1 SeuratProject 8111 3276 5 Default 10.393293 9.517939
GTATCTTGTTGCTCAC-1_1 SeuratProject 7407 3118 6 Default 9.693533 10.476576
p1 <- map(c('percent.mito','percent.ribo') ,function(x) ST.FeaturePlot(se, features = x, dark.theme = TRUE, cols = c("dark blue", "cyan", "yellow", "red", "dark red")))
cowplot::plot_grid(plotlist = c(p1))
这里来了一个分析空转数据值得思考的问题：在空间上去掉一个Sopt意味着什么？
# Keep spots with more than 500 unique genes and less than 10% mitochondrial content
se.subset <- SubsetSTData(se, expression = nFeature_RNA > 500 & nCount_RNA < 20000 & percent.mito < 10)
cat("Spots removed: ", ncol(se) - ncol(se.subset), "\n")
Spots removed: 3617
p2 <- ST.FeaturePlot(se.subset, features = c("nFeature_RNA",'nCount_RNA','percent.mito','percent.ribo'), dark.theme = T,
cols = c("black", "blue", "cyan", "yellow", "red", "darkred"))+ ggtitle(" subset ")
p3 <- ST.FeaturePlot(se, features = c("nFeature_RNA",'nCount_RNA','percent.mito','percent.ribo'), dark.theme = T,
cols = c("black", "blue", "cyan", "yellow", "red", "darkred")) + ggtitle("no subset ")
cowplot::plot_grid(p2,p3)
这上面可能看的不清楚，我用PPT组合了另一个项目的（也是10X的官方数据），subset前后的图像，可以感受到数量与质量的矛盾：考虑到一个切片3000多个spot就要好几万，这样贸然去掉是不是不太好？所以在我们的分析中需要慎重考虑，另外需要注意的是我们的spot 并不是一个细胞。
其实不好的亚群，说到底还是该亚群的基因不make sense。如果你有很好的理由去除掉某种类型的基因，这也很容易做到。例如，您可能希望在数据集中只保留蛋白质编码基因。这里我们作为演示，做一下，我们的宗旨还是不先入为主地去除基因和细胞。
ensids <- read.table(file = list.files(system.file("extdata", package = "STutility"), full.names = T, pattern = "mouse_genes"), header = T, sep = "\t", stringsAsFactors = F)
# 注意这里是小鼠的不是人，这里我们演示了如何使用预定义的covernersion表将Seurat对象划分为子集，只包含蛋白质编码基因，但是您也可以从其他地方获得这一信息，例如bioMart。
head(ensids)
gene_id gene_type gene_name
1 ENSMUSG00000102693.1 TEC 4933401J01Rik
2 ENSMUSG00000064842.1 snRNA Gm26206
3 ENSMUSG00000051951.5 protein_coding Xkr4
4 ENSMUSG00000102851.1 processed_pseudogene Gm18956
5 ENSMUSG00000103377.1 TEC Gm37180
6 ENSMUSG00000104017.1 TEC Gm37363
# Print available biotypes
unique(ensids$gene_type)
keep.genes <- subset(ensids, gene_type %in% "protein_coding")$gene_name
keep.genes
# Subset Seurat object
se.subset <- se[intersect(rownames(se), keep.genes), ]
cat("Number of genes removed : ", nrow(se) - nrow(se.subset), "\n")
图像处理
下面我认为是STUtility 的核心了：方便地操纵我们的图像。在看到STUtility之前，我们可能不知道原来图像还可以这样处理。创建Seurat对象之后，我们就可以从infoTable中提供的文件路径中加载H&E映像。LoadImages()函数允许您将图像加载到Seurat对象中，并自动保存每个图像的压缩版本，您可以用于绘图。
载入图像
se <- LoadImages(se, time.resolve = F, verbose = T)
ImagePlot(se, method = "raster", darken = TRUE, type = "raw")
我们可以画一个marker看看：
FeatureOverlay(se,
sampleids = 1:2,
features = "CD27",
pt.size = 1.5,
cols = c("dark blue", "cyan", "yellow", "red", "dark red"),
dark.theme = T,
type = "raw",ncols.samples = 2)
出于可视化目的，另一个有用的特性是屏蔽HE图像的背景。函数MaskImages()可以移除背景，目前对于边界清晰的组织来说效果很好。掩蔽是一个重要的问题，由于染色，切片等原因，在某些组织类型中它可能会失败。如果发生这种情况，你可以尝试修改MaskImages()中的参数，或者创建你自己的屏蔽函数，看看你是否可以获得更好的结果。
所提供的掩蔽（masking）算法使用SLIC方法，属于简单的线性迭代聚类(更多信息，参见（https://www.r-bloggers.com/superpixels-in-imager/）)。该算法根据图像平面上的颜色相似度和接近度对像素进行聚类，生成超像素。在对HE图像运行SLIC方法之前，我们注意到一些预处理可以显著改善掩蔽。
除了刚才看到的指未经任何修改的原始图像（ method = "raster"），还有以下四种：
se <- MaskImages(object = se)
ImagePlot(se, ncols = 2, method = "raster", type = "masked", darken = T) # Masked image
ImagePlot(se, ncols = 2, method = "raster", type = "masked.masks") # Mask
自动对齐方法(AlignImages())首先尝试从每个图像中检测组织的边缘。默认情况下将第一个图像(reference.index = 1)作为参考，但是你可以用reference.index参数指定哪个图像作为参考。然后对于每幅图像，学习一个变换矩阵，来映射坐标到参考图像。这种对齐方法有时会失败，在这种情况下，你可以使用ManualAlignImages()函数手动对齐图像。在以下情况下通常需要这样做:
(1)组织比图像背景大，因此部分组织在帧边缘之外
(2)组织具有对称的形状(例如，如果组织是圆形，使用组织边缘将很难找到最佳对齐)
(3)masking 失败
实际的转换是使用imager R包中的imwarp()函数使用“反向”转换策略完成的。这种方法确保每个像素都是使用线性插值绘制，注意，这样对齐后的图像将经历一些质量损失。
se <- se %>% MaskImages() %>% AlignImages()
ImagePlot(se, method = "raster", darken = TRUE, type = "processed", ncols = 2)
ImagePlot(se, method = "raster", darken = TRUE, type = "processed.masks", ncols = 2)
图像的旋转。因为切片的时候或者贴玻片的时候可能并没对，这时候就可以用这个来做了。
transforms <- list("2" = list("angle" = 90))
se.rotate90 <- WarpImages(se, transforms)
ImagePlot(se.rotate90, method = "raster", darken = T)
transforms <- list("2" = list("mirror.x" = T))
se.mirrorx <- WarpImages(se, transforms)
ImagePlot(se.mirrorx, method = "raster", darken = T)
transforms <- list("2" = list("mirror.y" = T))
se.mirrory <- WarpImages(se, transforms)
ImagePlot(se.mirrory, method = "raster", darken = T)
旋转后的对象可以用来绘制基因表达，分群等。
FeatureOverlay(se.mirrorx, features = "CD27",
sampleids = 1:2,
cols = c("dark blue", "cyan", "yellow", "red", "dark red"),
method = "raster",
dark.theme = T, ncols.samples = 2)
手动绘制感兴趣的区域
在STUtility中包含了一个用于手动注释的shiny应用程序。它允许用户选择并给出一个/几个特定的捕捉点标签，这可以用于可视化或DEA目的。默认情况下，应用程序将以浏览器模式打开。注释完成后，只需关闭浏览器窗口并返回R。这里我们不再演示。
se <- ManualAnnotation(se)
3D 可视化
STutility允许通过使用Create3DStack()函数在对齐后的HE图像中检测到的（细胞）核来实现特征的3D可视化。坐标是根据颜色强度提取的，因为核的颜色通常比周围组织的颜色深。这不是一个精确的细胞分割，但它可以很好地捕捉到细胞核的密度，从而确定各种形态结构。
一旦核坐标从每个对齐的部分被提取出来，一个z值将被分配到每个部分来创建一个3D堆栈。特征值可以在堆栈的各个点上插值，然后通过将值映射到一个颜色尺度来实现3D可视化。下面是一些标准，必须满足：
切片必须来自相同的组织类型，每个切片的形态必须相似 （切的方向要一致） 图像必须是对齐的，也就是说你必须先运行LoadImages()， MaskImages()和AlignImages()(或者ManualAlignImages()) 图片必须以高于默认400像素的分辨率加载进来。如果图像低于400像素，Create3DStack()会自动重新加载更高分辨率的图像，或者在运行Create3DStack()之前可以运行SwitchResolution()来重新加载更高分辨率的图像。 细胞分割是基于颜色强度的，因此图像中有任何的噪音如毛发，褶皱，气泡或灰尘等将会带来影响，同时切片厚度和染色不均匀也会影响分割效果。 我们假定组织切片已用苏木精和伊红染色
se <- Create3DStack(se)
stack_3d <- setNames(GetStaffli(se)@scatter.data, c("x", "y", "z", "grid.cell"))
ggplot(stack_3d, aes(x, 2e3 - y)) +
geom_point(size = 0.1, color = "lightgray") +
facet_wrap(~z, ncol = 1) +
theme_void() +
theme(plot.background = element_rect(fill = "black"),
plot.title = element_text(colour = "white"),
legend.text = element_text(colour = "white"))
ggplot(interpolated.data, aes(x, 2e3 - y, color = val)) +
geom_point(size = 0.1) +
facet_wrap(~z, ncol = 2) +
theme_void() +
ggtitle("CD27") +
scale_color_gradientn(colours = c("black", "dark blue", "cyan", "yellow", "red", "dark red")) +
theme(plot.background = element_rect(fill = "black"),
plot.title = element_text(colour = "white"),
legend.text = element_text(colour = "white"))
下面是多张切片对齐后的空间表达，这是一张由plotly绘制的3D交互动图：
FeaturePlot3D(se, features = "CD27", dark.theme = TRUE, pt.size = 0.6, pts.downsample = 5e4)
References
[1] Bergenstråhle, J., Larsson, L. & Lundeberg, J. Seamless integration of image and molecular analysis for spatial transcriptomics workflows. BMC Genomics 21, 482 (2020).: https://doi.org/10.1186/s12864-020-06832-3
[2] https://rdcu.be/ccSyI
[3] https://ludvigla.github.io/STUtility_web_site/index.html
[4] https://www.spatialresearch.org/
[5] https://support.10xgenomics.com/spatial-gene-expression/datasets
看完记得顺手点个“在看”哦！