晨曦的分析笔记：空间转录组实战（二）

Hi，大家好，我是晨曦
上期我们讲解了空间转录组的常规流程，大家可以很直观的发现，其中的大部分内容其实和我们scRNA-seq分析的常规流程是类似的，也就是说我们在有了scRNA-seq的基础上，可以更好的入门空间转录组分析。
疯狂暗示：快去看以前的推文~
拿去耍！！空间转录组实战来了！你在实验室的装逼利器！

那么本期内容我们继续聚焦在空间转录组，这次我们通过一篇文献来尝试进行空间转录组部分图片的可视化

注意：这篇推文并不是一个标准答案，而是晨曦针对空间转录组文献的一种exercise，其目的在于希望和各位小伙伴一起走一遍空间转录组分析流程，一起探索空间转录组的种种细节

那么，我们开始吧

ST 数据获取

我们选择的文献是一篇针对心脏的空间转录组学

第一步我们依旧是获取数据，但是在获取数据之前，我们首先需要对文章有一个大概的了解，所以我们可以从摘要入手快速获取文献的大致内容

晨曦解读

请各位小伙伴重点注意标黄色的内容，本篇文献选择了三个阶段的心脏做空间转录组学，所以从这里我们可以知道，后续的样本至少会有三个（Ps:但是绝对不会只有三个，因为对于这种类型的期刊，每一个种类的样本数量肯定会大于3），然后映射细胞类型到特定的转录组解剖区域，针对这里的理解则是把三个阶段的样本整合/合并以后进行细胞亚群的划分然后把这个大整体的细胞亚群映射到每一个样本上，最后则是作者生成了一个公开的网站可以供我们对数据进行进一步的探讨

（Ps:后期我们获得数据分析可能考虑的来源之一）

然后既然我们的目的其实就是为了获取数据来为我们的课题服务，我们就在文章中去寻找数据链接可能出现的地方，通常在文献的method或者Data区域就可以获得我们感兴趣的数据所在链接或者地址，这篇文献的数据内容如下：

从上面的内容可以知道，我们想要的表达矩阵储存在这个网址：

https://www.spatialresearch.org/resources-published-datasets/

然后我们进入这个网址找到我们想要的数据集，如下：

我们依旧是选择Filtered后的文件，这里为什么要选择这个文件呢？

ST数据是通过组织切片透化以后在每一个spot上标注位置信息来达到存储位置信息的目的，那么组织切片并不是满满当当的在一个切片内，所以筛选的含义就是扣除了那些没有位置信息的spot，这样就可以直观的展示原始切片的形状

然后我们把这个文件下载下来以后，剩下的步骤就是R语言中进行操作了，但是在操作之前，我们来看一下，我们下载的文件究竟是一个什么样子（一个文件中包含了两个文件，分别是matrix以及meta）

matrix文件

观测为Gene ID——Ensembl ID信息

变量为位置信息——因为每一个位置上按照道理只会出现一个spot或者细胞，所以只要指代唯一，那么就是可以的，所以这里变量信息给我们提供的就是Barcode ID，只不过是用位置信息代替的（17✖20可能代表X轴为17，Y轴为20的这个点的某Gene ID表达量为0）

上一期推文，晨曦已经说过，我们不用执着于文件的形式，而是要关注文件所提供的信息，因为万变不离其宗

scRNA-seq提供给我们的其实就是一个表达矩阵的信息，变量为细胞ID（barcode ID），观测为基因信息（Gene symbol），表达数值为稀疏矩阵的形式

ST提供给我们的其实就是表达矩阵的信息外加位置信息，观测为基因信息（Gene symbol/ID），变量为细胞ID（barcode ID或者位置信息），这里为什么可以用位置信息来代替细胞ID呢，

举个例子，一个教室内小明特指一个同学，同一个教室第三排第二座也特指一个同学，所以这两个信息其实是可以替换的，但是你用小明在R中就无法转换到第三排第二座，因为小明可以理解为是一个一维信息，而第三排第二座是一个二维信息，高维可以降维打击，低维无法向上跃迁（仅供理解，没有实际根据）

