植物空间转录组分析1：Seurat基本流程 - 简书 (jianshu.com)

植物空间转录组分析2：STEEL+Seurat - 简书 (jianshu.com)

空间转录组是2022生命科学十大创新产品名单，因此将来会在生物学领域有非常大的应用空间，目前植物类的相关文章较少，在本次的系列教程中，将使用复旦大学戚继团队兰花空间转录组的数据，希望大家一起学习，掌握植物空间转录组基本的分析方法(*^▽^*)
数据连接OSF | Spatiotemporal atlas of organogenesis in development of orchid flowers,与单细胞的数据结构基本一致，多了spatial这个文件夹，主要包含的就是切片信息和spot定位信息

文章
数据

1：前言

在本篇文章中，作者一共制作了兰花空间转录组切片，并使用了STEEL方法进行聚类，出于学习考虑，本次分析先使用Seurat的常见流程分析其中一个切片，后续推文中将使用文章中的STELL方法进行聚类并使用Seurat的其他函数进行后续分析
本次只分析Slide1，感兴趣的可以试试其他切片

Slide1

2：数据载入

和单细胞一样，只是多了spatial这个文件夹需要输入

##=============================1.安装依赖包=====================================
##BiocManager::install(c('Seurat','ggplot2','cowplot','dplyr'))
##BiocManager::install("monocle",force = TRUE) 
library(Seurat)
library(ggplot2)
library(cowplot)
library(dplyr)
library(magrittr)
library(gtools)
library(stringr)
library(Matrix)
library(tidyverse)
library(patchwork)
##=============================2.载入数据=======================================
orc1<- Load10X_Spatial("Slide_1",filename = "filtered_feature_bc_matrix.h5",assay = "spatial")

3: 质量控制

不只是能体现小提琴图，还能将每个spot的测序质量体现在切片上，帮我们确定哪些组织可能会出现问题

##=============================3.质量控制=======================================
dir.create("QC")
##UMI统计画图
plot1 <- VlnPlot(orc1, features = "nCount_Spatial", pt.size = 0.1) + NoLegend()
plot2 <- SpatialFeaturePlot(orc1, features = "nCount_Spatial") + theme(legend.position = "right")
pearplot <- plot_grid(plot1,plot2) 
ggsave("QC/Slide1_UMI.pdf", plot = pearplot, width = 12, height = 5) 
##Feature/Gene统计画图
plot1 <- VlnPlot(orc1, features = "nFeature_Spatial", pt.size = 0.1) + NoLegend()
plot2 <- SpatialFeaturePlot(orc1, features = "nFeature_Spatial") + theme(legend.position = "right")
pearplot <- plot_grid(plot1,plot2) 
ggsave("QC/Slide1_Feature.pdf", plot = pearplot, width = 12, height = 5)

UMI
feature

可以明显的看到，不管是UMI还是Feature，在spot中差异主要是有与切片组织的细胞数目以及细胞类型导致的。比如，在花瓣中可以看到非常少的UMI以及Feature，在中间的分生组织中，UMI和Feature的数目非常多。因此，单细胞的标准方法(如LogNormalize函数)可能会有问题，在空间转录组分析中，一般使用其他方法进行标准化。

4. 数据标准化

目前空间转录组数据推荐使用sctransform，具体的原理我还没有看懂，所以先根据流程走一遍，看看效果
为了探究标准化方法的不同,sctransform和log规范化结果如何与UMIs的数量相关。为了进行比较，首先运行sctransform来存储所有基因的值，并通过NormalizeData运行一个log规范化过程。

##============================4.数据标准化======================================
## 目前空间转录组推荐使用SCTransform，集合了单细胞标准化的NormalizeData()，ScaleData()，FindVariableFeatures()，
dir.create("Normalize")
orc1 <- SCTransform(orc1, assay = "Spatial", return.only.var.genes = FALSE, verbose = FALSE)
## 比较SCTransform与log Normalization的区别
orc1 <- NormalizeData(orc1, verbose = FALSE, assay = "Spatial")
## 计算UMI与Feature的相关性，分为5组
orc1 <- GroupCorrelation(orc1, group.assay = "Spatial", assay = "Spatial", slot = "data", do.plot = FALSE)
orc1 <- GroupCorrelation(orc1, group.assay = "Spatial", assay = "SCT", slot = "scale.data", do.plot = FALSE)
p1 <- GroupCorrelationPlot(orc1, assay = "Spatial", cor = "nCount_Spatial_cor") + ggtitle("Log Normalization") + 
  theme(plot.title = element_text(hjust = 0.5))
p2 <- GroupCorrelationPlot(orc1, assay = "SCT", cor = "nCount_Spatial_cor") + ggtitle("SCTransform Normalization") + 
  theme(plot.title = element_text(hjust = 0.5))
p3 <- plot_grid(p1, p2)
ggsave("Normalize/Normalization_cor_Slide1.pdf", plot = p3, width = 12, height = 5)  

标准化

