晨曦的空间转录组笔记(三):scRNA-seq联合Visium
Hi,大家好,我是晨曦
今天这期的空间转录组笔记,让我们来一起探索一下如何把我们的单细胞数据和空间转录组数据进行一个联合分析,好让我们的文章层次可以更上一个档次。
晨曦的空间转录组笔记系列传送门
Ps:本篇推文为晨曦学习空间转录组学分析的笔记,里面的观点如果和各位小伙伴的观点不同,欢迎各位小伙伴在评论区进行讨论~
那么,我们开始吧
引言
首先,让我们用简练的话语来回答一个问题:为什么我们要把空间转录组和单细胞转录组联合起来?
这个问题的答案就是,空间转录组优于单细胞转录组的一方面在于提供了空间位置信息,但是目前来看,空间转录组并没有做到完全的单细胞,其每一个spot大概会有1-10个细胞左右,所以并不是完全的“单细胞”,这个时候,我们通过算法,可以把单细胞的精度和空间转录组的位置信息联合起来,这样我们就可以在一定程度上达到——单细胞+空间转录组
那么,解决了上一个问题,我们再来回答另一个问题:空间转录组联合单细胞转录组目前来说都有哪些方法?
在回答这个问题的之前,其实晨曦的观点一直秉承着“一枪在手”的原则,也就是说如果没有一种方法可以压倒性的取得绝对的优势,或者我们所用的方法已经被证明是错误的或者不正确的,只有在这个时候,晨曦才会考虑去学习其它的方法
Ps:可能这就是懒惰吧QAQ
后期的推文可能也会更新几种不同的联合方法,方便各位小伙伴进一步进行学习,目前来看,常见的空间转录组联合单细胞的方法有:
1. 
AddModuleScore

