本次使用GSE217414数据，包含4个结直肠癌肝转移患者样本数据，这个数据集的数据整理的挺好的，都做了归类。
1.导入
rm(list = ls())
library(Seurat) # 4.4.0
library(SeuratData) # 5.0.2
library(tidyverse)
library(sloop)
library(ggplot2)
library(patchwork)
library(dplyr)
# 10x Visium
sce_s <- Load10X_Spatial(data.dir = "./RawData/GSE6716963/",
filename = "GSM6716963_19G081_filtered_feature_bc_matrix.h5")
# # 高像素
# image <- Read10X_Image(image.dir = "./RawData/GSM5494475/spatial/",
# image.name = "tissue_hires_image.png")
# sce_s <- Load10X_Spatial(data.dir = "./RawData/GSM5494475", image = image,
# filename = "GSM5494475_HNSCC201125T04_filtered_feature_bc_matrix.h5")
2.check数据
head(sce_s[["slice1"]]@coordinates)
# tissue row col imagerow imagecol
# AAACAAGTATCTCCCA-1 1 50 102 1363 3426
# AAACACCAATAACTGC-1 1 59 19 4005 3980
# AAACAGAGCGACTCCT-1 1 14 94 1658 1430
# AAACATTTCCCGGATT-1 1 61 97 1510 4041
# AAACCGGGTAGGTACC-1 1 42 28 3736 3029
# AAACCGTTCGTCCAGG-1 1 52 42 3278 3576
dim(sce_s[["slice1"]]@image)
# [1] 564 600 3
# 564 和 600：表示图像的高度(564像素)和宽度(600像素)。
# 3：表示图像的颜色通道数，即RGB(红、绿、蓝)三通道。每个像素位置有三个值，分别对应于红色、绿色和蓝色的强度，用于合成彩色图像。
head(GetTissueCoordinates(sce_s))
# imagerow imagecol
# AAACAAGTATCTCCCA-1 148.1522 372.3913
# AAACACCAATAACTGC-1 435.3261 432.6087
# AAACAGAGCGACTCCT-1 180.2174 155.4348
# AAACATTTCCCGGATT-1 164.1304 439.2391
# AAACCGGGTAGGTACC-1 406.0869 329.2391
# AAACCGTTCGTCCAGG-1 356.3043 388.6956
head(GetTissueCoordinates(sce_s[["slice1"]]))
# imagerow imagecol
# AAACAAGTATCTCCCA-1 148.1522 372.3913
# AAACACCAATAACTGC-1 435.3261 432.6087
# AAACAGAGCGACTCCT-1 180.2174 155.4348
# AAACATTTCCCGGATT-1 164.1304 439.2391
# AAACCGGGTAGGTACC-1 406.0869 329.2391
# AAACCGTTCGTCCAGG-1 356.3043 388.6956
@coordinates：
这里的 @coordinates 提供的是点在 全分辨率 下（即加载的原始HE图像分辨率）的坐标。因此，这些坐标直接反映在原始高分辨率图像上的点的像素位置。
GetTissueCoordinates()：
GetTissueCoordinates() 返回的坐标通常是基于 缩放后的图像，即使用了低分辨率图像进行缩放或处理后的点的坐标。Seurat 会将这些缩放后的图像应用到显示或处理操作中。
纠正上一个推文的错误，这里的数据不一致是正确的，是因为缩放的原因
3.数据预处理
plot1 <- VlnPlot(sce_s, features = "nCount_Spatial", pt.size = 0.2) +
NoLegend()
plot2 <- SpatialFeaturePlot(sce_s, features = "nCount_Spatial",
pt.size.factor = 2) +
theme(legend.position = "right")
wrap_plots(plot1, plot2)
# SCTransform归一化
library(glmGamPoi)
sce_s <- SCTransform(sce_s, assay = "Spatial",
return.only.var.genes = FALSE, verbose = FALSE)
SpatialDimPlot(sce_s,alpha = 0)
4.基因表达可视化
SpatialFeaturePlot(sce_s, features = c("TP53", "CD3E"))
# pt.size.factor = 2)
library(ggplot2)
plot <- SpatialFeaturePlot(sce_s, features = c("TP53")) +
# alpha = c(0.1, 1)调整透明度 pt.size.factor = 5)
theme(legend.text = element_text(size = 0),
legend.title = element_text(size = 20),
legend.key.size = unit(1, "cm"));plot
ggsave(filename = "./output/images/spatial_tp53.pdf",
height = 9, width = 7)
5.降维聚类可视化
sce_s <- RunPCA(sce_s, assay = "SCT", verbose = FALSE)
sce_s <- FindNeighbors(sce_s, reduction = "pca", dims = 1:30)
sce_s <- FindClusters(sce_s, verbose = FALSE)
sce_s <- RunUMAP(sce_s, reduction = "pca", dims = 1:30)
p1 <- DimPlot(sce_s, reduction = "umap", label = TRUE)
p2 <- SpatialDimPlot(sce_s, label = TRUE, label.size = 3)
#pt.size.factor = 2)
p1 + p2
6.提取特定clusters
SpatialDimPlot(sce_s,
cells.highlight = CellsByIdentities(object = sce_s,
idents = c(2,4,6)),
facet.highlight = TRUE, ncol = 3)
7.交互式绘图
# 设置interactive为True
SpatialDimPlot(sce_s, interactive = TRUE)
SpatialFeaturePlot(sce_s, features = "TP53", interactive = TRUE)
8.空间可变特征识别
一、根据组织内预先注释的解剖区域进行差异表达,比如上面通过FindClusters找到的不同cluster
# 指定比较clust5和6
de_markers <- FindMarkers(sce_s, ident.1 = 5, ident.2 = 6)
# 在空间位置上可视化差异表达基因
SpatialFeaturePlot(object = sce_s,
features = rownames(de_markers)[1:3],
alpha = c(0.1, 1), #pt.size.factor = 5,
ncol = 3)
二、FindSpatiallyVariables中实现的另一种方法是在没有预注释的情况下搜索表现出空间模式的特征
# 默认method = 'markvariogram', 其他有moransi等
# 差异基因从高变基因中获取
sce_s <- FindSpatiallyVariableFeatures(sce_s,
assay = "SCT",
features = VariableFeatures(sce_s)[1:1000],
selection.method = "moransi")
# 提取6个差异基因看一看
top.features <- head(SpatiallyVariableFeatures(sce_s,selection.method = "moransi"),6)
# 异基因可视化/笔者这里自定义了
top.features <- c("CASQ1", "AEBP1", "MME")
SpatialFeaturePlot(sce_s, features = top.features,
ncol = 3, , alpha = c(0.1, 1)) #pt.size.factor = 5
9.对解剖区域(clusters)子集化
修整图片范围
# 根据切片图像坐标(前缀需和Key保持一致),进一步去除多余的spots数据:
datSub <- subset(sce_s, idents = c(1, 2, 3))
SpatialDimPlot(datSub, crop = TRUE, label = TRUE,label.size = 3) # pt.size.factor = 5,
# 获得坐标信息
img<- GetTissueCoordinates(datSub)
head(img)
# 获得cluster的信息
ident<- Idents(datSub)
head(ident)
# 合并切片图像坐标、分群信息；
imgId<- data.frame(img,ident)
head(imgId)
# 赵小明777老师指出显微照片的像素坐标原点在左上角，而图表的原点在左下角，因此需要做垂直翻转y=540-imagerow：
library(ggplot2)
ggplot(img, aes(x=imagecol, y=540-imagerow,color = ident)) +
geom_point(size = 1.6)+
theme_bw()
datSub <- subset(datSub, imagecol > 200 | imagecol < 400, invert = TRUE)
p1 <- SpatialDimPlot(sce_s, crop = TRUE, label = TRUE)
p2 <- SpatialDimPlot(datSub, crop = FALSE, label = TRUE, pt.size.factor = 1, label.size = 3)
p1 + p2
这里笔者根据10X官网和其他老师的代码做了尝试，但是出现了报错，目前还没有找到解决的办法hhh但不影响后面的分析。
10.整合单细胞+空转数据
library(dplyr)
library(qs)
reference <- qread("sc_dataset.qs")
reference <- subset(reference,downsample = 2000)
reference <- SCTransform(reference,ncells = 3000,
verbose = FALSE) %>%
RunPCA(verbose = FALSE) %>%
RunUMAP(dims = 1:30)
DimPlot(reference, group.by = "celltype", label = TRUE)
这次换一个新演员整合两者的数据
anchors <- FindTransferAnchors(reference = reference,
query = sce_s,
normalization.method = "SCT")
predictions.assay <- TransferData(anchorset = anchors,
refdata = reference$celltype,
prediction.assay = TRUE,
weight.reduction = sce_s[["pca"]],
dims = 1:30)
sce_s[["predictions"]] <- predictions.assay
整合进去了可视化-细胞亚群映射
DefaultAssay(sce_s) <- "predictions"
SpatialFeaturePlot(sce_s, features = c("mast cells", "C3"),
pt.size.factor = 1.6, ncol = 2, crop = TRUE)
根据这些预测分数，使用者可以预测位置受空间限制的细胞类型(这里就是能够观察不同细胞亚群在组织上的空间位置)。
使用基于marker point的方法来定义空间高变特征,但使用细胞类型预测分数作为“标记”而不是基因表达。
# 同样需要获取高变异特征，但映射的时候使用的是celltype的分数
sce_s <- FindSpatiallyVariableFeatures(sce_s,
assay = "predictions",
selection.method = "moransi",
features = rownames(sce_s),
r.metric = 5, slot = "data")
# 提取4个差异clsutes看一看
top.clusters <- head(SpatiallyVariableFeatures(sce_s, selection.method = "moransi"), 4)
SpatialPlot(object = sce_s, features = top.clusters, ncol = 2)
这个不展示了，跟上面的图是差不多的。
今天暂时学习到这儿，走完这个流程就已经可以增加基因表达情况和细胞亚群在组织中分布的情况了~ 但没有解决怎么剪裁HE图片hhh
参考资料：
10XGENOMICS：https://www.10xgenomics.com/cn Seurat空转：https://satijalab.org/seurat/articles/spatial_vignette.html 生信技能树：https://mp.weixin.qq.com/s/C57QHpg_uD_YkgrmnvRMPA https://mp.weixin.qq.com/s/xVYCwfJI4MUAuka1W99efQ 生信随笔：https://mp.weixin.qq.com/s/X8eP58iDEx4j7pSyGbze9Q
注：若对内容有疑惑或者有发现明确错误的朋友，请联系后台(欢迎交流)。更多内容可关注公众号：生信方舟
- END -