最近一直在学习单细胞方面的知识，无意中看到别人的帖子在分享如何画堆叠的小提琴图(StackedVlnPlot)，顺便就学习了一下。发现StackedVlnPlot图确实看起来比原始的Seurat函数VlnPlot画的小提琴要清晰不少，可以更美观地有展示多个marker基因在不同cluster中的表达量分布情况。StackedVlnPlot图也经常出现在质量比较好的文章里面。于是乎我就把实现StackedVlnPlot图的方式汇总了一下，写成此帖方便大家学习。画StackedVlnPlot图方法汇总下来共有四种方式，其中第四种是用别人写好的扩展包来实现，最为方便快捷。
现在来明确一下具体需求：如何将Seurat包中VlnPlot函数中画出的marker基因的小提琴图(下面第一个图)实现重排，以到达将同一个marker基因的不同cluster规则地排列在一起(下面第二张图)的目的。
完成这个需求有以下几种实现方法：
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)
最终效果图：
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')
最终效果图：
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')
最终效果图：
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)
最终效果图：
最后
scRNAseq目前来说很火热，也都知道其主要目的之一就是细胞分群，所以展示细胞群的marker基因表达情况就很常见了，这里我总结了4种展示的方式方便大家学习和使用。今天就分享到这里，我会不断地分享学习笔记给大家，也欢迎志同道合的朋友来交流学习。😉😉😉
