之前一个图的做法，是一篇《cell reports》上的文章，展示marker基因的表达的时候，添加上了等高线密度图，（参考：单细胞marker基因可视化的补充---密度图与等高线图）。但是这个图只有表达的细胞上面添加了等高线密度。所以这里我们修改一下，其实很简单，只需要调整密度的表达量阈值即可。
image.png
原文作者有很多的函数上传到github，自行下载原文查看，还是有很多好代码的，但是它的代码并不能满足我们复杂的作图。而且数据完全不一样。不过它图的修饰和颜色可以参考。所以这里我们重新写了一个作图函数，可以展示UMAP、TSNE降维的单细胞seurat对象单个或者多个marker基因的展示。
接下来我们看看具体的效果。演示视频参见B站： https://www.bilibili.com/video/BV1YF41117Py/?spm_id_from=333.999.0.0&vd_source=05b5479545ba945a8f5d7b2e7160ea34首先做一下TSNE降维的单个marker展示：
#TSNE降维单个基因
KS_plot_density(obj=uterus,
marker= "CCL5",
dim = "TSNE", size =1)
image.png
然后是多个marker的展示：增加ncol参数。设置组图的列数。
#TSNE降维单多个基因
KS_plot_density(obj=uterus,
marker=c("ACTA2", "RGS5","CCL5", "STK17B"),
dim = "TSNE", size =1, ncol = 2)
image.png
最后看一下UMAP降维的数据：
#UMAP降维marker基因展示
KS_plot_density(obj=human_data,
marker=c("CD3E", "S100A2"),
dim = "UMAP", size =1, ncol =2)
image.png
具体函数如下：
KS_plot_density <- function(obj,#单细胞seurat对象
marker,#需要展示的marker基因
dim=c("TSNE","UMAP"),
size,#点的大小
ncol=NULL#多个marker基因表达图排序
){
require(ggplot2)
require(ggrastr)
require(Seurat)
cold <- colorRampPalette(c('#f7fcf0','#41b6c4','#253494'))
warm <- colorRampPalette(c('#ffffb2','#fecc5c','#e31a1c'))
mypalette <- c(rev(cold(11)), warm(10))
if(dim=="TSNE"){
xtitle = "tSNE1"
ytitle = "tSNE2"
}
if(dim=="UMAP"){
xtitle = "UMAP1"
ytitle = "UMAP2"
}
if(length(marker)==1){
plot <- FeaturePlot(obj, features = marker)
data <- plot$data
if(dim=="TSNE"){
colnames(data)<- c("x","y","ident","gene")
}
if(dim=="UMAP"){
colnames(data)<- c("x","y","ident","gene")
}
p <- ggplot(data, aes(x, y)) +
geom_point_rast(shape = 21, stroke=0.25,
aes(colour=gene,
fill=gene), size = size) +
geom_density_2d(data=data[data$gene>0,],
aes(x=x, y=y),
bins = 5, colour="black") +
scale_fill_gradientn(colours = mypalette)+
scale_colour_gradientn(colours = mypalette)+
theme_bw()+ggtitle(marker)+
labs(x=xtitle, y=ytitle)+
theme(
plot.title = element_text(size=12, face="bold.italic", hjust = 0),
axis.text=element_text(size=8, colour = "black"),
axis.title=element_text(size=12),
legend.text = element_text(size =10),
legend.title=element_blank(),
aspect.ratio=1,
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()
)
return(p)
}else{
gene_list <- list()
for (i in 1:length(marker)) {
plot <- FeaturePlot(obj, features = marker[i])
data <- plot$data
if(dim=="TSNE"){
colnames(data)<- c("x","y","ident","gene")
}
if(dim=="UMAP"){
colnames(data)<- c("x","y","ident","gene")
}
gene_list[[i]] <- data
names(gene_list) <- marker[i]
}
plot_list <- list()
for (i in 1:length(marker)) {
p <- ggplot(gene_list[[i]], aes(x, y)) +
geom_point_rast(shape = 21, stroke=0.25,
aes(colour=gene,
fill=gene), size = size) +
geom_density_2d(data=gene_list[[i]][gene_list[[i]]$gene>0,],
aes(x=x, y=y),
bins = 5, colour="black") +
scale_fill_gradientn(colours = mypalette)+
scale_colour_gradientn(colours = mypalette)+
theme_bw()+ggtitle(marker[i])+
labs(x=xtitle, y=ytitle)+
theme(
plot.title = element_text(size=12, face="bold.italic", hjust = 0),
axis.text=element_text(size=8, colour = "black"),
axis.title=element_text(size=12),
legend.text = element_text(size =10),
legend.title=element_blank(),
aspect.ratio=1,
panel.grid.major = element_blank(),
panel.grid.minor = element_blank()
)
plot_list[[i]] <- p
}
Seurat::CombinePlots(plot_list, ncol = ncol)
}
}
这个展示非常完美，觉得分享有用的，点个赞再走呗
网友评论