meta文件

上面的信息晨曦看了一下，应该是样本的一些额外信息，比如说线粒体比例、UMI数量以及样本的额外信息等等，这些信息再我们后续分析的时候，如果有需要可以再重点探索一下，目前来看我们后续的分析可能并不会需要上面这个文件

注意：不知道各位小伙伴有没有发现，我们用wps读取meta文件的时候居然“串列”了！！简单的用WPS操作了一下，下面这个才是meta应该有的样子

晨曦画箭头的地方就是我们的空间位置信息，这时候可能有小伙伴有疑问，这里的信息为什么和观测不一样？

可以很清楚的看到观测信息是位置信息省略小数后的产物，晨曦在这里说的目的是因为有的作者观测时barcode ID，那么位置信息我们就需要从另一个文件中获取，请牢记不需要纠结数据的形式，需要注意数据的内容，毕竟万变不离其宗

提问：晨曦，为什么二维空间的位置信息是一个三维表示的形式，这个是什么意思？

回答：其实我们后期可视化或者阅读文献可以知道，这篇文献的空转数据有多个样本，所以观测所表现出来的形式应该解读为：样本ID✖X坐标✖Y坐标，然后我们通过对文献的解析我们发现，这篇文献作者采用了心脏不同时期、不同数量的心脏切片作为数据来源

既然我们已经获取了文件，那么我们接下来就转战到R语言进行后续相关分析的操作

ST 数据分析

#准备工作library(Seurat)#scRNA-seq以及ST分析的主体library(ggplot2)#绘图library(patchwork)#拼图library(data.table)#读取数据library(tidyverse)#整理数据library(dplyr)#整理数据

晨曦解读

工欲善其事必先利其器，一套组合拳走起~

#读取数据（这一块我们主要探索一下数据，就是玩~）#探索meta数据（空间位置信息）position_ST_ninth <-fread("meta_data.tsv")#读取空间信息文件head(position_ST_ninth[1:10,1:15])position_ST_ninth <-as.data.frame(position_ST_ninth)#数据框rownames(position_ST_ninth) <- position_ST_ninth[,1]position_ST_ninth <- position_ST_ninth[,-1]#因为是多样本所以不能放在一张图上展示，我们选择第九周的其中一个样本作为研究position_ST_ninth <- position_ST_ninth[position_ST_ninth$weeks==9,]position_ST_ninth <- position_ST_ninth[position_ST_ninth$Sample == "FH9_1000L3_CN31_E2",]#可视化p1=ggplot(position_ST_ninth, aes(x = position_ST_ninth$new_x, y = position_ST_ninth$new_y)) + geom_point(data = position_ST_ninth,color="red")+ labs(x = "x",y = "y")+ theme_bw()p1

晨曦解读

左边的图是文章中的可视化结果，右边的图则是我们单纯通过空间信息得到的可视化结果，可以很清楚的看到基本上除了配色可以说是一摸一样，所以也证明了我们通过空间位置信息绘制可视化的准确性，以及这个文件究竟给我们带来了什么信息——空间位置坐标，但是所谓的HE图片信息则是看作者有没有提供，但是很遗憾，本篇文献的作者并没有提供HE的一些信息

#那么接下来，我们已经明确了空间位置信息，那么接下来我们就可以来处理表达矩阵，进而完成空间转录组的标准流程

