（视频教程）单细胞marker基因展示值等高线密度图

之前微信VIP群小伙伴问道一个图的做法，是一篇《cell reports》上的文章，展示marker基因的表达的时候，添加上了等高线密度图，之前我们提到过单细胞作图在上面添加等高线（参考：单细胞marker基因可视化的补充---密度图与等高线图）。但是小伙伴发现这个图只有表达的细胞上面添加了等高线密度。所以这里我们修改一下，其实很简单，只需要调整密度的表达量阈值即可。

图片

（reference：Dissecting Cell Lineage Specification and Sex Fate Determination in Gonadal Somatic Cells Using Single-Cell Transcriptomics）

原文作者有很多的函数上传到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)

图片

然后是多个marker的展示：增加ncol参数。设置组图的列数。

#TSNE降维单多个基因KS_plot_density(obj=uterus,                 marker=c("ACTA2", "RGS5","CCL5", "STK17B"),                 dim = "TSNE", size =1, ncol = 2)

图片

最后看一下UMAP降维的数据：​​​​​​​

#UMAP降维marker基因展示KS_plot_density(obj=human_data,                 marker=c("CD3E", "S100A2"),                 dim = "UMAP", size =1, ncol =2)

图片

具体函数如下：

​​​​​​​​​​​​​​

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)          }  
}

这个展示非常完美，觉得分享有用的，点个赞再走呗！