scRNAseq || 堆叠小提琴图(StackedVlnPlot)的4种实现方式

生信店小二 · 简书 · · 2020-07-24 19:51

最近一直在学习单细胞方面的知识，无意中看到别人的帖子在分享如何画堆叠的小提琴图(StackedVlnPlot)，顺便就学习了一下。发现StackedVlnPlot图确实看起来比原始的Seurat函数VlnPlot画的小提琴要清晰不少，可以更美观地有展示多个marker基因在不同cluster中的表达量分布情况。StackedVlnPlot图也经常出现在质量比较好的文章里面。于是乎我就把实现StackedVlnPlot图的方式汇总了一下，写成此帖方便大家学习。画StackedVlnPlot图方法汇总下来共有四种方式，其中第四种是用别人写好的扩展包来实现，最为方便快捷。
现在来明确一下具体需求：如何将Seurat包中VlnPlot函数中画出的marker基因的小提琴图(下面第一个图)实现重排，以到达将同一个marker基因的不同cluster规则地排列在一起(下面第二张图)的目的。

Seurat VInPlot

StackedVlnPlot

完成这个需求有以下几种实现方法：
1、通过Seuart->scanpy来实现，第一张是Seurat包VlnPlot函数画的图，第二张是scanpy中stacked_violin函数画的图，那么现在问题就变成为Seurat对象到scanpy对象的转换
2、用R原生函数实现StackedVlnPlot的方法
3、使用基于scanpy包衍生的scanyuan包来说实现
4、使用R包MySeuratWrappers来实现

下面来看看每种方式的具体实现过程
1、第1种实现方法的具体代码
该方法来自 seurat结果转为scanpy可处理对象。
此过程用到的包，R：seurat、hdf5r、loomR、scater；python：scanpy、loompy

#读入seurat处理后的rds文件，以下是R代码
library(Seurat)
library(loomR)
sdata <- readRDS(file = "/Users/yuanzan/Desktop/tmp/seurat_project.rds")

# seurat对象转换为loop文件
sdata.loom <- as.loom(x = sdata, filename = "/Users/yuanzan/Desktop/tmp/sdata.loom", verbose = FALSE)
# Always remember to close loom files when done
sdata.loom$close_all()

#Scanpy读取loom文件转换为能够操作的anndata对象，以下是python代码
import scanpy as sc
adata = sc.read_loom("/Users/yuanzan/Desktop/tmp/sdata.loom", sparse=True, cleanup=False, X_name='spliced', obs_names='CellID', var_names='Gene', dtype='float32')
marker_genes = ['Stfa1', 'Ngp', 'Ccl5', 'Ccl4', 'BC100530', 'Gzma', 'Gata2', 'Cd74']
ax = sc.pl.stacked_violin(adata, marker_genes, groupby='ClusterName', rotation=90)

#转换X轴和Y轴标签：
ax = sc.pl.stacked_violin(sdata, marker_genes, groupby='ClusterName', rotation=90,swap_axes=True, save="_tmp.png")

#旋转图片角度
from PIL import Image
import matplotlib.pyplot as plt
img = Image.open('./figures/stacked_violin_tmp.png')
img.transpose(Image.ROTATE_270)

最终效果图：

stacked_violin

2、第2种实现方法的具体代码
该方法来自 R函数实现单细胞StackedVlnPlot 。

library(Seurat)
library(ggplot2)
modify_vlnplot <- function(obj, feature, pt.size = 0, plot.margin = unit(c(-0.75, 0, -0.75, 0), "cm"),...) {
       p <- VlnPlot(obj, features = feature, pt.size = pt.size, ... ) +
               xlab("") + ylab(feature) + ggtitle("") +
               theme(legend.position = "none",
               axis.text.x = element_blank(),
               axis.text.y = element_blank(),
               axis.ticks.x = element_blank(),
               axis.ticks.y = element_line(),
               axis.title.y = element_text(size = rel(1), angle = 0, vjust = 0.5),
               plot.margin = plot.margin )
       return(p)
}

## main function
StackedVlnPlot <- function(obj, features, pt.size = 0, plot.margin = unit(c(-0.75, 0, -0.75, 0), "cm"), ...) {
       plot_list <- purrr::map(features, function(x) modify_vlnplot(obj = obj,feature = x, ...))
            plot_list[[length(plot_list)]]<- plot_list[[length(plot_list)]] +
            theme(axis.text.x=element_text(), axis.ticks.x = element_line())
       p <- patchwork::wrap_plots(plotlist = plot_list, ncol = 1)
       return(p)
}

StackedVlnPlot(sdata, c('Retnlg', 'Pygl', 'Anxa1', 'Igf1r', 'Stfa2l1'), pt.size=0, cols=my36colors)

