一、Seurat v3.2 对空间转录组Visium的结果分析中实现的功能：

• 归一化
• 降维和聚类
• 检测空间可变特征(spatially-variable features)
• 互动式的可视化
• 与单细胞RNA-seq数据整合
• 处理多个切片

二、安装Seurat v3.2 ：

#For learning Seurat v2.3
#Author:Robin 2019.12.22
# Enter commands in R (or R studio, if installed)

# Install the devtools package from Hadley Wickham
install.packages("devtools")

#devtools::install_github("satijalab/seurat", ref = "spatial")
devtools::install_github('satijalab/seurat-data')

成功安装！

> library(Seurat)
> packageVersion("Seurat")
[1] ‘3.2.2’

加载包：

library(Seurat)
library(SeuratData)
library(ggplot2)
library(cowplot)
library(dplyr)

三、下载数据集

• 从10X Visium流程中生成结果文件的方法：
  先来回顾一下跑完Space Ranger得到哪些结果文件：

├── analysis
│   ├── clustering
│   ├── diffexp
│   ├── pca
│   ├── tsne
│   └── umap
├── cloupe.cloupe
├── filtered_feature_bc_matrix
│   ├── barcodes.tsv.gz
│   ├── features.tsv.gz
│   └── matrix.mtx.gz
├── filtered_feature_bc_matrix.h5  #实际上需要的文件①
├── metrics_summary.csv
├── molecule_info.h5
├── possorted_genome_bam.bam
├── possorted_genome_bam.bam.bai
├── raw_feature_bc_matrix
│   ├── barcodes.tsv.gz
│   ├── features.tsv.gz
│   └── matrix.mtx.gz
├── raw_feature_bc_matrix.h5
├── spatial  # 实际上需要的文件②；记录空间信息:用户提供的原始全分辨率brightfield图像的下采样版本。下采样是通过box滤波实现的，它对全分辨率图像中像素块的RGB值进行平均，得到下采样图像中一个像素点的RGB值。
│   ├── aligned_fiducials.jpg 
│   ├── detected_tissue_image.jpg
│   ├── scalefactors_json.json
│   ├── tissue_hires_image.png
│   ├── tissue_lowres_image.png
│   └── tissue_positions_list.csv
└── web_summary.html

9 directories, 20 files

用seurat进行下游分析主要用到两个结果文件。一个是filtered_feature_bc_matrix.h5文件，一个是spatial镜像结果目录。

• 从10x官网下载测试数据

这里使用从10x官网下载的小鼠脑组织样本MouseBrain Serial Section 1 (Sagittal-Posterior)。

下载网址：https://support.10xgenomics.com/spatial-gene-expression/datasets

点击选择的样本，下载两个数据就行
cell matrix HDF5 (filtered)和Spatial imaging data

image.png

导入R包读取文件

library("Seurat")
library("ggplot2")
library("cowplot")
library("dplyr")
library("hdf5r")
##读取矩阵文件
name='Posterior1'
expr <- "/data/Sagittal-Posterior1/V1_Mouse_Brain_Sagittal_Posterior_filtered_feature_bc_matrix.h5"
expr.mydata <- Seurat::Read10X_h5(filename =  expr )
mydata <- Seurat::CreateSeuratObject(counts = expr.mydata, project = 'Posterior1', assay = 'Spatial')
mydata$slice <- 1
mydata$region <- 'Posterior1' #命名
#读取镜像文件
imgpath <- "/data/Sagittal-Posterior1/spatial"
img <- Seurat::Read10X_Image(image.dir = imgpath)
Seurat::DefaultAssay(object = img) <- 'Spatial'
img <- img[colnames(x = mydata)]
mydata[['image']] <- img
mydata  #查看数据
An object of class Seurat
32285 features across 3355 samples within 1 assay
Active assay: Spatial (32285 features)

从mydata的输出信息我们可以知道，这个样本包含3355个spot点、32285个基因。

也可以直接加载自己的数据：

brain <-  Load10X_Spatial(data.dir = a,
                  filename = 'V1_Mouse_Brain_Sagittal_Anterior_Section_2_filtered_feature_bc_matrix.h5',
                  assay = "Spatial",
                  slice = 'spatial',
                  filter.matrix = T)

四、基础统计作图

##UMI统计画图
plot1 <- VlnPlot(mydata, features = "nCount_Spatial", pt.size = 0.1) + NoLegend()
plot2 <- SpatialFeaturePlot(mydata, features = "nCount_Spatial") + theme(legend.position = "right")
plot_grid(plot1, plot2)

image.png

UMI数大部分集中到10000-20000区间，最高不超过80000，并且组织中高UMI数的区域主要集中在左下角。后面可以特意性关注一下左下角区域的基因的表达和主要的细胞类型。

##gene数目统计画图
plot1 <- VlnPlot(mydata, features = "nFeature_Spatial", pt.size = 0.1) + NoLegend()
plot2 <- SpatialFeaturePlot(mydata, features = "nFeature_Spatial") + theme(legend.position = "right")
plot_grid(plot1, plot2)

