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

舊號無故被封，小號再發(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中會有所幫助克锣，這些包括：

一次查看一個或多個樣本的多個基因。
一次查看多個樣本的特征腔长，包括：Genes, UMIs, Clusters

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

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

探索性分析（個性化分析）：

導(dǎo)入庫：

讀取h5格式的稀疏矩陣

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

定義函數(shù)：定義geom_spatial函數(shù)使在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, ...)
  )
}

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

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

定義路徑
路徑應(yīng)與相應(yī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é)合聚類和組織信息以輕松繪制：
在這一點上您没，我們還需要根據(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

讀入矩陣碌廓，條形碼和基因:
對于最簡單的方法传轰，我們正在使用Seurat包讀入我們的filtered_feature_bc_matrix.h5。但是谷婆，如果您無權(quán)訪問該程序包慨蛙，則可以從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异袄。您可以進(jìn)行上述修改以編寫循環(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庫并行執(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:
每個點的總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")

每個點的基因總數(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，每個樣本針對sample_names進(jìn)行過濾

我們還將使用給定于每個樣本的圖像尺寸讽营，以確保我們的繪圖具有正確的x和y限制虎忌，如下所示僚匆。

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

注意：斑點不按比例縮放

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

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的子集綁定在一起
此處是海馬區(qū)特異性基因Hpca
注意：這是小鼠的一個示例宙址，對于人類來說，基因符號將是HPCA
與使用dplyr::select()之類的函數(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

最后編輯于：2020.03.25 11:02:54

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者

人面猴
序言：七十年代末莉兰，一起剝皮案震驚了整個濱河市挑围，隨后出現(xiàn)的幾起案子，更是在濱河造成了極大的恐慌糖荒，老刑警劉巖杉辙，帶你破解...
沈念sama閱讀 206,126評論 6贊 481
死咒
序言：濱河連續(xù)發(fā)生了三起死亡事件，死亡現(xiàn)場離奇詭異寂嘉，居然都是意外死亡奏瞬，警方通過查閱死者的電腦和手機(jī)枫绅，發(fā)現(xiàn)死者居然都...
沈念sama閱讀 88,254評論 2贊 382
救了他兩次的神仙讓他今天三更去死
文/潘曉璐我一進(jìn)店門，熙熙樓的掌柜王于貴愁眉苦臉地迎上來硼端，“玉大人并淋，你說我怎么就攤上這事≌渥颍” “怎么了县耽？”我有些...
開封第一講書人閱讀 152,445評論 0贊 341
道士緝兇錄：失蹤的賣姜人
文/不壞的土叔我叫張陵，是天一觀的道長镣典。經(jīng)常有香客問我兔毙，道長，這世上最難降的妖魔是什么兄春？我笑而不...
開封第一講書人閱讀 55,185評論 1贊 278
?港島之戀（遺憾婚禮）
正文為了忘掉前任澎剥，我火速辦了婚禮，結(jié)果婚禮上赶舆，老公的妹妹穿的比我還像新娘哑姚。我一直安慰自己，他們只是感情好芜茵，可當(dāng)我...
茶點故事閱讀 64,178評論 5贊 371
惡毒庶女頂嫁案：這布局不是一般人想出來的
文/花漫我一把揭開白布叙量。她就那樣靜靜地躺著，像睡著了一般九串。火紅的嫁衣襯著肌膚如雪绞佩。梳的紋絲不亂的頭發(fā)上，一...
開封第一講書人閱讀 48,970評論 1贊 284
城市分裂傳說
那天猪钮，我揣著相機(jī)與錄音品山，去河邊找鬼。笑死躬贡，一個胖子當(dāng)著我的面吹牛谆奥，可吹牛的內(nèi)容都是我干的。我是一名探鬼主播拂玻，決...
沈念sama閱讀 38,276評論 3贊 399
雙鴛鴦連環(huán)套：你想象不到人心有多黑
文/蒼蘭香墨我猛地睜開眼，長吁一口氣：“原來是場噩夢啊……” “哼宰译！你這毒婦竟也來了檐蚜？” 一聲冷哼從身側(cè)響起，我...
開封第一講書人閱讀 36,927評論 0贊 259
萬榮殺人案實錄
序言：老撾萬榮一對情侶失蹤沿侈，失蹤者是張志新（化名）和其女友劉穎闯第，沒想到半個月后，有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體缀拭，經(jīng)...
沈念sama閱讀 43,400評論 1贊 300
?護(hù)林員之死
正文獨居荒郊野嶺守林人離奇死亡咳短，尸身上長有42處帶血的膿包…… 初始之章·張勛以下內(nèi)容為張勛視角年9月15日...
茶點故事閱讀 35,883評論 2贊 323
?白月光啟示錄
正文我和宋清朗相戀三年填帽，在試婚紗的時候發(fā)現(xiàn)自己被綠了。大學(xué)時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片咙好。...
茶點故事閱讀 37,997評論 1贊 333
活死人
序言：一個原本活蹦亂跳的男人離奇死亡篡腌，死狀恐怖，靈堂內(nèi)的尸體忽然破棺而出勾效，到底是詐尸還是另有隱情嘹悼，我是刑警寧澤，帶...
沈念sama閱讀 33,646評論 4贊 322
?日本核電站爆炸內(nèi)幕
正文年R本政府宣布层宫，位于F島的核電站杨伙，受9級特大地震影響，放射性物質(zhì)發(fā)生泄漏萌腿。R本人自食惡果不足惜限匣，卻給世界環(huán)境...
茶點故事閱讀 39,213評論 3贊 307
男人毒藥：我在死后第九天來索命
文/蒙蒙一、第九天我趴在偏房一處隱蔽的房頂上張望毁菱。院中可真熱鬧米死，春花似錦、人聲如沸鼎俘。這莊子的主人今日做“春日...
開封第一講書人閱讀 30,204評論 0贊 19
一樁弒父案，背后竟有這般陰謀
文/蒼蘭香墨我抬頭看了看天上的太陽贸伐。三九已至勘天，卻和暖如春，著一層夾襖步出監(jiān)牢的瞬間捉邢，已是汗流浹背脯丝。一陣腳步聲響...
開封第一講書人閱讀 31,423評論 1贊 260
情欲美人皮
我被黑心中介騙來泰國打工，沒想到剛下飛機(jī)就差點兒被人妖公主榨干…… 1. 我叫王不留伏伐，地道東北人宠进。一個月前我還...
沈念sama閱讀 45,423評論 2贊 352
代替公主和親
正文我出身青樓，卻偏偏與公主長得像藐翎，于是被迫代替她去往敵國和親材蹬。傳聞我的和親對象是個殘疾皇子，可洞房花燭夜當(dāng)晚...
茶點故事閱讀 42,722評論 2贊 345

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

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

1. 新版10X Visium

2. 舊版Sptial

推薦閱讀更多精彩內(nèi)容