隨著生物學(xué)背景知識(shí)的增加震糖，單細(xì)胞圖譜的可視化直接用10X的Loup或者seurat的Dimplot函數(shù)直接繪制的umap/tsne圖往往很難達(dá)到要求了女揭，這就要求我們提高繪圖技能蚤假。我們都知道ggplot2是一款很好的繪圖R包，甚至可以說(shuō)在語(yǔ)法上擴(kuò)展了R語(yǔ)言本身吧兔。那么勤哗，當(dāng)我們需要繪圖的時(shí)候，自然我們會(huì)想到它及其周邊掩驱。今天我們就主要地看一下ggforce這個(gè)包帶給我們的可能性芒划。
首先，我載入數(shù)據(jù)：
library(Seurat)
library(ggplot2)
library(tidyverse)
pbmc <- readRDS('G:\\Desktop\\Desktop\\RStudio\\single_cell\\filtered_gene_bc_matrices\\hg19pbmc_tutorial.rds')
pbmc
An object of class Seurat
13714 features across 2638 samples within 1 assay
Active assay: RNA (13714 features)
3 dimensional reductions calculated: pca, umap, tsne···
為了提供更多的分群結(jié)果欧穴，我們?cè)倥芤淮蜦indClusters,當(dāng)然你也可自己構(gòu)建分組方式民逼，比如不同的樣本，VDJ不同的克隆型等
pbmc<-FindClusters(pbmc,resolution = c(.4,.8,1.2,1.6,2))
head(pbmc@meta.data)
orig.ident nCount_RNA nFeature_RNA percent.mt RNA_snn_res.0.5
AAACATACAACCAC pbmc3k 2419 779 3.0177759 1
AAACATTGAGCTAC pbmc3k 4903 1352 3.7935958 3
AAACATTGATCAGC pbmc3k 3147 1129 0.8897363 1
AAACCGTGCTTCCG pbmc3k 2639 960 1.7430845 2
AAACCGTGTATGCG pbmc3k 980 521 1.2244898 6
AAACGCACTGGTAC pbmc3k 2163 781 1.6643551 1
seurat_clusters RNA_snn_res.0.4 RNA_snn_res.0.8 RNA_snn_res.1.2
AAACATACAACCAC 1 2 1 5
AAACATTGAGCTAC 0 3 2 2
AAACATTGATCAGC 2 2 1 1
AAACCGTGCTTCCG 3 1 4 4
AAACCGTGTATGCG 8 6 7 8
AAACGCACTGGTAC 1 2 1 1
RNA_snn_res.1.6 RNA_snn_res.2
AAACATACAACCAC 9 1
AAACATTGAGCTAC 2 0
AAACATTGATCAGC 1 2
AAACCGTGCTTCCG 4 3
AAACCGTGTATGCG 8 8
AAACGCACTGGTAC 1 1
為了調(diào)用ggplot2我們把UMAP的坐標(biāo)放到metadata中：
pbmc<-AddMetaData(pbmc,pbmc@reductions$umap@cell.embeddings,col.name = colnames(pbmc@reductions$umap@cell.embeddings))
head(pbmc@meta.data)
讀入一套我珍藏多年的顏色列表：
allcolour=c("#DC143C","#0000FF","#20B2AA","#FFA500","#9370DB","#98FB98","#F08080","#1E90FF","#7CFC00","#FFFF00",
"#808000","#FF00FF","#FA8072","#7B68EE","#9400D3","#800080","#A0522D","#D2B48C","#D2691E","#87CEEB","#40E0D0","#5F9EA0",
"#FF1493","#0000CD","#008B8B","#FFE4B5","#8A2BE2","#228B22","#E9967A","#4682B4","#32CD32","#F0E68C","#FFFFE0","#EE82EE",
"#FF6347","#6A5ACD","#9932CC","#8B008B","#8B4513","#DEB887")
我用ggplot畫一個(gè)帶有標(biāo)簽的umap圖：
class_avg <- pbmc@meta.data %>%
group_by(RNA_snn_res.2) %>%
summarise(
UMAP_1 = median(UMAP_1),
UMAP_2 = median(UMAP_2)
)
umap <- ggplot(pbmc@meta.data ,aes(x=UMAP_1,y=UMAP_2))+
geom_point(aes(color=RNA_snn_res.2))+
scale_color_manual(values = allcolour)+
geom_text(aes(label = RNA_snn_res.2), data = class_avg)+
theme(text=element_text(family="Arial",size=18)) +
theme(panel.background = element_rect(fill='white', colour='black'),
panel.grid=element_blank(), axis.title = element_text(color='black',
family="Arial",size=18),axis.ticks.length = unit(0.4,"lines"),
axis.ticks = element_line(color='black'),
axis.ticks.margin = unit(0.6,"lines"),
axis.line = element_line(colour = "black"),
axis.title.x=element_text(colour='black', size=18),
axis.title.y=element_text(colour='black', size=18),
axis.text=element_text(colour='black',size=18),
legend.title=element_blank(),
legend.text=element_text(family="Arial", size=18),
legend.key=element_blank())+
theme(plot.title = element_text(size=22,colour = "black",face = "bold")) +
guides(colour = guide_legend(override.aes = list(size=5)))
umap
學(xué)會(huì)了嗎?
為了使我們的圖層不要那么復(fù)雜涮帘，還是先畫一個(gè)簡(jiǎn)單的：
umap <- ggplot(pbmc@meta.data ,aes(x=UMAP_1,y=UMAP_2,color=RNA_snn_res.2))+
geom_point()
umap
好奇的我們來(lái)看一下umap這個(gè)圖都有什么：
head(umap$data)
umap$theme
umap$layers
umap$scales
umap$mapping
umap$coordinates
umap$facet
umap$plot_env
umap$labels
然后拼苍，我們請(qǐng)出ggforce這個(gè)包，看看第一次的驚喜调缨。
library(ggforce)
umap + facet_zoom(x = RNA_snn_res.2 == "14")
想要細(xì)致刻畫某個(gè)亞群疮鲫，這不失是一個(gè)方法：
umap + facet_zoom(xlim = c(-15, -10), ylim = c(0, 2.5))
umap + geom_mark_rect(aes(label = RNA_snn_res.2), show.legend = FALSE) +
theme_void()
給每個(gè)群加框加標(biāo)簽，優(yōu)雅：
library(concaveman)
umap +
geom_mark_hull(aes(label = RNA_snn_res.2)) +
theme_void()
可以根據(jù)自己的數(shù)據(jù)格式換一換 ：
umap +
geom_mark_hull(aes(label = RNA_snn_res.2, fill = RNA_snn_res.0.4), show.legend = FALSE) +
theme_void()
如果有需要特別化為一類的可以用背景色來(lái)圈紫乙丁：
umap +
geom_mark_hull(aes(label = RNA_snn_res.2, fill = RNA_snn_res.2), show.legend = FALSE, expand = unit(3, "mm")) +
theme_void()
umap +
geom_mark_hull(aes(label = RNA_snn_res.2, fill = RNA_snn_res.2), show.legend = FALSE, expand = unit(3, "mm")) +
theme_no_axes()
desc <- 'I am a luck dog'
umap +
geom_mark_ellipse(aes(filter = RNA_snn_res.2 == '14', label = '14',
description = desc))
想對(duì)某一亞群做進(jìn)一步的注釋：
你也可以：
umap +
geom_mark_hull(aes(filter = RNA_snn_res.2 == '14', label = '14',
description = desc))
umap +
geom_voronoi_tile(aes(fill = RNA_snn_res.2, group = -1L)) +
geom_voronoi_segment()
Voronoi圖背后的想法是將圖的區(qū)域分割成盡可能多的部分俊犯。與網(wǎng)格或熱圖不同，Voronoi根據(jù)與其他點(diǎn)的接近程度為每個(gè)點(diǎn)繪制自定義形狀伤哺。它返回一個(gè)看起來(lái)像彩色玻璃的圖燕侠。這可以很好地確定每個(gè)區(qū)域內(nèi)的最近點(diǎn)。例如立莉，零售商可以使用它來(lái)查看他們的商店位置所覆蓋的區(qū)域绢彤，并可以幫助他們做出決策，根據(jù)每個(gè)Voronoi形狀的大小來(lái)優(yōu)化他們的位置蜓耻。
其實(shí)這并看不清
umap +
geom_voronoi_tile(aes(fill = RNA_snn_res.2), max.radius = 1,colour = 'black')
看一看出哪些地方的細(xì)胞比較密集茫舶，這一點(diǎn)當(dāng)然需要好的降維結(jié)構(gòu)了，細(xì)胞密集與否分別代表什么刹淌？越密集的區(qū)域細(xì)胞距離越近饶氏，說(shuō)明異質(zhì)性較低。當(dāng)然芦鳍，這和降維結(jié)構(gòu)有關(guān)嚷往。
umap +
geom_voronoi_tile(aes(fill = RNA_snn_res.2), max.radius = .1,colour = 'black')
umap +
geom_mark_hull(aes(fill = RNA_snn_res.2), expand = unit(3, "mm")) +
coord_cartesian(xlim = c(-15, -10), ylim = c(0, 2.5)) +
geom_voronoi_segment()
有種細(xì)胞的感覺(jué)了嗎?
最后，作為附贈(zèng)：
pbmc@meta.data %>%
gather_set_data(6:11) %>%
ggplot(aes(x, id = id, split = y, value = 1)) +
geom_parallel_sets(aes(fill = RNA_snn_res.0.4), show.legend = FALSE, alpha = 0.3) +
geom_parallel_sets_axes(axis.width = 0.1, color = "lightgrey", fill = "white") +
geom_parallel_sets_labels(angle = 0) +
theme_no_axes()
pbmc@meta.data %>%
count(RNA_snn_res.0.4) %>%
ggplot() +
geom_arc_bar(aes(x0 = 0, y0 = 0, r0 = 0.7, r = 1, amount = n, fill = RNA_snn_res.0.4), alpha = 0.3, stat = "pie")
p1 <- pbmc@meta.data %>%
count(RNA_snn_res.0.4) %>%
mutate(focus = ifelse(RNA_snn_res.0.4 == "0", 0.2, 0)) %>%
ggplot()+
geom_arc_bar(aes(x0 = 0, y0 = 0, r0 = 0.7, r = 1, amount = n,
fill = RNA_snn_res.0.4, explode = focus),
alpha = 1, stat = "pie") +
scale_fill_manual(values = allcolour)
- sc-RAN-seq 數(shù)據(jù)分析||Seurat新版教程:Guided Clustering Tutorial
- sc-RAN-seq 數(shù)據(jù)分析||Seurat新版教程: Integrating datasets to learn cell-type specific responses
- sc-RAN-seq 數(shù)據(jù)分析||Seurat新版教程: Using sctransform in Seurat
- 單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析||Seurat新版教程:Differential expression testing
- 單細(xì)胞轉(zhuǎn)錄組 數(shù)據(jù)分析||Seurat新版教程:New data visualization methods in v3.0
- 單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析||Seurat并行策略
- Seurat Weekly NO.0 || 開刊詞
- Seurat Weekly NO.1 || 到底分多少個(gè)群是合適的柠衅？皮仁！
- Seurat Weekly NO.2 || 我該如何取子集
- 你到底想要什么樣的umap/tsne圖？
- scRNA-seq擬時(shí)分析 || Monocle2 踩坑教程
- scRNA-seq數(shù)據(jù)分析 || Monocle3
https://www.r-bloggers.com/the-ggforce-awakens-again/
https://rviews.rstudio.com/2019/09/19/intro-to-ggforce/