2. Seurat自带的联合分析方法(FindTransferAnchors and TransferData)
3. MIA
其中第一种和第二种方法都是有R包和相关的函数可以帮助我们完成,但是第三种方法目前来看并没有大佬把这个方法写成函数或者相关的R包,所以晨曦在推文中也只能重点介绍前两种方法
Ps:毕竟没有达到大佬的程度可以手写算法TUT
那么,我们这期推文就首先来介绍第一种方法——AddModuleScore
AddModuleScore_联合单细胞转录组和空间转录组数据
AddModuleScore函数是Seurat包中的一个函数,主要功能就是为基因集来进行打分,提到打分是不是觉得很熟悉,我们经常使用的GSEA或者GSVA不也是进行打分,免疫浸润分析中ssGSEA也是进行打分,那么我就可以简单对打分进行一个简单的定义——通过算法给基因集一个分数
然后要想学习一个函数,首先我们就需要去了解这个函数
#调取帮助文档library(Seurat)?AddModuleScore
晨曦解读
其实如果单纯看描述,大多数人应该不清楚这个函数究竟在做什么,我们从这块能得到的直观信息只有一个——输入数据为:Seurat对象和一个gene list(感兴趣的基因列表)
所以我们直接进入Seurat官网查看函数的源码了解其具体的作用,通过查看源码我们知道了AddModuleScore函数执行的大概思路为先计算基因集中所有基因在不同细胞的平均值,再根据平均值把表达矩阵切割成若干个bin,然后从切割后的每一个bin随机抽取对照基因(基因集外的基因,默认100个)作为背景值。最后所有的目标基因算一个平均值,所有的背景基因算一个平均值,两者相减就是该gene set的score值。因此,在整合不同样本的情况下,即使使用相同基因集为相同细胞打分,也会产生不同的富集评分
上面算法的大致思路如果觉得还是难以理解,其实我们只需要知道这个函数可以把空间转录组和单细胞转录组联合起来达到最基本的会用也是可以的~
简单来说其实就是,我们空间转录组一个sopt可以认为是混合成分,我们想要做的其实就是一个spot里面只有一类细胞或者说我们可以做到把这一个spot内我们感兴趣的基因的表达通过算法给其进行凝练,这让就可以让我们的精度得到提升,那么打分其实就是通过单细胞测序的数据通过算法生成某个细胞亚群的打分,然后把这个打分当作表达量然后映射到空间转录组上,这样我们就相当于把空间转录组的数据精度进行了提高,因为其映射到的表达是来源于单细胞转录组的
那么,理论部分了解这些应该就足够了,我们接下来通过代码来进一步了解。
代码操作
#准备工作library(Seurat)library(tidyHeatmap)library(tidyverse)library(ggplot2)library(ggpubr)
晨曦解读
工欲善其事必先利其器~
#标准scRNA-seq流程获取基因列表#这里我们直接拿外周血的单细胞数据集来作为演示,请注意,我们这里的scRNA-seq数据和Visium数据本身并没有联系,主要是为了演示整体的流程,大家在后续使用的时候,需要确保单细胞和空转数据取自同一个组织library(dplyr)#处理数据library(Seurat)#scRNA-seqlibrary(patchwork)#拼图library(tidyverse)#加载数据pbmc.data <- Read10X(data.dir = "pbmc3k_filtered_gene_bc_matrices/filtered_gene_bc_matrices/hg19/")#创建Seurat对象pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)pbmc#QCpbmc <- PercentageFeatureSet(pbmc, pattern = "^MT-", col.name = "percent.mt")#SCTpbmc <- SCTransform(pbmc, vars.to.regress = "percent.mt", verbose = FALSE)#SCT步骤可以替代这三个步骤NormalizeData(), ScaleData(), and FindVariableFeatures()#标准流程pbmc <- RunPCA(pbmc, verbose = FALSE)pbmc <- FindNeighbors(pbmc, dims = 1:30, verbose = FALSE)pbmc <- FindClusters(pbmc, verbose = FALSE,resolution = 0.4)pbmc <- RunUMAP(pbmc, dims = 1:30, verbose = FALSE)DimPlot(pbmc, label = TRUE) + NoLegend()#可视化查看数据VlnPlot(pbmc, features = c("CD8A", "GZMK", "CCL5", "S100A4", "ANXA1", "CCR7", "ISG15", "CD3D"),       pt.size = 0.2, ncol = 4)FeaturePlot(pbmc, features = c("CD8A", "GZMK", "CCL5", "S100A4", "ANXA1", "CCR7"), pt.size = 0.2,           ncol = 3)#定义细胞亚群new.cluster.ids <- c("Naive CD4 T", "CD14+ Mono", "Memory CD4 T", "B", "CD8 T", "FCGR3A+ Mono", "NK", "DC", "Platelet")names(new.cluster.ids) <- levels(pbmc)pbmc <- RenameIdents(pbmc, new.cluster.ids)DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()#寻找差异基因pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)pbmc.markers %>% group_by(cluster) %>% slice_max(n = 2, order_by = avg_log2FC)
晨曦解读
我们通过scRNA-seq标准流程获得两个数据,分别是pbmc.markers以及pbmc,这两个分别是marker基因以及经过标准流程到细胞注释后的pbmc Seurat对象
#空转数据我们选择的是上期推文的空转数据#为什么没有采用官网的h5格式文件,是因为我们在数据挖掘的时候基本上很少会看到整齐的h5文件#当然了解了杂乱的格式,整齐的格式也就会处理了#准备工作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()#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))
晨曦解读
上面我们就走完了空间转录组的标准流程,那么接下来就是我们的重头戏,我们想要把单细胞数据映射到空间转录组上
注意:本期推文并没有实际生物学意义,因为本身数据来源都是不同的,只作为演示
#AddModuleScore#选择基因列表,这里因为我们并没有感兴趣的基因,所以就常规选择了#如果有感兴趣的基因列表(比如免疫、转移、炎症等等)可以提取感兴趣的基因列表然后映射到空间转录组上table(pbmc.markers$cluster)#获取细胞亚群情况markers <- pbmc.markers %>% dplyr::select(gene,everything()) %>% subset(p_val_adj < 0.05)#基因列表不适宜过多,大概数量可以控制在100-1000左右marker_10 = markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_log2FC) %>% filter(cluster=="NK")#进行AddModuleScore打分只需要两个文件:Genelist以及Visium_Seurat对象#转换成列表gene_list=list(gene=marker_10$gene)#打分heart_ST <-AddModuleScore(heart_ST,features = gene_list,name = "NK")#计算结果保存在[email protected][["NK1"]]head([email protected][["NK1"]])str([email protected][["NK1"]])#空间转录组有3111个barcode ID#这个打分的含义就是这个亚群在每个细胞上的一个分数,可以进行平行比较#至于具体有什么生物学意义就仁者见仁,智者见智了
晨曦解读
这个时候我们就获得了打分的信息,是一个向量的形式储存在Seurat对象中
这个时候我们需要做的就是把我们的打分数据映射到空间位置信息上,完成单细胞和空间转录组的联合分析,我们需要思考我们应该如何进行这一步的操作:
思考:我们有了坐标信息就可以通过ggplot2绘制空间可视化图形,有了打分信息可以把其视为表达信息映射到空间可视化图形上,但是我们为什么可以这么操作呢?
其实这里晨曦的理解是:打分的目的其实就是把表达量或者某种信息抽象为一个具体的数值,然后这个数值就可以进行比较和映射到别的信息上去,拿我们经常见到的免疫浸润来说,我们打分以后可以让其在相同样本下不同的免疫细胞之间可比,但是富集分数究竟是什么其实并不重要,我们只需要知道哪一个免疫细胞类型明显富集就可以,那么我们这里打分也是这个道理,因为空间转录组每一个spot内含有多个细胞(混合结构),这样划分的细胞亚群肯定不是十分准确,那么我就进行打分,每一个spot通过算法进行打分,也就是把混合变成单一的过程,这样就完成了一个“纯化”的过程
这里多说一句,按照这个上面的逻辑来看,万物皆可打分,免疫浸润打分最常见,甚至于这个算法其实只要有一个要计算的对象和一个基因列表就可以做到通过基因列表完成到对象的打分,基因列表为纯粹,对象为混合,纯粹替代混合完成打分,其目的还是只有一个就是可以做到不同条件下的可比
#映射#需要两个信息:打分(充当表达量)和空间位置信息(可视化)#获取空间位置信息(通过调用Seurat中也可以获取)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",]#打分信息我们已经有了,然后接下来就可以映射了metadata <[email protected]#获取空间位置信息metadata <-metadata[rownames(position_ST_ninth),]metadata$x=position_ST_ninth$new_xmetadata$y=position_ST_ninth$new_yggplot(metadata, aes(x = x, y = y)) + geom_point(data = metadata, aes(color= NK1))+ theme_bw()+ xlab(NULL)+ ylab(NULL)+ scale_color_continuous(high="red",low="white")

