【10X空間轉(zhuǎn)錄組Visium】（四）R下游分析的探索性代碼示例

舊號(hào)無(wú)故被封，小號(hào)再發(fā)一次

更多空間轉(zhuǎn)錄組文章：

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

官網(wǎng)地址：https://support.10xgenomics.com/spatial-gene-expression/software/pipelines/latest/rkit
將Visium數(shù)據(jù)加載到R中會(huì)有所幫助，這些包括：

  • 一次查看一個(gè)或多個(gè)樣本的多個(gè)基因。
  • 一次查看多個(gè)樣本的特征，包括：Genes, UMIs, Clusters

下面的示例顯示如何繪制此信息以構(gòu)成以下組合的圖形：

  • Tissue - Total UMI.
  • Tissue - Total Gene.
  • Tissue - Cluster.
  • Tissue - Gene of interest.

探索性分析（個(gè)性化分析）：

導(dǎo)入庫(kù)：

讀取h5格式的稀疏矩陣

library(ggplot2)
library(Matrix)
library(rjson)
library(cowplot)
library(RColorBrewer)
library(grid)
library(readbitmap)
library(Seurat)
library(dplyr)

定義函數(shù)：定義geom_spatial函數(shù)使在ggplot中繪制組織圖像變得簡(jiǎn)單。

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, ...)
  )
}

讀取數(shù)據(jù)：
定義樣品

sample_names <- c("Sample1", "Sample2")
sample_names

定義路徑
路徑應(yīng)與相應(yīng)樣品名稱(chēng)的順序相同。

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:
確定圖像的高度和寬度，以便最終進(jìn)行正確的繪圖。

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])
}

結(jié)合聚類(lèi)和組織信息以輕松繪制：
在這一點(diǎn)上，我們還需要根據(jù) scale factor 調(diào)整正在使用的圖像的光斑位置。在這種情況下，我們使用的是低分辨率圖像，該圖像已被Space Ranger調(diào)整為600像素（最大尺寸），但也保持了proper aspec ratio。

例如，如果您的圖像為12000 x 11000，則圖像大小將調(diào)整為600 x550。如果您的圖像為11000 x 12000，則圖像大小將調(diào)整為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

讀入矩陣，條形碼和基因:
對(duì)于最簡(jiǎn)單的方法，我們正在使用Seurat包讀入我們的filtered_feature_bc_matrix.h5。但是，如果您無(wú)權(quán)訪(fǎng)問(wèn)該程序包，則可以從filtered_feature_be_matrix目錄中讀取文件，并以條形碼作為行名，基因作為列名來(lái)重建data.frame。請(qǐng)參見(jiàn)下面的代碼示例。

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。您可以進(jìn)行上述修改以編寫(xiě)循環(huán)以讀取這些內(nèi)容。

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庫(kù)并行執(zhí)行此步驟。

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:
每個(gè)點(diǎn)的總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")

每個(gè)點(diǎn)的基因總數(shù):

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")

合并所有必要數(shù)據(jù)

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"))

繪圖：
將大量圖形組合在一起的最便捷方法是將它們構(gòu)造成列表并利用cowplot包進(jìn)行排布

在這里，我們將使用bcs_merge，每個(gè)樣本針對(duì)sample_names進(jìn)行過(guò)濾

我們還將使用給定于每個(gè)樣本的圖像尺寸，以確保我們的繪圖具有正確的x和y限制，如下所示。

xlim(0,max(bcs_merge %>% 
          filter(sample ==sample_names[i]) %>% 
          select(width)))+

注意：斑點(diǎn)不按比例縮放

定義要繪制的調(diào)色板

myPalette <- colorRampPalette(rev(brewer.pal(11, "Spectral")))

每個(gè)組織覆蓋點(diǎn)的總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

每個(gè)組織覆蓋點(diǎn)的總基因：

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

每個(gè)組織覆蓋點(diǎn)的聚類(lèi)分布

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的子集綁定在一起
  • 此處是海馬區(qū)特異性基因Hpca
  • 注意：這是小鼠的一個(gè)示例，對(duì)于人類(lèi)來(lái)說(shuō)，基因符號(hào)將是HPCA
  • 與使用dplyr::select()之類(lèi)的函數(shù)相比，轉(zhuǎn)換為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
最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡(jiǎn)書(shū)系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容