Seurat新版教程:分析空间转录组数据(上)

男,

一个长大了才会遇到的帅哥,

稳健,潇洒,大方,靠谱。

一段生信缘,一棵技能树,

一枚大型测序工厂的螺丝钉,

一个随机森林中提灯觅食的津门旅客。

1思考题:

2

3+ 如何将空间数据与表达数据关联在一起?

4+ 有了空间转录组数据,如何与单细胞转录组数据联用?

5+ 做了多层切片如何展示真实的三维空间的转录本信息?

随着转录组技术的发展,空间转录组已经正式走向商业化时代,作为单细胞数据分析的工具箱的Seurat与时俱进,也相应地开发了空间转录组分析的一套函数,让我们跟随卑微小王看看Seurat官网教程吧。

本教程演示如何使用Seurat v3.2分析空间解析的RNA-seq数据。虽然分析流程类似于Seurat的单细胞RNA-seq分析流程,但我们引入了交互可视化工具,特别强调了空间和分子信息的集成。本教程将介绍以下任务,我们相信这些任务在许多空间分析中都很常见:

归一化

降维与聚类

检测spatially-variable特性

交互式可视化

与单细胞RNA-seq数据集成

处理多个片(multiple slices)

我们分析了使用来自10x Genomics 的Visium技术【Visium technology(https://www.10xgenomics.com/spatial-transcriptomics/)】生成的数据集。我们将在不久的将来扩展Seurat以处理其他数据类型,包括SLIDE-Seq(https://science.sciencemag.org/content/363/6434/1463)、STARmap(https://science.sciencemag.org/content/361/6400/eaat5691)和MERFISH(https://science.sciencemag.org/content/362/6416/eaau5324)。

安装1devtools::install_github("satijalab/seurat", ref = "spatial")

这种方法,只能好运。直接下载安装包,本地安装。在遇到问题的时候,假装这是你写的R包。参考《R包开发》这本书,尝试本地重构:1devtools::load_all('H:\\singlecell\\Seurat\\seurat-master\\seurat')

但是,一般用不到这么复杂,下载安装包,本地安装基本就可以了:1[https://codeload.github.com/satijalab/seurat/legacy.tar.gz/spatial](https://codeload.github.com/satijalab/seurat/legacy.tar.gz/spatial)  # 下载地址

2install.packages("H:/singlecell/Seurat/satijalab-seurat-v3.1.1-302-g1cb8a3d.tar.gz", repos = NULL, type = "source")

注意,别和之前的包安装冲突啦,这个3.2还处在开发中(2019-12-19)1library(Seurat)

2library(SeuratData)

3library(ggplot2)

4library(cowplot)

5library(dplyr)

数据集

和以往一样,seurat为了使其教程的可用,会提供测试的已发表的数据集。在这里,我们将使用最新发布的使用Visium v1化学生成的sagital小鼠大脑切片数据集。有两个连续的前段和两个(匹配的)连续的后段。

您可以在这(https://support.10xgenomics.com/spatial-gene-expression/datasets

)下载数据,并使用Load10X_Spatial函数将其加载到Seurat。这将读取spaceranger管道的输出,并返回Seurat对象,该对象包含spot级别的表达数据以及相关的组织切片图像。您还可以使用我们的SeuratData包方便地访问数据,如下所示。安装数据集之后,可以键入?stxBrain以了解更多信息。1InstallData("stxBrain")

每次我在国内直接用这种方法下载数据集都没有成功,如上,下载安装包,本地安装:1 install.packages("H:/singlecell/Seurat/stxBrain.SeuratData_0.1.1.tar.gz", repos = NULL, type = "source")

2brain <- LoadData("stxBrain", type = "anterior1")

当然,我们不推荐这种读取数据的方法,毕竟没有人把我们的数据也打包成一个R包的样子,我们拿到的是10 X Space Ranger的输出结果: 1├── analysis

2│?? ├── clustering

3│?? ├── diffexp

4│?? ├── pca

5│?? ├── tsne

6│?? └── umap

7├── cloupe.cloupe

8├── filtered_feature_bc_matrix

9│?? ├── barcodes.tsv.gz

10│?? ├── features.tsv.gz

11│?? └── matrix.mtx.gz

12├── filtered_feature_bc_matrix.h5

13├── metrics_summary.csv

14├── molecule_info.h5

15├── possorted_genome_bam.bam

16├── possorted_genome_bam.bam.bai

17├── raw_feature_bc_matrix

18│?? ├── barcodes.tsv.gz

19│?? ├── features.tsv.gz

20│?? └── matrix.mtx.gz

21├── raw_feature_bc_matrix.h5

22├── spatial  # 空间信息全在这 :这些文件是用户提供的原始全分辨率brightfield图像的下采样版本。

23下采样是通过box滤波实现的,它对全分辨率图像中像素块的RGB值进行平均,得到下采样图像中一个像素点的RGB值。

24│?? ├── aligned_fiducials.jpg   这个图像的尺寸是tissue_hires_image.png。由基准对齐算法发现的基准点用红色高亮显示。此文件对于验证基准对齐是否成功非常有用。

25│?? ├── detected_tissue_image.jpg

26│?? ├── scalefactors_json.json

27│?? ├── tissue_hires_image.png 图像的最大尺寸为2,000像素

28│?? ├── tissue_lowres_image.png 图像的最大尺寸为600像素。

29│?? └── tissue_positions_list.csv

30└── web_summary.html

其实只要 Space Ranger的输出结果 filtered_feature_bc_matrix.h5文件和一个spatial文件夹就可以了,一如如?stxBrain中给出的示例: 1setwd("H:\\singlecell\\spaceranger\\V1_Mouse_Brain_Sagittal_Anterio/")

2  # Load the expression data

3  expr.url <- 'H:\\singlecell\\spaceranger\\V1_Mouse_Brain_Sagittal_Anterio/V1_Mouse_Brain_Sagittal_Anterior_filtered_feature_bc_matrix.h5'

4  expr.data <- Seurat::Read10X_h5(filename =  expr.url )

5  anterior1 <- Seurat::CreateSeuratObject(counts = expr.data, project = 'anterior1', assay = 'Spatial')

6  anterior1$slice <- 1

7  anterior1$region <- 'anterior'

8  # Load the image data

9  img.url <- 'H:\\singlecell\\spaceranger\\V1_Mouse_Brain_Sagittal_Anterio/V1_Mouse_Brain_Sagittal_Anterior_spatial.tar.gz'

10  untar(tarfile =  img.url)

11  img <- Seurat::Read10X_Image(image.dir = 'H:\\singlecell\\spaceranger\\V1_Mouse_Brain_Sagittal_Anterio/spatial')

12  Seurat::DefaultAssay(object = img) <- 'Spatial'

13  img <- img[colnames(x = anterior1)]

14  anterior1[['image']] <- img

15

16  anterior1

17

18An object of class Seurat

1931053 features across 2696 samples within 1 assay

20Active assay: Spatial (31053 features)

21

22}

如果这两个文件在同一个文件下的话,也可以这样读入: 1list.files("H:\\singlecell\\spaceranger\\V1_Mouse_Brain_Sagittal_Anterio") # 注意命名要正确

2

3# [1] "filtered_feature_bc_matrix.h5"                                  "spatial"

4# [3] "V1_Mouse_Brain_Sagittal_Anterior_filtered_feature_bc_matrix.h5" "V1_Mouse_Brain_Sagittal_Anterior_spatial.tar.gz"

5

6brain<-Seurat::Load10X_Spatial("H:\\singlecell\\spaceranger\\V1_Mouse_Brain_Sagittal_Anterio")

7?Load10X_Spatial

8brain

9

10An object of class Seurat

1131053 features across 2696 samples within 1 assay

12Active assay: Spatial (31053 features)

空间数据如何存储在Seurat中?

来自10x的visium数据包括以下数据类型:

通过基因表达矩阵得到一个点(spot )

组织切片图像(采集数据时H&E染色)

用于显示的原始高分辨率图像与低分辨率图像之间的比例因子。

在Seurat对象中,spot by基因表达矩阵与典型的“RNA”分析类似,但包含spot水平,而不是单细胞水平的数据。图像本身存储在Seurat对象中的一个images 槽(slot )中。图像槽还存储必要的信息,以将斑点与其在组织图像上的物理位置相关联。

数据预处理

在spot中基因表达数据进行初始的预处理步骤与典型的scRNA-seq相似。我们首先需要对数据进行规范化,以考虑数据点之间测序深度的差异。我们注意到,对于空间数据集来说,分子数/点的差异可能是巨大的,特别是在整个组织的细胞密度存在差异的情况下。我们在这里看到大量的异质性,这需要有效的标准化。

熟悉的函数名~1plot1 <- VlnPlot(brain, features = "nCount_Spatial", pt.size = 0.1) + NoLegend()

2plot2 <- SpatialFeaturePlot(brain, features = "nCount_Spatial") + theme(legend.position = "right")

3plot_grid(plot1, plot2)

这些图表明,分子计数(molecular counts)在点间的差异不仅是技术性的,而且还依赖于组织解剖。例如,组织中神经元消耗的区域(如皮层白质),可再生地显示出较低的分子计数。因此,标准方法(如LogNormalize函数)可能会有问题,因为它会强制每个数据点在标准化之后具有相同的底层“大小”。

作为一种替代方法,我们推荐使用sctransform (Hafemeister和Satija,已出版),它构建了基因表达的正则化负二项模型,以便在保留生物差异的同时考虑技术因素。有关sctransform的更多信息,请参见 here(https://www.biorxiv.org/content/10.1101/576827v2)的预印和here(https://satijalab.org/seurat/sctransform_vignette.html)的Seurat教程。sctransform将数据归一化,检测高方差特征,并将数据存储在SCT分析中。1brain <- SCTransform(brain, assay = "Spatial", return.only.var.genes = FALSE, verbose = FALSE)

sct 与log-归一化相比,结果如何?

为了探究规范化方法的不同,我们研究了sctransform和log规范化结果如何与UMIs的数量相关。为了进行比较,我们首先重新运行sctransform来存储所有基因的值(这将会比较慢),并通过NormalizeData运行一个log规范化过程。1# rerun normalization to store sctransform residuals for all genes

2brain <- SCTransform(brain, assay = "Spatial", return.only.var.genes = FALSE, verbose = FALSE)

3# also run standard log normalization for comparison

4brain <- NormalizeData(brain, verbose = FALSE, assay = "Spatial")

1# Computes the correlation of the log normalized data and sctransform residuals with the number

2# of UMIs

3brain <- GroupCorrelation(brain, group.assay = "Spatial", assay = "Spatial", slot = "data", do.plot = FALSE)

4brain <- GroupCorrelation(brain, group.assay = "Spatial", assay = "SCT", slot = "scale.data", do.plot = FALSE)

5p1 <- GroupCorrelationPlot(brain, assay = "Spatial", cor = "nCount_Spatial_cor") + ggtitle("Log Normalization") +

6    theme(plot.title = element_text(hjust = 0.5))

7p2 <- GroupCorrelationPlot(brain, assay = "SCT", cor = "nCount_Spatial_cor") + ggtitle("SCTransform Normalization") +

8    theme(plot.title = element_text(hjust = 0.5))

9plot_grid(p1, p2)

对于上面的箱形图,我们计算每个特征(基因)与UMIs数量(这里的nCount_Spatial变量)的相关性。然后,我们根据基因的平均表达将它们分组,并生成这些相关性的箱形图。您可以看到,log-normalization未能充分规范化前三组的基因,这表明技术因素影响高表达基因的规范化表达估计值。相反,sctransform规范化降低了这种效果。差别真的很大么?读者诸君自行判断。

基因表达可视化

在Seurat v3.2中,我们加入了新的功能来探索和与空间数据固有的可视化特性。Seurat的SpatialFeaturePlot功能扩展了FeaturePlot,可以将表达数据覆盖在组织组织上。例如,在这组小鼠大脑数据中,Hpca基因是一个强的海马marker ,Ttr是一个脉络丛marker 。1SpatialFeaturePlot(brain, features = c("Hpca", "Ttr"))

Seurat的默认参数强调分子数据的可视化。然而,你也可以调整斑点的大小(和它们的透明度)来改善组织学图像的可视化,通过改变以下参数:

pt.size。因素-这将比例大小的斑点。默认为1.6

alpha -最小和最大透明度。默认是c(1,1)

尝试设置为alpha c(0.1, 1),以降低表达式较低的点的透明度

降维、聚类和可视化

然后,我们可以使用与scRNA-seq分析相同的工作流,对RNA表达数据进行降维和聚类。我们可以在UMAP空间(使用DimPlot)或使用SpatialDimPlot将分群的结果显示在图像上。1brain <- RunPCA(brain, assay = "SCT", verbose = FALSE)

2brain <- FindNeighbors(brain, reduction = "pca", dims = 1:30)

3brain <- FindClusters(brain, verbose = FALSE)

4brain <- RunUMAP(brain, reduction = "pca", dims = 1:30)

5

1Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric

2To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'

3This message will be shown once per session

422:08:28 UMAP embedding parameters a = 0.9922 b = 1.112

522:08:28 Read 2696 rows and found 30 numeric columns

622:08:28 Using Annoy for neighbor search, n_neighbors = 30

722:08:28 Building Annoy index with metric = cosine, n_trees = 50

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

9[----|----|----|----|----|----|----|----|----|----|

10**************************************************|

1122:08:30 Writing NN index file to temp file C:\Users\ADMINI~1\AppData\Local\Temp\Rtmp08OC2k\file3ddc1e802bb6

1222:08:30 Searching Annoy index using 1 thread, search_k = 3000

1322:08:31 Annoy recall = 100%

1422:08:32 Commencing smooth kNN distance calibration using 1 thread

1522:08:33 Initializing from normalized Laplacian + noise

1622:08:33 Commencing optimization for 500 epochs, with 106630 positive edges

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

18[----|----|----|----|----|----|----|----|----|----|

19**************************************************|

2022:08:50 Optimization finished

1p1 <- DimPlot(brain, reduction = "umap", label = TRUE)

2p2 <- SpatialDimPlot(brain, label = TRUE, label.size = 3)

3plot_grid(p1, p2)

4

image

因为有许多颜色,所以可视化哪个颜色属于哪个簇是很有挑战性的。我们有一些策略来帮助解决这个问题。通过设置label参数,可以在每个集群的中间位置放置一个彩色框(参见上面的图)以及do.hover。SpatialDimPlot的悬停参数允许交互式查看当前的spot标识。 1# move your mouse

2SpatialDimPlot(brain, do.hover = TRUE)

3

4Warning messages:

51: In if (is.na(col)) { :

6  the condition has length > 1 and only the first element will be used

72: In if (is.na(col)) { :

8  the condition has length > 1 and only the first element will be used

93: `error_y.color` does not currently support multiple values.

104: `error_x.color` does not currently support multiple values.

115: `line.color` does not currently support multiple values.

126: The titlefont attribute is deprecated. Use title = list(font = ...) instead.

你也可以使用cells.highlight,用于在空间坐标图上划分感兴趣的特定单元格。这对于区分单个集群的空间定位非常有用,如下所示:

LinkedDimPlot和LinkedFeaturePlot函数支持交互式可视化。这些图将UMAP表示与组织图像表示联系起来,并允许交互选择。例如,您可以在UMAP图中选择一个区域,图像表示中相应的点将突出显示。

教程官网(https://links.jianshu.com/go?to=https%3A%2F%2Fsatijalab.org%2Fseurat%2Fv3.1%2Fspatial_vignette.html)

(0)

相关推荐