到这里,我们AddModuleScore进行空间转录组和单细胞转录组联合就结束了,让我们再来回顾一下其中比较重要的内容
1.
AddModuleScore函数需要的输入数据为Seurat对象以及基因列表(这个列表可以自己进行定义,如:感兴趣的亚群、炎症、转移、铁死亡、铜死亡等等)

2. 打分的目的其实就是做到让彼此可比(空间转录组因为每一个spot内的细胞数量都不相等,相互比较其实不是很正确,通过单细胞转录组获取到的单一亚群的marker基因对其进行打分,得到的分数就可以代表这个亚群在每一个spot的情况,这个时候绘制的可视化就比较准确)
3. AddModuleScore函数的核心代码比较简单,重点需要理解逻辑
4. 更多内容欢迎大家阅读下面的参考文献以及参考教程!!
那么,这期推文就给大家介绍到这里啦,欢迎大家在评论区留言,互相学习,共同进步QAQ
我是晨曦,我们下期再见

参考文献以及教程:
A.Seurat包的打分函数AddModuleScore - 简书 (jianshu.com)
B.Seurat的打分函数AddMouduleScore - 简书 (jianshu.com)
C.Tools for Single Cell Genomics • Seurat (satijalab.org)
D.Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq - PubMed (nih.gov)
Ps:回复“空间转录组03”获取本期推文的示例代码以及示例文件哦~
晨曦单细胞文献阅读系列

非肿瘤单细胞分析模板已到位!眼馋单细胞的小伙伴快来看!手把手教你产出第一篇单细胞SCI!

晨曦碎碎念系列传送门(未完待续...)
晨曦单细胞笔记系列传送门
晨曦从零开始学画图系列传送门
晨曦单细胞数据库系列传送门

END

撰文丨晨   曦
排版丨四金兄
主编丨小雪球
欢迎大家关注解螺旋生信频道-挑圈联靠公号~
继续阅读
阅读原文