前面我们有介绍了利用10x Space Ranger软件分析空间转录组原始数据得到可用于下游分析的矩阵和镜像文件。今天来介绍一下怎么利用Space Ranger的结果文件进行后续分析，这里主要使用Seurat在进行下游分析。

先来回顾一下跑完Space Ranger得到哪些结果文件：

Outputs:
- Run summary HTML:                         /opt/sample345/outs/web_summary.html
- Outputs of spatial pipeline:              /opt/sample345/outs/spatial
- Run summary CSV:                          /opt/sample345/outs/metrics_summary.csv
- BAM:                                      /opt/sample345/outs/possorted_genome_bam.bam
- BAM index:                                /opt/sample345/outs/possorted_genome_bam.bam.bai
- Filtered feature-barcode matrices MEX:    /opt/sample345/outs/filtered_feature_bc_matrix
- Filtered feature-barcode matrices HDF5:   /opt/sample345/outs/filtered_feature_bc_matrix.h5
- Unfiltered feature-barcode matrices MEX:  /opt/sample345/outs/raw_feature_bc_matrix
- Unfiltered feature-barcode matrices HDF5: /opt/sample345/outs/raw_feature_bc_matrix.h5
- Secondary analysis output CSV:            /opt/sample345/outs/analysis
- Per-molecule read information:            /opt/sample345/outs/molecule_info.h5
- Loupe Browser file:                       /opt/sample345/outs/cloupe.cloupe

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

安装R包

由于seurat分析空间转录组的R包satijalab-seurat是在GitHub 上的，如果我们需要直接安装的话需要先安装R包devtools，然后利用devtools工具中的install_github来安装GitHub上的R包。

安装devtools

install.packages('devtools')

安装satijalab-seurat

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

考虑到直接安装github上的R包速度是很慢的，非常考验网速，可能需要多次才能安装成功，我们也可以直接下载安装包，本地安装。

#下载
https://codeload.github.com/satijalab/seurat/legacy.tar.gz/spatial 
#安装
install.packages("satijalab-seurat-v3.1.5-351-g85610bc.tar.gz", repos = NULL, type = "source")

注意该R包还在开发中，不要和之前安装的seurat包冲突。

数据准备

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

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

image.png

点击选择的样本，下载两个数据就行
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个基因。

基础统计作图

##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三个函数的功能。

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

基因表达可视化

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

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

降维、聚类和可视化

先进行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)
mydata <- RunTSNE(mydata, reduction = "pca",dims = 1:30)

tsne展示结果：

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

由于亚群的颜色比较接近，有时候不太好判断，我们可以是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函数进行交互式可视化，这里不再展示。

image.png

扫码关注公众号，内容更精彩哦！