#准备工作library(AnnotationDbi)library(org.Hs.eg.db)library(tidyverse)library(Seurat)library(ggplot2)library(patchwork)library(data.table)library(tidyverse)library(dplyr)library(ggsci)#读取数据（位置信息）position_ST <-fread("meta_data.tsv")#读取空间信息文件head(position_ST[1:10,1:15])position_ST <-as.data.frame(position_ST)#数据框rownames(position_ST) <- position_ST[,1]position_ST <- position_ST[,-1]position_ST_ninth <- position_ST[position_ST$weeks==9,]position_ST_ninth <- position_ST_ninth[position_ST_ninth$Sample == "FH9_1000L3_CN31_E2",]#读取数据（表达矩阵）heart_ex <- fread("filtered_matrix.tsv")head(heart_ex[1:10,1:15])heart_ex <-as.data.frame(heart_ex)#数据框rownames(heart_ex) <- heart_ex[,1]heart_ex <- heart_ex[,-1]#数据预处理(去重+探针注释)dat <- data.frame(str_split(rownames(heart_ex), '[.]', simplify = T))rownames(heart_ex) <- dat$X1dat$symbol <- mapIds(org.Hs.eg.db, keys=dat$X1, column="SYMBOL", keytype="ENSEMBL", multiVals="first")dat <- dat[!is.na(dat$symbol),]dat <- dat[!duplicated(dat$symbol),]eq <-intersect(dat$X1,rownames(heart_ex))heart_ex=heart_ex[eq,]identical(dat$X1,rownames(heart_ex))rownames(heart_ex) <- dat$symbol#创建Visium对象并进行标准流程heart_ST <- CreateSeuratObject(heart_ex)heart_ST <- SCTransform(heart_ST, assay = "RNA", verbose = FALSE)heart_ST <- RunPCA(heart_ST, assay = "SCT", verbose = FALSE) plot1 <- DimPlot(heart_ST , reduction = "pca", group.by="orig.ident")plot2 <- ElbowPlot(heart_ST , ndims=20, reduction="pca") plot1+plot2pc.num=1:30heart_ST <- FindNeighbors(heart_ST , reduction = "pca", dims = pc.num)heart_ST <- FindClusters(heart_ST , verbose = FALSE,resolution = 0.6)#resolution参数越小亚群越少heart_ST <- RunTSNE(heart_ST, reduction = "pca", dims = pc.num)p1 <- DimPlot(heart_ST, reduction = "tsne", label = TRUE)p1

#因为我们这个数据的格式并不是标准的H5对象，所以我们只能通过另外的位置信息文件来获得位置信息

#进而标注我们的亚群在空间位置上heart_ST_meta <[email protected]position_ST_ninth <-position_ST_ninth[rownames(heart_ST_meta),]heart_ST_meta$x=position_ST_ninth$new_xheart_ST_meta$y=position_ST_ninth$new_yp2 <-ggplot(heart_ST_meta, aes(x = x, y = y)) + geom_point(data = heart_ST_meta, aes(color=seurat_clusters))+ theme_bw()p2 + scale_color_aaas()

晨曦解读

配色以及亚群分类的情况稍微有所不同，但是生信毕竟复现不可能是100%的，但是做到这里，我本人却产生了两个问题：

问题1：各位小伙伴仔细看代码就可以知道，晨曦这里是把表达矩阵直接导入生成的ST对象，但是实际上表达矩阵是多个样本多个时期合并后的产物，直接这么导入是否存在问题，这一点我们并没有考虑到？

问题2：空间位置信息如果没有导入进ST对象，后期寻找空间高可变基因的流程会直接出现报错，这一点又该如何解决呢？

那么，我们继续探索来解决上面这两个问题！！

解决问题1的思考流程如下：

我们试图从文献中寻找相关答案：

晨曦解读

文献中提到作者分别收集三个时期的心脏组织，每个时期获取的组织切片分别是四、九和六，所以总共应该是存在十九个样本，然后其实我们纠结的点在于，我们这些样本再不去除批次效应的前提下是否具有可比性，因为单细胞在进行多数据集联合分析的时候，其实我们并不是每次都需要去除批次效应的，一般来说如果说相同样本之间存在明显的差异，那么我们需要去除批次效应来平衡非生物学差异，但是如果样本与样本之间高度相关即可以认为其差异可能并非其他因素造成的，那么这个时候，我们可能只需要使用SCT联合merge函数来进行数据集的合并，那么情况究竟如何，我们继续从文献中探索：

晨曦解读

