前面测试过monocle2/3,但是在画图的时候,尤其是用默认软件画图的时候总是觉得不够好看。比如下面这个paper中的数据(https://pubmed.ncbi.nlm.nih.gov/35896115/),monocle默认的函数就实现不了这个效果了。
所以总是想自己去美化,这时候就需要自己提取画图相关的数据了。
我好多学习也是参考下面这个帖子学习的,很多单细胞相关的学习也是向这个优秀群主学习的(KS科研分享与服务):
https://mp.weixin.qq.com/s/NWvz4vGlLeoEZ1zuM0HmIQ
======加载数据=====
还是用我们经常用的pbmc数据,加载以前测试monocle2已经分析好的cds数据。
https://www.jianshu.com/p/380b9071e244
library(monocle)
library(RColorBrewer)
library(ggplot2)
library(Seurat)
library(ggpubr)
library(tidyverse)
library(dplyr)
library(ggsignif)
library(patchwork)
library(tidydr)
library(ggforce)
library(ggrastr)
library(viridis)
load("cds.rda")
plot_cell_trajectory(cds, color_by = "Pseudotime")
plot_cell_trajectory(cds, color_by = "State")
plot_cell_trajectory(cds, color_by = "cell_type")
和前面美化热图一样,要想提取相应的画图数据,还是看下plot_cell_trajectory函数是如何实现的,是如何准备ggplot的输入数据的,定位到monocle2的plotting.R函数,查看plot_cell_trajectory函数里面的内容,其实最后也是通过ggplot画图实现的。
其实这个图分解来看的话,也就2个要素:散点和轨迹边。plot_cell_trajectory将这2个量分别存在了data_df和edge_df中。
所以data_df相当于把相应的cell,orig.ident, component,State, Pseudotime, celltype等放在了一起,然后一起用ggplot画图。
data_df <- t(reducedDimS(cds)) %>%
as.data.frame() %>%
select_('Component 1' = 1, 'Component 2' = 2) %>%
rownames_to_column("Cells") %>%
mutate(pData(cds)$State,
pData(cds)$Pseudotime,
pData(cds)$orig.ident,
pData(cds)$celltype)
然后改个列名。
colnames(data_df) <- c("cells","Component_1","Component_2","State",
"Pseudotime","orig.ident","celltype")
下面获得图上的轨迹线,也是从这个函数中获得。
reduced_dim_coords <- reducedDimK(cds)
ca_space_df <- Matrix::t(reduced_dim_coords) %>%
as.data.frame() %>%
select_(prin_graph_dim_1 = 1, prin_graph_dim_2 = 2) %>%
mutate(sample_name = rownames(.), sample_state = rownames(.))
dp_mst <- minSpanningTree(cds)
edge_df <- dp_mst %>%
igraph::as_data_frame() %>%
select_(source = "from", target = "to") %>%
left_join(ica_space_df %>% select_(source="sample_name", source_prin_graph_dim_1="prin_graph_dim_1", source_prin_graph_dim_2="prin_graph_dim_2"), by = "source") %>%
left_join(ica_space_df %>% select_(target="sample_name", target_prin_graph_dim_1="prin_graph_dim_1", target_prin_graph_dim_2="prin_graph_dim_2"), by = "target")
ggplot() +
geom_point_rast(data = data_df, aes(x = Component_1,
y = Component_2,
color =Pseudotime))
基本散点图的样子已经出来了,群主用的是 geom_point_rast,用geom_point也是可以的。
再做点美化。
ggplot() +
geom_point_rast(data = data_df, aes(x = Component_1,
y = Component_2,
color =Pseudotime)) +
scale_color_viridis()+
theme_bw()+
theme_dr(arrow = grid::arrow(length = unit(0, "inches")))+#坐标轴主题修改
theme(
panel.background = element_blank(),
panel.border = element_blank(),
panel.grid = element_blank(),
axis.ticks.length = unit(0.8, "lines"),
axis.ticks = element_blank(),
axis.line = element_blank(),
axis.title = element_text(size=15),
)
下面添加轨迹线。
ggplot() +
geom_point_rast(data = data_df, aes(x = Component_1,
y = Component_2,
color =Pseudotime)) +
scale_color_viridis()+
theme_bw()+
theme_dr(arrow = grid::arrow(length = unit(0, "inches")))+#坐标轴主题修改
theme(
panel.background = element_blank(),
panel.border = element_blank(),
panel.grid = element_blank(),
axis.ticks.length = unit(0.8, "lines"),
axis.ticks = element_blank(),
axis.line = element_blank(),
axis.title = element_text(size=15),
) +
geom_segment(aes_string(x="source_prin_graph_dim_1",
y="source_prin_graph_dim_2",
xend="target_prin_graph_dim_1",
yend="target_prin_graph_dim_2"),
size=0.75, linetype="solid", na.rm=TRUE, data=edge_df)
我们自己添加一个箭头:
ggplot() +
geom_point_rast(data = data_df, aes(x = Component_1,
y = Component_2,
color =Pseudotime)) +
scale_color_viridis()+
theme_bw()+
theme_dr(arrow = grid::arrow(length = unit(0, "inches")))+
theme(
panel.background = element_blank(),
panel.border = element_blank(),
panel.grid = element_blank(),
axis.ticks.length = unit(0.8, "lines"),
axis.ticks = element_blank(),
axis.line = element_blank(),
axis.title = element_text(size=15),
) +
geom_segment(aes_string(x="source_prin_graph_dim_1",
y="source_prin_graph_dim_2",
xend="target_prin_graph_dim_1",
yend="target_prin_graph_dim_2"),
size=0.75, linetype="solid", na.rm=TRUE, data=edge_df)+
geom_arc(arrow = arrow(length = unit(0.15, "inches"), type = "closed",angle=30),
aes(x0=0.5,y0=-2,r=3, start=-1.5, end=0.5),
lwd=1.5,color="red")
下面就是学习把饼图添加上去了,可以查看每个state下面不同celltype的比例。
Cellratio <- prop.table(table(data_df$State, data_df$celltype), margin = 2)#计算各组样本不同细胞群比例
Cellratio <- as.data.frame(Cellratio)
colnames(Cellratio) <- c('State',"celltype","Freq")
ggplot() +
geom_point_rast(data = data_df, aes(x = Component_1,
y = Component_2,
color =Pseudotime)) +
scale_color_viridis()+
theme_bw()+
theme_dr(arrow = grid::arrow(length = unit(0, "inches")))+
theme(
panel.background = element_blank(),
panel.border = element_blank(),
panel.grid = element_blank(),
axis.ticks.length = unit(0.8, "lines"),
axis.ticks = element_blank(),
axis.line = element_blank(),
axis.title = element_text(size=15),
) +
geom_segment(aes_string(x="source_prin_graph_dim_1",
y="source_prin_graph_dim_2",
xend="target_prin_graph_dim_1",
yend="target_prin_graph_dim_2"),
size=0.75, linetype="solid", na.rm=TRUE, data=edge_df)+
geom_arc(arrow = arrow(length = unit(0.15, "inches"), type = "closed",angle=30),
aes(x0=0.5,y0=-2,r=3, start=-1.5, end=0.5),
lwd=1.5,color="red")+
geom_arc_bar(data=subset(Cellratio,State=='1'),stat = "pie",
aes(x0=2,y0=3.5,r0=0,r=0.8,amount=Freq,fill=celltype))
其它的箭头和饼图,也可以通过同样的方式添加上去。