旧号无故被封，小号再发一次
更多空间转录组文章：
1. 新版10X Visium
- 【10X空间转录组Visium】（一）Space Ranger 1.0.0(更新于20191205）
- 【10X空间转录组Visium】（二）Loupe Browser 4.0.0
- 【10X空间转录组Visium】（三）跑通Visium全流程记录
- 【10X空间转录组Visium】（四）R下游分析的探索性代码示例
- 【10X空间转录组Visium】（五）Visium原理、流程与产品
- 【10X空间转录组Visium】（六）新版Seurat v3.2分析Visium空间转录组结果的代码实操
- 【10X空间转录组Visium】（七）思考新版Seurat V3.2作者在Github给予的回答
2. 旧版Sptial
- 【旧版空间转录组Spatial】（一）ST Spot Detector使用指南
- 【旧版空间转录组Spatial】（二）跑通流程试验记录
- 【旧版空间转录组Spatial】（三）ST Spot Detector实操记录
官网地址：https://support.10xgenomics.com/spatial-gene-expression/software/pipelines/latest/rkit
将Visium数据加载到R中会有所帮助，这些包括：
- 一次查看一个或多个样本的多个基因。
- 一次查看多个样本的特征，包括：Genes, UMIs, Clusters
下面的示例显示如何绘制此信息以构成以下组合的图形：
- Tissue - Total UMI.
- Tissue - Total Gene.
- Tissue - Cluster.
- Tissue - Gene of interest.
探索性分析（个性化分析）：
导入库：
读取h5格式的稀疏矩阵
library(ggplot2)
library(Matrix)
library(rjson)
library(cowplot)
library(RColorBrewer)
library(grid)
library(readbitmap)
library(Seurat)
library(dplyr)
定义函数：定义geom_spatial函数使在ggplot中绘制组织图像变得简单。
geom_spatial <- function(mapping = NULL,
data = NULL,
stat = "identity",
position = "identity",
na.rm = FALSE,
show.legend = NA,
inherit.aes = FALSE,
...) {
GeomCustom <- ggproto(
"GeomCustom",
Geom,
setup_data = function(self, data, params) {
data <- ggproto_parent(Geom, self)$setup_data(data, params)
data
},
draw_group = function(data, panel_scales, coord) {
vp <- grid::viewport(x=data$x, y=data$y)
g <- grid::editGrob(data$grob[[1]], vp=vp)
ggplot2:::ggname("geom_spatial", g)
},
required_aes = c("grob","x","y")
)
layer(
geom = GeomCustom,
mapping = mapping,
data = data,
stat = stat,
position = position,
show.legend = show.legend,
inherit.aes = inherit.aes,
params = list(na.rm = na.rm, ...)
)
}
读取数据：
定义样品
sample_names <- c("Sample1", "Sample2")
sample_names
定义路径
路径应与相应样品名称的顺序相同。
image_paths <- c("/path/to/Sample1-spatial/tissue_lowres_image.png",
"/path/to/Sample2-spatial/tissue_lowres_image.png")
scalefactor_paths <- c("/path/to/Sample1-spatial/scalefactors_json.json",
"/path/to/Sample2-spatial/scalefactors_json.json")
tissue_paths <- c("/path/to/Sample1-spatial/tissue_positions_list.txt",
"/path/to/Sample2-spatial/tissue_positions_list.txt")
cluster_paths <- c("/path/to/Sample1/outs/analysis_csv/clustering/graphclust/clusters.csv",
"/path/to/Sample2/outs/analysis_csv/clustering/graphclust/clusters.csv")
matrix_paths <- c("/path/to/Sample1/outs/filtered_feature_bc_matrix.h5",
"/path/to/Sample2/outs/filtered_feature_bc_matrix.h5")
Read in Down Sampled Images:
确定图像的高度和宽度，以便最终进行正确的绘图。
images_cl <- list()
for (i in 1:length(sample_names)) {
images_cl[[i]] <- read.bitmap(image_paths[i])
}
height <- list()
for (i in 1:length(sample_names)) {
height[[i]] <- data.frame(height = nrow(images_cl[[i]]))
}
height <- bind_rows(height)
width <- list()
for (i in 1:length(sample_names)) {
width[[i]] <- data.frame(width = ncol(images_cl[[i]]))
}
width <- bind_rows(width)
Convert the Images to Grobs:
此步骤提供与ggplot2的兼容性
grobs <- list()
for (i in 1:length(sample_names)) {
grobs[[i]] <- rasterGrob(images_cl[[i]], width=unit(1,"npc"), height=unit(1,"npc"))
}
images_tibble <- tibble(sample=factor(sample_names), grob=grobs)
images_tibble$height <- height$height
images_tibble$width <- width$width
scales <- list()
for (i in 1:length(sample_names)) {
scales[[i]] <- rjson::fromJSON(file = scalefactor_paths[i])
}
Read in Clusters:
clusters <- list()
for (i in 1:length(sample_names)) {
clusters[[i]] <- read.csv(cluster_paths[i])
}
结合聚类和组织信息以轻松绘制：
在这一点上，我们还需要根据 scale factor 调整正在使用的图像的光斑位置。在这种情况下，我们使用的是低分辨率图像，该图像已被Space Ranger调整为600像素（最大尺寸），但也保持了proper aspec ratio。
例如，如果您的图像为12000 x 11000，则图像大小将调整为600 x550。如果您的图像为11000 x 12000，则图像大小将调整为550 x 600。
bcs <- list()
for (i in 1:length(sample_names)) {
bcs[[i]] <- read.csv(tissue_paths[i],col.names=c("barcode","tissue","row","col","imagerow","imagecol"), header = FALSE)
bcs[[i]]$imagerow <- bcs[[i]]$imagerow * scales[[i]]$tissue_lowres_scalef # scale tissue coordinates for lowres image
bcs[[i]]$imagecol <- bcs[[i]]$imagecol * scales[[i]]$tissue_lowres_scalef
bcs[[i]]$tissue <- as.factor(bcs[[i]]$tissue)
bcs[[i]] <- merge(bcs[[i]], clusters[[i]], by.x = "barcode", by.y = "Barcode", all = TRUE)
bcs[[i]]$height <- height$height[i]
bcs[[i]]$width <- width$width[i]
}
names(bcs) <- sample_names
读入矩阵，条形码和基因:
对于最简单的方法，我们正在使用Seurat包读入我们的filtered_feature_bc_matrix.h5。但是，如果您无权访问该程序包，则可以从filtered_feature_be_matrix目录中读取文件，并以条形码作为行名，基因作为列名来重建data.frame。请参见下面的代码示例。
matrix <- list()
for (i in 1:length(sample_names)) {
matrix[[i]] <- as.data.frame(t(Read10X_h5(matrix_paths[i])))
}
可选：如果您希望从filtered_feature_bc_matrix目录中读取而不是使用Seurat。您可以进行上述修改以编写循环以读取这些内容。
matrix_dir = "/path/to/Sample1/outs/filtered_feature_bc_matrix/"
barcode.path <- paste0(matrix_dir, "barcodes.tsv.gz")
features.path <- paste0(matrix_dir, "features.tsv.gz")
matrix.path <- paste0(matrix_dir, "matrix.mtx.gz")
matrix <- t(readMM(file = matrix.path))
feature.names = read.delim(features.path,
header = FALSE,
stringsAsFactors = FALSE)
barcode.names = read.delim(barcode.path,
header = FALSE,
stringsAsFactors = FALSE)
rownames(matrix) = barcode.names$V1
colnames(matrix) = feature.names$V2
可选：如果要分析大量样本，也可以使用doSNOW库并行执行此步骤。
library(doSNOW)
cl <- makeCluster(4)
registerDoSNOW(cl)
i = 1
matrix<- foreach(i=1:length(sample_names), .packages = c("Matrix", "Seurat")) %dopar% {
as.data.frame(t(Read10X_h5(matrix_paths[i])))
}
stopCluster(cl)
Make Summary data.frames:
每个点的总UMI
umi_sum <- list()
for (i in 1:length(sample_names)) {
umi_sum[[i]] <- data.frame(barcode = row.names(matrix[[i]]),
sum_umi = Matrix::rowSums(matrix[[i]]))
}
names(umi_sum) <- sample_names
umi_sum <- bind_rows(umi_sum, .id = "sample")
每个点的基因总数:
gene_sum <- list()
for (i in 1:length(sample_names)) {
gene_sum[[i]] <- data.frame(barcode = row.names(matrix[[i]]),
sum_gene = Matrix::rowSums(matrix[[i]] != 0))
}
names(gene_sum) <- sample_names
gene_sum <- bind_rows(gene_sum, .id = "sample")
合并所有必要数据
In this final data.frame, we have information about your spot barcodes, spot tissue category (in/out), scaled spot row and column position, image size, and summary data.
bcs_merge <- bind_rows(bcs, .id = "sample")
bcs_merge <- merge(bcs_merge,umi_sum, by = c("barcode", "sample"))
bcs_merge <- merge(bcs_merge,gene_sum, by = c("barcode", "sample"))
绘图：
将大量图形组合在一起的最便捷方法是将它们构造成列表并利用cowplot包进行排布
在这里，我们将使用bcs_merge，每个样本针对sample_names进行过滤
我们还将使用给定于每个样本的图像尺寸，以确保我们的绘图具有正确的x和y限制，如下所示。
xlim(0,max(bcs_merge %>%
filter(sample ==sample_names[i]) %>%
select(width)))+
注意：斑点不按比例缩放
定义要绘制的调色板
myPalette <- colorRampPalette(rev(brewer.pal(11, "Spectral")))
每个组织覆盖点的总UMI
plots <- list()
for (i in 1:length(sample_names)) {
plots[[i]] <- bcs_merge %>%
filter(sample ==sample_names[i]) %>%
ggplot(aes(x=imagecol,y=imagerow,fill=sum_umi)) +
geom_spatial(data=images_tibble[i,], aes(grob=grob), x=0.5, y=0.5)+
geom_point(shape = 21, colour = "black", size = 1.75, stroke = 0.5)+
coord_cartesian(expand=FALSE)+
scale_fill_gradientn(colours = myPalette(100))+
xlim(0,max(bcs_merge %>%
filter(sample ==sample_names[i]) %>%
select(width)))+
ylim(max(bcs_merge %>%
filter(sample ==sample_names[i]) %>%
select(height)),0)+
xlab("") +
ylab("") +
ggtitle(sample_names[i])+
labs(fill = "Total UMI")+
theme_set(theme_bw(base_size = 10))+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
axis.text = element_blank(),
axis.ticks = element_blank())
}
plot_grid(plotlist = plots)
image.png
每个组织覆盖点的总基因：
plots <- list()
for (i in 1:length(sample_names)) {
plots[[i]] <- bcs_merge %>%
filter(sample ==sample_names[i]) %>%
ggplot(aes(x=imagecol,y=imagerow,fill=sum_gene)) +
geom_spatial(data=images_tibble[i,], aes(grob=grob), x=0.5, y=0.5)+
geom_point(shape = 21, colour = "black", size = 1.75, stroke = 0.5)+
coord_cartesian(expand=FALSE)+
scale_fill_gradientn(colours = myPalette(100))+
xlim(0,max(bcs_merge %>%
filter(sample ==sample_names[i]) %>%
select(width)))+
ylim(max(bcs_merge %>%
filter(sample ==sample_names[i]) %>%
select(height)),0)+
xlab("") +
ylab("") +
ggtitle(sample_names[i])+
labs(fill = "Total Genes")+
theme_set(theme_bw(base_size = 10))+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
axis.text = element_blank(),
axis.ticks = element_blank())
}
plot_grid(plotlist = plots)
image.png
每个组织覆盖点的聚类分布
plots <- list()
for (i in 1:length(sample_names)) {
plots[[i]] <- bcs_merge %>%
filter(sample ==sample_names[i]) %>%
filter(tissue == "1") %>%
ggplot(aes(x=imagecol,y=imagerow,fill=factor(Cluster))) +
geom_spatial(data=images_tibble[i,], aes(grob=grob), x=0.5, y=0.5)+
geom_point(shape = 21, colour = "black", size = 1.75, stroke = 0.5)+
coord_cartesian(expand=FALSE)+
scale_fill_manual(values = c("#b2df8a","#e41a1c","#377eb8","#4daf4a","#ff7f00","gold", "#a65628", "#999999", "black", "grey", "white", "purple"))+
xlim(0,max(bcs_merge %>%
filter(sample ==sample_names[i]) %>%
select(width)))+
ylim(max(bcs_merge %>%
filter(sample ==sample_names[i]) %>%
select(height)),0)+
xlab("") +
ylab("") +
ggtitle(sample_names[i])+
labs(fill = "Cluster")+
guides(fill = guide_legend(override.aes = list(size=3)))+
theme_set(theme_bw(base_size = 10))+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
axis.text = element_blank(),
axis.ticks = element_blank())
}
plot_grid(plotlist = plots)
image.png
绘制感兴趣的基因:
- 将
bcs_merge的data.frame与包含我们感兴趣的基因的矩阵matrix的子集绑定在一起 - 此处是海马区特异性基因Hpca
- 注意：这是小鼠的一个示例，对于人类来说，基因符号将是HPCA
- 与使用dplyr::select()之类的函数相比，转换为data.table允许极快的取子集方法：
plots <- list()
for (i in 1:length(sample_names)) {
plots[[i]] <- bcs_merge %>%
filter(sample ==sample_names[i]) %>%
bind_cols(as.data.table(matrix[i])[, "Hpca", with=FALSE]) %>%
ggplot(aes(x=imagecol,y=imagerow,fill=Hpca)) +
geom_spatial(data=images_tibble[i,], aes(grob=grob), x=0.5, y=0.5)+
geom_point(shape = 21, colour = "black", size = 1.75, stroke = 0.5)+
coord_cartesian(expand=FALSE)+
scale_fill_gradientn(colours = myPalette(100))+
xlim(0,max(bcs_merge %>%
filter(sample ==sample_names[i]) %>%
select(width)))+
ylim(max(bcs_merge %>%
filter(sample ==sample_names[i]) %>%
select(height)),0)+
xlab("") +
ylab("") +
ggtitle(sample_names[i])+
theme_set(theme_bw(base_size = 10))+
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
axis.text = element_blank(),
axis.ticks = element_blank())
}
plot_grid(plotlist = plots)
image.png
