单细胞集大成者的seurat包的可视化本质上是ggplot2语法
单细胞数据分析里面最基础的就是降维聚类分群,参考前面的例子:人人都能学会的单细胞聚类分群注释 ,这个大家基本上问题不大了,使用seurat标准流程即可,不过它默认出图并不好看,详见以前我们做的投票:可视化单细胞亚群的标记基因的5个方法,下面的5个基础函数相信大家都是已经烂熟于心了:
VlnPlot(pbmc, features = c("MS4A1", "CD79A")) FeaturePlot(pbmc, features = c("MS4A1", "CD79A")) RidgePlot(pbmc, features = c("MS4A1", "CD79A"), ncol = 1) DotPlot(pbmc, features = unique(features)) + RotatedAxis() DoHeatmap(subset(pbmc, downsample = 100), features = features, size = 3)
前些天我们演示了:各个单细胞亚群的差异基因数量投射到umap图,然后就有留言是不希望看到两个图来说明对应的单细胞亚群的差异基因数量,而是在差异基因数量的FeaturePlot上面就直接标记好单细胞亚群即可。
其实这个需求,跟我们一直以来的马拉松授课常规练习题类似,就是:R语言绘图练习——ggplot2画tSNE的聚类点图(带圈带阴影),其实我们看到的:可视化单细胞亚群的标记基因的5个方法,也就是说seurat包的的5个基础函数,本质上就是针对数据,进行ggplot2语法的可视化。那我们当然是可以自己写代码啦,先回顾一下:各个单细胞亚群的差异基因数量投射到umap图,里面的两个图:
library(SeuratData) #加载seurat数据集
getOption('timeout')
options(timeout=10000)
#InstallData("pbmc3k")
data("pbmc3k")
sce <- pbmc3k.final
library(Seurat)
table(Idents(sce))
DimPlot(sce,label = T)
library(future)
# check the current active plan
plan()
plan("multiprocess", workers = 4)
plan()
sce.markers <- FindAllMarkers(object = sce, only.pos = TRUE,
min.pct = 0.25,
thresh.use = 0.25)
DT::datatable(sce.markers)
pro='markers'
write.csv(sce.markers,file=paste0(pro,'_sce.markers.csv'))
save(sce.markers,file = paste0(pro, 'sce.markers.Rdata'))
pro='markers'
load(file = paste0(pro, 'sce.markers.Rdata'))
sce$celltype = Idents(sce)
df = as.data.frame(table(sce.markers$cluster))
colnames(df) = c('celltype','NofDEGs')
sce$NofDEGs = df[match(sce$celltype,df$celltype),2]
library(patchwork)
DimPlot(sce,label = T,repel = T)+FeaturePlot(sce,'NofDEGs')
接下来我们演示如何使用ggplot2语法自定义绘图:
第一步是模仿FeaturePlot出图
如下所示:
dat = as.data.frame(sce@reductions$umap@cell.embeddings)
dat$cluster=as.character(Idents(sce))
dat$NofDEGs = sce$NofDEGs
head(dat)
library(ggplot2)
library(cowplot)
library(ggrepel)
library(patchwork)
my_FeaturePlot = ggplot()+
geom_point(data=dat,mapping=aes(x=UMAP_1,
y=UMAP_2,
col= NofDEGs),
size=AutoPointSize(data=dat))+
scale_colour_gradient2(low="lightgrey", high = "blue",guide = 'colourbar')+
theme_cowplot() +
theme(
plot.title = element_text(hjust = 0.5),
legend.title=element_blank()
) +ggtitle("NofDEGs")
my_FeaturePlot
library(patchwork)
FeaturePlot(sce,'NofDEGs') + my_FeaturePlot
可以看到,两个图的差异非常小了:
第二步加上各个单细胞亚群的轮廓
这个比较简单, 就是 stat_ellipse 函数,它是正宗的ggplot2语法,如下所示:
table(dat$cluster)
my_FeaturePlot2 = my_FeaturePlot +
stat_ellipse(data=dat,mapping=aes(x= UMAP_1,
y= UMAP_2, group= cluster),
colour ='black',
geom = "path",
linetype = 2, ###圆圈线的类型
size=1, ###圆圈线的粗细
alpha=0.5)
my_FeaturePlot2
出图如下:
最终代码是;
library(tidyverse)
median_df <- dat %>% group_by(cluster) %>%
summarise(median.1 = median(UMAP_1),
median.2 = median(UMAP_2))
head(median_df)
my_FeaturePlot2
my_FeaturePlot2 + ggplot2::geom_text(data = median_df,
aes(median.1,
median.2, label =cluster) )
my_FeaturePlot2 + ggrepel::geom_text_repel(data = median_df,
aes(median.1,
median.2, label =cluster),
min.segment.length = 0)
出图如下所示:
该如何系统性学习ggplot呢
如果你要从ggplot2开始一步步调制成为它这样的美图,需要下很深的功夫,一张统计图就是从数据到几何对象(点、线、条形等)的图形属性(颜色、形状、大小等)的一个映射。
✦ 数据(Data),最基础的是可视化的数据和一系列图形映射(aesthetic mappings),该映射描述了数据中的变量如何映射到可见的图形属性。 ✦ 几何对象(Geometric objects, geoms)代表在图中实际看到的点、线、多边形等。 ✦ 统计转换(Statistical trassformations, stats)是对数据进行某种汇总,例如将数据分组创建直方图,或将一个二维的关系用线性模型进行解释。 ✦ 标度(Scales)是将数据的取值映射到图形空间,例如用颜色、大小或形状来表示不同的取值,展现标度的常见做法是绘制图例和坐标轴。 ✦ 坐标系(Coordinate system, coord)描述数据是如何映射到图形所在的平面,同时提供看图所需的坐标轴和网格线。 ✦ 分面(faceting)如何将数据分解为子集,以及如何对子集作图并展示。 ✦ 主题(theme)控制细节显示,例如字体大小和图形的背景色。
前面我们介绍了绘图小白神包:
另外推荐5个ggplot2资源
ggplot2作者亲自写的书
链接:https://ggplot2-book.org/facet.html
书名是:ggplot2: Elegant Graphics for Data Analysis 作者:Hadley Wickham
This is the online version of work-in-progress 3rd edition of “ggplot2: elegant graphics for data analysis”
虽然这本书有对应的中文译本,但是时间上相对滞后,建议直接看这个在线实时更新版本。
知识点参考卡片(速记表,小抄)
链接:https://ggplot2.tidyverse.org/reference/
sthda网站的ggplot核心图表示例
链接:http://www.sthda.com/english/wiki/ggplot2-essentials
书籍本身提供售卖,价格是17欧元,不过内容都是电子化了,大家直接网页浏览,就是免费的哈!
绘图菜谱
链接:http://www.cookbook-r.com/Graphs/
这个有中文翻译版本,务必直接下单购买,放在书桌旁边随时翻阅。
最后一个是 https://stackoverflow.com/
你会发现,你想实现的各种稀奇古怪的绘图需求,只需要你能使用英文描述出来,就是能找到答案的!
Science: 微生物单细胞、高通量、菌株分辨率,我全都要!| 深度长文
如果你对单细胞转录组研究感兴趣,但又不知道如何入门,也许你可以关注一下下面的课程
看完记得顺手点个“在看”哦!
长按扫码可关注