美文网首页单细胞测序
【10X空间转录组Visium】（四）R下游分析的探索性代码示例

【10X空间转录组Visium】（四）R下游分析的探索性代码示例

作者: Geekero | 来源:发表于2020-03-25 10:20 被阅读0次

旧号无故被封，小号再发一次

更多空间转录组文章：

1. 新版10X Visium
2. 旧版Sptial

官网地址：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

相关文章

网友评论

    本文标题：【10X空间转录组Visium】（四）R下游分析的探索性代码示例

    本文链接：https://www.haomeiwen.com/subject/owhtuhtx.html

    延伸阅读

    深度阅读

    • 您也可以注册成为美文阅读网的作者，发表您的原创作品、分享您的心情！

    栏目导航

    热点阅读

    单细胞测序

    关于我们|服务条款|联系我们|【10X空间转录组Visium】（四）R下游分析的探索性代码示例|投稿指南|网站地图|RSS订阅|排版工具|手机版

    提供经典美文摘抄,优美散文欣赏,现代诗歌精选,短篇小说,心情随笔,表白情书范文,故事会在线阅读欣赏

    Copyright © 2014-2023 Haomeiwen.com All Rights Reserved. 好美文阅读网 版权所有

    备案信息：桂公网安备 45052102000051号 · 桂ICP备13007215号-3

    本站所收录作品、热点评论等信息部分来源互联网，目的只是为了系统归纳学习和传递资讯

    所有作品版权归原创作者所有，与本站立场无关，如不慎侵犯了你的权益，请联系我们告知，我们将做删除处理！