#配色方案
my36colors <- c('#E5D2DD', '#53A85F', '#F1BB72', '#F3B1A0', '#D6E7A3', '#57C3F3', '#476D87',
         '#E95C59', '#E59CC4', '#AB3282', '#23452F', '#BD956A', '#8C549C', '#585658',
         '#9FA3A8', '#E0D4CA', '#5F3D69', '#C5DEBA', '#58A4C3', '#E4C755', '#F7F398',
         '#AA9A59', '#E63863', '#E39A35', '#C1E6F3', '#6778AE', '#91D0BE', '#B53E2B',
         '#712820', '#DCC1DD', '#CCE0F5', '#CCC9E6', '#625D9E', '#68A180', '#3A6963',
         '#968175')

最终效果图：

StackedVlnPlot

3、第3种实现方法的具体代码
该方法来自生信工程师的自我修养。

import scanpy as sc
import scanyuan as scy

# 此处示例为读取loom文件，也可以是其他scanpy支持的数据格式
adata = sc.read_loom("/Users/yuanzan/Desktop/tmp/sdata.loom", sparse=True, cleanup=False, X_name='spliced', obs_names='CellID', var_names='Gene', dtype='float32')
marker_genes = ['Stfa1', 'Ngp', 'Ccl5', 'Ccl4', 'BC100530', 'Gzma', 'Gata2', 'Cd74']
ax = scy.stacked_violin_t(adata, marker_genes, figsize=[8,4], groupby='ClusterName')

最终效果图：

StackedVlnPlot

4、第4种实现方法的具体代码
目前已经有人扩展了Seuart包的功能，弄了一个扩展包MySeuratWrappers，该包的github源码链接为： https://github.com/lyc-1995/MySeuratWrappers ，下面我用10x提供的测试数据集pbmc_1k_v3( http://cf.10xgenomics.com/samples/cell-exp/3.0.0/pbmc_1k_v3/pbmc_1k_v3_fastqs.tar )来完成这个扩展包的实现StackedVlnPlot图的功能。

#安装MySeuratWrappers
>install.packages("remotes")
>remotes::install_github("lyc-1995/MySeuratWrappers")
>library(MySeuratWrappers)

#前面比对等过程忽略，直接从seurat对象开始，例如我已经保存了比对并过来好的表达矩阵到R镜像文件，这里直接load使用
>load('pbmc_umi_counts.RData') #加载矩阵到环境
>str(mat)    #查看矩阵结构
Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
  ..@ i       : int [1:3157185] 15 29 61 72 81 87 91 92 94 97 ...
  ..@ p       : int [1:1141] 0 3197 5709 7791 9346 11851 14464 16679 20764 25434 ...
  ..@ Dim     : int [1:2] 25744 1140
  ..@ Dimnames:List of 2
  .. ..$ : chr [1:25744] "7SK" "A1BG" "A1CF" "A2M" ...
  .. ..$ : chr [1:1140] "AAACCCAAGGAGAGTA" "AAACGCTTCAGCCCAG" "AAAGAACAGACGACTG" "AAAGAACCAATGGCAG" ...
  ..@ x       : num [1:3157185] 2 1 1 1 1 1 2 1 3 1 ...
  ..@ factors : list()

>mat[1:5,1:4]   #查看矩阵前5行5列
5 x 5 sparse Matrix of class "dgCMatrix"
        AAACCCAAGGAGAGTA AAACGCTTCAGCCCAG AAAGAACAGACGACTG AAAGAACCAATGGCAG
7SK                    .                2                .                .
A1BG                   1                .                1                .
A1CF                   .                .                .                .
A2M                    2                .                .                2
A2M-AS1                .                .                3                .

#下面是Seurat的一系列操作
obj <- CreateSeuratObject(counts=mat,project='pbmc')
obj <- NormalizeData(obj, normalization.method = "LogNormalize", scale.factor = 10000)
obj <- FindVariableFeatures(obj, selection.method = "vst", nfeatures = 2000)
obj <- ScaleData(obj)
obj <- RunPCA(obj, features = VariableFeatures(object = obj))
obj <- FindNeighbors(obj, dims = 1:10, k.param = 20)
obj <- FindClusters(obj, resolution = 0.6)
find_markers <- FindAllMarkers(obj, test.use = "wilcox", only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
plot_markers <- sort(c('IL7R', 'CD79A', 'MS4A1', 'CD8A', 'CD8B', 'LYZ')
VlnPlot(obj, features = maker,stacked=T,pt.size=0.5)

最终效果图：

StackedVlnPlot

最后

今天就分享到这里，我会不断地分享学习笔记给大家，也欢迎志同道合的朋友来交流学习。各位看官们帮忙点个赞吧！！！😉😉😉

原文地址：访问原文地址
快照地址：访问文章快照

分享到微博