对于上面的箱形图，计算每个特征(基因)与UMIs数量(这里的nCount_Spatial变量)的相关性。然后，根据基因的平均表达将它们分组，并生成这些相关性的箱形图。
同时这里需要注意，SCT标准化后有个单独的assay(SCT)，里面包含三个矩阵，分别为count,data以及scale.data，log标准化也有一个assay(Spatial)，里面也包含三个矩阵，但是scale.data的数据为0
标准化之后进行高变基因的筛选，这部分没什么可说的，与单细胞是一致的

top10 <- head(VariableFeatures(orc1), 10) 
plot1 <- VariableFeaturePlot(orc1) 
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE, size=2.5) 
ggsave("Normalize/VariableFeature_Slide1.pdf", plot = plot2, width = 9, height = 6)

高变基因筛选

5.数据降维与聚类

由于文章中使用了STEEL方法，所以这次的分析只是走走流程，之后会按照文章中的STEEL方法进行测试

##============================5.数据降维与聚类==================================
dir.create("cluster")
## PCA降维并提取主成分
orc1 <- RunPCA(orc1, features = VariableFeatures(orc1_new)) 
plot1 <- DimPlot(orc1, reduction = "pca", group.by="orig.ident") 
plot2 <- ElbowPlot(orc1, ndims=40, reduction="pca") 
pearplot <- plot_grid(plot1,plot2)
ggsave("cluster/PCA_Slide1.pdf", plot = pearplot, width = 13, height = 6)

pdf("cluster/PC_heatmap_Slide1.pdf",height = 12,width = 12)
DimHeatmap(orc1, dims = 1:9, nfeatures=10, cells = 200, balanced = TRUE)
dev.off()

## 细胞聚类
## 此步利用 细胞-PC值 矩阵计算细胞之间的距离，
## 然后利用距离矩阵来聚类。其中有两个参数需要人工选择，
## 第一个是FindNeighbors()函数中的dims参数，需要指定哪些pc轴用于分析，选择依据是之前介绍的cluster/pca.png文件中的右图。
## 第二个是FindClusters()函数中的resolution参数，需要指定一个数值，用于决定clusters的相对数量，数值越大cluters越多。
orc1 <- FindNeighbors(orc1, dims = 1:20) 
orc1 <- FindClusters(orc1, resolution = 0.9)

metadata <- orc1@meta.data
table(orc1@meta.data$seurat_clusters)
cell_cluster <- data.frame(cell_ID=rownames(metadata), cluster_ID=metadata$seurat_clusters)
write.csv(cell_cluster,'cluster/cell_cluster_Slide1.csv',row.names = F)

PCA

使用resolution = 0.9只生成了12个聚类，原文中有40个，大家在设置的时候可以把resolution调到2甚至3左右看看效果

## 非线性降维
## tsne
orc1 <- RunTSNE(orc1, dims = 1:20)
embed_tsne <- Embeddings(orc1, 'tsne')
write.csv(embed_tsne,'cluster/embed_tsne_Slide1.csv')

p1 <- DimPlot(orc1, reduction = "tsne", label = TRUE)
p2 <- SpatialDimPlot(orc1, label.size = 3, pt.size.factor = 1.3)
pearplot <- plot_grid(p1,p2)
ggsave("cluster/tsne_Slide1.pdf", plot = p2, width = 6, height = 6)

## UMAP
orc1 <- RunUMAP(orc1,dims = 1:20)
embed_umap <- Embeddings(orc1, 'umap')
write.csv(embed_umap,'cluster/embed_umap_Slide1.csv')

p1 <- DimPlot(orc1, reduction = "umap", label = TRUE)
p2 <- SpatialDimPlot(orc1, label = TRUE, label.size = 3, pt.size.factor = 1.3)
pearplot <- plot_grid(p1,p2)
ggsave("cluster/umap_Slide1.pdf", plot = pearplot, width = 13, height = 6)

TSNE
UMAP

6. 部分基因展示

此处是我比较纠结的地方，因为一开始得到的数据和原文中有点不符，后来发现可能的原因是文章中使用的展示数据是counts

## featureplot
dir.create("featureplot")
p1 <- SpatialFeaturePlot(orc1, features = "PAXXG070630",slot="counts",pt.size.factor = 1.6,min.cutoff = 0.1)
p2 <- SpatialFeaturePlot(orc1, features = "PAXXG070630",slot="data",pt.size.factor = 1.6,min.cutoff = 0.1)
p3 <- SpatialFeaturePlot(orc1, features = "PAXXG070630",slot="scale.data",pt.size.factor = 1.6,min.cutoff = 0.1)
pearplot <- plot_grid(p1,p2,p3,align = "v",axis = "b",ncol = 3)
ggsave("featureplot/SCTSlide1.pdf", plot = pearplot, width = 18, height = 7)

三种展示方式
原文

总结

目前先把前期的工作的进行了一下，原文中进行了两个拟时分析，但因为聚类不一样所以先不进行分析，研究完STEEL之后再进行拟时分析。
总体来讲与单细胞变化不大，有优势也有劣势吧，对于细胞分型来说，空间转录组有先天优势，可以直接根据切片信息判断细胞类型，但是空间转录组最小的单位spot，不能代表单个细胞，分辨率可能还存在问题

转载请注明>>>周小钊的博客, 打赏推荐博客