image.png

基因数目大部分处于2500-7500之间，结合UMI数据的分布可以发现UMI数目高的区域基因数也高，说明基因数和UMI数基本上是呈正相关的。

#线粒体统计
mydata[["percent.mt"]] <- PercentageFeatureSet(mydata, pattern = "^mt[-]")
plot1 <- VlnPlot(mydata, features = "percent.mt", pt.size = 0.1) + NoLegend()
plot2 <- SpatialFeaturePlot(mydata, features = "percent.mt") + theme(legend.position = "right")
plot_grid(plot1, plot2)

注意如果是人的数据pattern = "^mt[-]改成pattern = "^MT[-]

image.png

总体来说，这个样本的线粒体比例不高，左边中上区域有一处线粒体比例稍微高一点，后面也可以仔细研究一下这一块区域到底是特定的细胞类型引起的还是组织活性的差异引起的。不过从这张图我们还可以发现一个有意思的现象，基因和UMI高表达的区域往往线粒体比例更低。

五、数据过滤

做单细胞RNAseq我们都会根据UMI、基因数、线粒体比例等进行过滤，那么做空间转录组数据分析其实我们也可以按这样的方式来过滤。具体的过滤条件需要根据具体样本数据来定，没有固定的标准。
比如这个样本我们可以设置过滤条件：
① 基因数大于200，小于7500
② UMI数大于1000，小于60000
③ 线粒体比例小于25%

mydata2 <- subset(mydata, subset = nFeature_Spatial > 200 & nFeature_Spatial <7500 & nCount_Spatial > 1000 & nCount_Spatial < 60000 & percent.mt < 25)
mydata2
An object of class Seurat
32285 features across 2977 samples within 1 assay
Active assay: Spatial (32285 features)

过滤后还剩2977个spot点。过滤后我们在绘制一下UMI分布图。

plot1 <- VlnPlot(mydata2, features = "nCount_Spatial", pt.size = 0.1) + NoLegend()
plot2 <- SpatialFeaturePlot(mydata2, features = "nCount_Spatial") + theme(legend.position = "right")
plot_grid(plot1, plot2)

image.png

那么现在问题来了，过滤之后组织图像里面缺了几块，显得特别丑。那么我们到底应不应该过滤呢？过滤数据可以减少利群的点，减少对后面聚类结果的影响，不过滤数据可以让组织图像保持完整性，绘图更好看一点，所以这个还真不好决断。

六、数据归一化

Seurat官方推荐使用sctransform 归一化方法，它构建了基因表达的正则化负二项模型，以便在保留生物差异的同时考虑技术因素。
Sctransform函数同时实现了NormalizeData、ScaleData、FindVariableFeatures三个函数的功能。

有关sctransform的更多信息，请参见 here的预印和here的Seurat教程。sctransform将数据归一化，检测高方差特征high-variance features，并将数据存储在SCTassay中。

mydata <- SCTransform(mydata, assay = "Spatial", verbose = FALSE)

七、基因表达可视化

Seurat的SpatialFeaturePlot功能扩展了FeaturePlot，可以将表达数据覆盖在组织组织上。这里展示的Hpca基因是一个强的海马marker ，Ttr是一个脉络丛marker 。可以通过基因的表达分布来初步判断一下海马区和脉络丛区处于组织切片的哪个位置。

SpatialFeaturePlot(mydata, features = c("Hpca", "Ttr"))

image.png

从结果的展示来看，这两个marker基因的分布还是挺集中的，这也说明理由空间转录组数据来分析小鼠脑的不同区域的表达差异应该还是比较准确的。另外，海马区的分布可以大概分成3大块，从上之下第一块弧形区域似乎处于线粒体高表达区域，而最下面一块弧形区处于基因高表达区。后面可以把这三个不同区域的数据进行差异基因和功能的比较也许会发现一些有意思的东西。

7.1 Seurat中的默认参数强调分子数据的可视化。

但是，您还可以通过更改以下参数来调整斑点的大小（及其透明度），以改善组织学图像的可视化：

• pt.size.factor-这将缩放斑点的大小。默认值为1.6
• alpha-最小和最大透明度。默认值为c（1，1）。
• 尝试设置为alphac（0.1，1），以降低具有较低表达的点的透明度

p1 <- SpatialFeaturePlot(mydata, features = "Ttr", pt.size.factor = 1)
p2 <- SpatialFeaturePlot(mydata, features = "Ttr", alpha = c(0.1, 1))
plot_grid(p1, p2)

image.png

八、降维、聚类和可视化

先进行PCA降维，再选择前30个维度进行聚类和umap降维。

降维、聚类和可视化