文献中指出三个样本之间的表达非常相似，聚类实际上来说是根据机器学习无监督的方式，其本质我们可以参考热图的观测聚类，是通过表达模式，那么这里文献中所说的表达非常相似，就可以等同于样本与样本之间的聚类并不存在太大的非生物差异，那么从这里来看，我们并不需要使用去除批次效应的算法（例如CCA等），而是可以直接使用SCT的方法即可，但是从一句话来看未免太过武断，我们继续往下看

晨曦解读

文献的最后列出了数据处理的流程（PS：不愧是Cell，真的贴心），对妊娠后三个时间点的ST数据联合（看成一个整体）进行归一化、降维和聚类，说明了把数合并起来以后进行归一化、降维、聚类，而我们使用SCT是替换了NormalizeData(), ScaleData(), and FindVariableFeatures()这三个过程，所以整体的流程是合理的，而文献中对于数据中来源的技术性以及生物学差异的平衡时直接使用ScaleData()来进行出处理，这一个步骤我们也已经进行过，所以通过对文献的解析，我们可以证明我们的分析流程是合理的

解决问题2的思考流程如下：

解决了问题1继续来看第二个问题，我们这里运行下面的代码看看问题的error：

GetTissueCoordinates(heart_ST)#表明我们的对象并不是10×Visium打包好的h5对象#Error: No images present in this Seurat object

说明我们的对象并不是10×标准的Visium对象，所以我们需要往Visium对象内添加内容，好完成后续的空间高可变基因的筛选

#h5对象内无空间位置信息GetTissueCoordinates(heart_ST)#表明我们的对象并不是10×Visium打包好的h5对象#Error: No images present in this Seurat object#我们来探索一下标准的10×Visium信息Mouse_Brain <- Load10X_Spatial(data.dir ="spatial_RNAseq/", filename = "Visium_FFPE_Mouse_Brain_filtered_feature_bc_matrix.h5", slice ="Mouse_Brain")head(GetTissueCoordinates(Mouse_Brain),3)# imagerow imagecol#AAACAGAGCGACTCCT-1 412.7678 391.5173#AAACCCGAACGAAATC-1 479.1430 220.3323#AAACCGGGTAGGTACC-1 203.2138 237.2540#可以很清楚的看到，我们标准的10×Visium数据中的位置信息是以这样的方式储存的#也就是说只要我们把位置信息整理成这样的形式就可以补充到我们缺失的h5对象中去head(Mouse_Brain@images$Mouse_Brain@coordinates)#会发现这个数据就很全，包括了组织标志以及坐标信息等等#既然知道了组成信息，那么我们接下来就尝试组建位置信息position_ST <- position_ST[,c(3,11,12)]colnames(position_ST) <- c("tissue","imagerow","imagecol")head(position_ST,3)position_ST$tissue <- str_split_fixed(a$positon, "x", 2)[,1]#数据整理完毕，我们开始进行数据的导入heart_ST@images <- list(position_ST)GetTissueCoordinates(heart_ST)#依旧报错（说明我们的数据还是没有整理妥当）#继续排查发现正统的10×Visium文件提取图片信息会有一个描述信息，考虑端口设置错误Mouse_Brain@images#$Mouse_Brain#Spatial data from the VisiumV1 technology for 2264 samples#Associated assay: Spatial #Image key: mousebrain_box <-Mouse_Brain@images$Mouse_Brainbox@coordinates <- position_STheart_ST@images <- list(heart = box)#这样我们的接口就完全匹配了，而且我们看一下样本3111也对应齐全

#唯一不完美的就是imagekey，但是无关紧要只不过是一个标志而已，如何修改其实也很简单，只需要把前面的导入的数据变量名称替换成heart即可

#为了教程直观没有替换，不影响GetTissueCoordinates(heart_ST)#搜寻空间高可变基因heart_ST <- FindSpatiallyVariableFeatures(heart_ST, assay = "SCT", features =VariableFeatures(heart_ST)[1:100], image = NULL, selection.method = "markvariogram", x.cuts = NULL, y.cuts = NULL)top.features <- head(SpatiallyVariableFeatures(heart_ST, selection.method = "markvariogram"), 2)SpatialFeaturePlot(heart_ST, features = top.features, ncol = 2, alpha = c(0.1, 1))