僵小鱼的故事

在那个遥远的上午，百无聊赖地做着项目，用Seurat^[1]画着一个又一个单细胞marker基因小提琴图，“这种图真是够了！”，她在心里想着，“是时候该换一种风格了，最好能够把多个基因的小提琴图规则美观的合并在一起，这样就能直接用在文章里了”。

果然在看过的顶级期刊文章里她找到了一个对胃口的图，“嗯，就是这个了”，想归想，真正要实现的时候发现：非常熟悉的Seurat不能直接产生这种图，“怎么办呢？”

“求助吧”，最近刚加了专业领域的微信群，群里大牛云集，“也许他们能给我点启发”，说着就试着发了自己的想法到群里。

对于群里的新成员，大家都还是比较热心的，有大佬直接指出“scanpy does it”，有的大佬指出“这张图来自Nature肝硬化的文章吧”，“都这么熟悉的么”，僵小鱼在群里回应道。“真的假的？”，吃瓜的我心里想着，然后默默打开搜索引擎，搜索“单细胞肝硬化”的新闻，找到原文“Resolving the fibrotic niche of human liver cirrhosis at single-cell level^[2]”，复制粘贴到sci-hub.tw，下载文章，这些都动作一气呵成，仿佛经历过多年的练习。果然，大佬诚不我欺，我从文章里找到了原图。

群里继续有人在出主意：“Seurat+ppt就能搞定”，“ggplot is able to do everything”，“PPT+1”，“PS食用更佳”，“要么AI，要么AI”，“AI”，“AI ...”。最终僵小鱼选择了AI，以下是成图：

AI才是真的解决一切！

僵小鱼的故事已经告一段落了，我的心里却久久未能平复，作为已经5年没向AI过低头的人，我始终秉承的信念是：python能解决一切。所以在僵小鱼声明“我目前只能接受seurat”后，众人不再理睬的那句“scanpy does it”始终萦绕在我的脑海。

scanpy是处理单细胞数据的python包，基本复现了seurat的主要功能，我曾经测试过，在处理大数据量的单细胞项目时，scanpy的速度和内存真是比seurat友好太多。今天再一次翻看了scanpy教程^[3]，scanpy堆叠小提琴图的风格和Nature文章里的那个图很像，下面是scanpy教程里的小提琴图，在风格上还真是很像：

但是怎样把seurat的对象转换成scanpy能够识别的数据格式呢，这一个是R S3对象，另一个是python的anndata^[4]对象。最初的想法是能不能把seurat对象的矩阵和分群信息导出到文件，再手动构建一个anndata对象，真要做的时候发现面临很多困难。

最终经过在google搜索，毫无意外的发现了同道中人，有相同需求的人在bioinformatics上提问Convert R RNA-seq data object to a Python object^[5]，通过查看这个页面提供的方案，我发现seurat官网^[6]提供了不同单细胞处理软件结果互通的转换方法。

seurat官网提供了seurat对象与SingleCellExperiment、loom、AnnData三种单细胞数据格式相互转换的方法。目前seurat（version 3.1）不支持写入scanpy要求的H5AD文件，所以目前的解决方案是：

1.Seurat对象转换为loom文件2.Scanpy读取loom文件转换为能够操作的anndata对象

要是实现上面的两个简单的步骤还需要安装一些R和python包，需要安装的有以下几个，如果已经安装了，忽略就好：

·R包：seurat^[7]·R包：hdf5r^[8]·R包：loomR^[9]·R包：scater^[10]·python包：scanpy^[11]·python包：loompy^[12]

安装好以上包之后，在R中执行以下代码 ，实现第一步：Seurat对象转换为loom文件

#读入seurat处理后的rds文件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 donesdata.loom$close_all()

接着 在jupyter中执行以下代码 ，实现第二步：Scanpy读取loom文件转换为能够操作的anndata对象

import scanpy as scadata = 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)

初步产生的图和scanpy教程里一样，挑选的marker基因在各个亚群中的表达小提琴图，规则的排布在了一起，基本实现了当初的想法。

再调一下参数，转换X轴和Y轴标签：

ax = sc.pl.stacked_violin(sdata, marker_genes, groupby='ClusterName', rotation=90,swap_axes=True, save="_tmp.png")

到这里和文章里的图只差顺时针旋转90度了，于是查看了scanpy的源代码，发现是一个叫stacked_violin的函数调用的seaborn.violinplot实现的这个小提琴图，那我就copy一下这个函数，修改一些设置，让其由纵向增加小提琴图变为横向增加小提琴图就可以了。转念又一想还是不浪费时间了，“用AI吧，简单操作一下就能搞定了”，不知道从哪里冒出一个声音，真是可恶，那就用代码简单实现一下旋转吧：

from PIL import Imageimport matplotlib.pyplot as pltimg = Image.open('./figures/stacked_violin_tmp.png')img.transpose(Image.ROTATE_270)

多年以后，再次读起这些文字，我一定会回想起此时自己没等看到最后一行便已明白，自己不会再走出这故事。从一个围观的吃瓜群众到讲述这个故事再到成为故事一部分的我，应该很庆幸此时的参与。一切过往自永远至永远不会再重复，因为激发起相同兴趣的情形不会有第二次机会在大地上出现。

---

---

References

[1] Seurat: https://satijalab.org/seurat/pancreas_integration_label_transfer.html
[2] Resolving the fibrotic niche of human liver cirrhosis at single-cell level: https://www.nature.com/articles/s41586-019-1631-3
[3] scanpy教程: https://scanpy-tutorials.readthedocs.io/en/latest/pbmc3k.html
[4] anndata: https://anndata.readthedocs.io/en/stable/api.html#module-anndata
[5] Convert R RNA-seq data object to a Python object: https://bioinformatics.stackexchange.com/questions/3943/convert-r-rna-seq-data-object-to-a-python-object
[6] seurat官网: https://satijalab.org/seurat/v3.1/conversion_vignette.html
[7] seurat: https://satijalab.org/seurat/install.html
[8] hdf5r: https://github.com/hhoeflin/hdf5r
[9] loomR: https://satijalab.org/loomR/loomR_tutorial.html
[10] scater: https://bioconductor.org/packages/release/bioc/html/scater.html
[11] scanpy: https://scanpy.readthedocs.io/en/stable/installation.html
[12] loompy: http://linnarssonlab.org/loompy/installation/index.html

---