mydata <- RunPCA(mydata, assay = "SCT", verbose = FALSE)
mydata <- FindNeighbors(mydata, reduction = "pca", dims = 1:30)
mydata <- FindClusters(mydata, verbose = FALSE)
mydata <- RunUMAP(mydata, reduction = "pca", dims = 1:30)

再然后，我们可以在UMAP空间(使用DimPlot)或使用SpatialDimPlot将聚类clustering的结果显示在图像上

p1 <- DimPlot(mydata, reduction = "umap", label = TRUE)
p2 <- SpatialDimPlot(mydata, label = TRUE, label.size = 3)
plot_grid(p1, p2)

image.png

由于颜色太多，因此可视化哪个立体像素voxel属于哪个cluster可能是一个挑战。
解决方案：

• 设置label参数会在每个cluster的中间处放置一个彩色框（请参见上图），
• 而SpatialDimPlot的do.hover参数允许交互式地查看当前点的标识。

# move your mouse
>SpatialDimPlot(mydata, do.hover = TRUE)
Error in SpatialDimPlot(mydata, do.hover = TRUE) : 
  参数没有用(do.hover = TRUE)

操作失败

image

由于亚群的颜色比较接近，有时候不太好判断，我们可以是cells.highlight来标记特定的亚群。

SpatialDimPlot(mydata, cells.highlight = CellsByIdentities(object = mydata, idents = c(1, 2, 3, 4, 
   5, 6)), facet.highlight = TRUE, ncol = 3)

image.png

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

LinkedDimPlot(brain)

image

tsne和umap两种展示方式在这次分析里差别不是特别大，tsne相对来说群于群之间分的更开，而umap则单个亚群位置更集中。这个时候我们也可以结合前面marker基因的表达分布图来大概判断一下每个亚群大概处于小鼠脑的那个区。

九、亚群marker基因分析

同时分析所有亚群的marker基因

data.markers <- FindAllMarkers(mydata, only.pos = FALSE, min.pct = 0.25, logfc.threshold = 0.25,test.use = "wilcox")

参数解析：
logfc.threshold: logfc阈值，阈值以上的基因才输出，默认0.25
only.pos： 是否只输出正的marker，也就是只输出在该亚群中高表达的marker，默认是FALSE
test.use： 算法选择，indAllMarkers提供8种差异分析算法，默认是wilcox
min.pct: 最低表达比例，指基因在两组细胞中至少有一组中有表达的细胞的比例等于或超过该阈值

查看每个亚群top3差异基因

topgene<-data.markers %>% group_by(cluster) %>% top_n(n = 3, wt = avg_logFC)
topgene
# A tibble: 48 x 7
# Groups:   cluster [16]
       p_val avg_logFC pct.1 pct.2 p_val_adj cluster gene 
       <dbl>     <dbl> <dbl> <dbl>     <dbl> <fct>   <chr>
 1 2.04e-245     2.03  1     0.911 3.73e-241 0       Plp1 
 2 3.04e-238     1.69  1     0.829 5.55e-234 0       Trf  
 3 1.13e-232     1.66  1     0.998 2.07e-228 0       Mbp  
 4 3.47e- 83     0.670 1     0.998 6.34e- 79 1       Mbp  
 5 3.81e- 77     0.643 1     0.818 6.96e- 73 1       Mobp 
 6 5.76e- 60     0.653 0.98  0.686 1.05e- 55 1       Agt  
 7 1.50e-161     2.14  1     0.655 2.74e-157 2       Pcp2 
 8 3.85e-155     1.81  1     0.808 7.04e-151 2       Pvalb
 9 6.11e-155     1.70  1     0.856 1.12e-150 2       Itpr1
10 3.04e-163     1.85  0.944 0.261 5.56e-159 3       Igf2 
# ... with 38 more rows

也可以单独两个或多个亚群比较分析marker基因分析

markers <- FindMarkers(mydata, ident.1 = 2,, ident.2 = 6)

绘制每个亚群top3差异基因热图

DoHeatmap(mydata, features = topgene$gene,size = 2) + NoLegend()

image.png

绘制marker基因小提琴图

VlnPlot(mydata,features = c('Plp1','Trf'), group.by = "seurat_clusters",pt.size = 0)

image.png

绘制marker基因tsne图

FeaturePlot(mydata, features = c('Plp1','Trf'), cols = c("grey", "red"),reduction = "tsne")

image.png

参考文章

• Seurat: Spatial Transcriptomics
• Seurat v3.1.2官网
• Seurat 新版教程：分析空间转录组数据
• Seurat4.0系列教程22：空间转录组的分析、可视化与整合
• Seurat(v4.0.4):空间转录组基本流程
• 重温seurat官方教程四：空间转录组分析
• 不可错过的单细胞转录组研究新维度：空间转录组
• 上课笔记|基于Space Ranger实现空间转录组数据分析

作者：邓老师呦
链接：https://www.jianshu.com/p/137376261b13
作者：Geekero
链接：https://www.jianshu.com/p/07593e4d99a9
来源：简书
著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。