单细胞分析——深入理解 AnnData 数据结构

Posted 2023-03-21

参考技术A

新手，纯搬运工
【内容出自 https://cloud.tencent.com/developer/article/1817039 】

一、环境准备：
搭建 Python 高效开发环境： Pycharm + Anaconda

二、安装 scanpy

三、AnnData
1、AnnData 介绍与结构
AnnData 是用于存储数据的对象，一般作为 scanpy 的数据存储格式。

下面我们动手构建一个用于创建 AnnoData 的虚拟数据

import numpy as np
import pandas as pd
import anndata as ad
from string import ascii_uppercase

n_obs = 1000

obs = pd.DataFrame()
obs[\'time\'] = np.random.choice([\'day 1\', \'day 2\', \'day 4\', \'day 8\'], n_obs)

var_names = [i*letter for i in range(1, 10) for letter in ascii_uppercase]

n_vars = len(var_names)

var = pd.DataFrame(index=var_names)

X = np.arange(n_obs*n_vars).reshape(n_obs, n_vars)

2、AnnoData 初始化

adata = ad.AnnData(X, obs=obs, var=var, dtype=\'int32\')

print(adata)

3、AnnoData 切片特性
可以看到 AnnData 具有和 dataframe 或 Array 相似的长相，同样具备相似的特性，比如切片：

print(adata.obs_names[:10].tolist())
print(adata.obs_names[-10:].tolist())
print(adata.var_names[:10].tolist())
print(adata.var_names[-10:].tolist())

print(X)

3、AnnoData 的 view 特性

AnnoData 可以实现与 numpy 中的 view 相似的功能。
换句话说就是，我们每次操作 AnnoData 时，并不是再新建一个 AnnoData 来存储数据，而是直接找到已经之前初始化好的 AnnoData 的内存地址，通过内存地址来直接改变 AnnoData 的值。这样做的好处是：

print(adata[:3, \'A\'].X)

adata[:3, \'A\'].X = [0, 0, 0]

print(adata[:5, \'A\'].X)

其实我们在调用 .[] 时，AnnoData已经在内部实现了该操作，也就是说该 view 会成为保存数据的 AnnoData 对象。

但是，如果将 AnnoData 对象的 view 中的一部分赋值，该内容会复制一份并生成新的数据存储对象。

adata_subset = adata[:5, [\'A\', \'B\']]
print(adata_subset)
adata_subset.obs[\'foo\'] = range(5)

可以看到，这时赋值会直接将 AnnoData 对象复制一份。现在 adata_subset 会重新得到一块内存用于存储实际数据，而不再仅仅是对 adata 的内存地址引用。

4、备份到本地

def print_size_in_MB(x):
print(\':.3 MB\'.format(x. sizeof ()/1e6))

print_size_in_MB(adata)

adata.isbacked

adata.filename = \'./write/test.h5ad\'

adata.isbacked

可以看到，我们的 adata 对象已经备份成功，而且就在本地 ‘./write/test.h5ad’ 目录。

前边提到的 view 特性在这里同样适用，我们来看看 adata_subset 是否备份成功。
adata_subset.isbacked
adata_subset.filename = \'./write/adata_subset_test.h5ad\'
adata_subset.isbacked

adata_subset 并没有被启用备份模式，重新设置备份模式。

需要注意的是：备份仅影响数据矩阵 X，所有注释信息都保留在内存中。如果想对全部数据的更改保存，则必须将导出到本地。

5、导出到本地
adata.write("./write/my_results.h5ad")
adata.write_csvs(\'./write/my_results_csvs\', )

6、读取数据
import scanpy as sc
import pandas as pd

adata = sc.read(filename)

anno = pd.read_csv(filename_sample_annotation)

adata.obs[\'cell_groups\'] = anno[\'cell_groups\'] # categorical annotation of type pandas.Categorical

adata.obs[\'time\'] = anno[\'time\'] # numerical annotation of type float

adata.obs = anno

官网： https://anndata.readthedocs.io/en/latest/

将R环境下的Seurat RDS格式数据转化成为到python环境下scanpy的anndata格式

参考技术A

无论是单细胞、空间组还是ATAC的数据，有时由于下游分析的需求或可视化的需求，同时由于python的运算速度的优势，目前越来越多单细胞分析的工具开始在python环境下开发（scanpy/spGCN/scVelo……），但是大家大多都习惯了R的分析环境（Seurat/Harmony/Monocle3……），所以我们经常需要在不同的环境中运行同一个分析对象，这所以涉及到的数据类型的转变就非常关键了。

想直接想找工具将RDS转为python可读数据对象的包，目前还没有……（如果有大佬可以开发一下）。

目前所以从数据本身出发有三种方式，总结自目前网络上一些可行的方法：

1，提取矩阵（稀疏/稠密）和特征信息（metadata），手动构筑 anndata （单细胞分析时python中的一种数据结构，具体了解可以看一下： https://www.jianshu.com/p/9b057e105c42 ，写得挺好）就好。前提是对R的S4对象和python的anndata对象有基础的认识，就可以搞定，这是最本质也是最万能的方法，除了门槛高。

2，存储的时候就注意，不要保存成rds，或者已经这样保存了也无所谓，可以读入再重新存：
（1）存储成 h5ad格式 。Seurat数据写成h5需要借助包 SeuratDisk : https://github.com/mojaveazure/seurat-disk
从R环境下Seurat的对象保存成h5ad的格式：

然后用python的anndata包/scanpy包直接读入就好，因为h5ad本来就是单细胞在跑一python环境中分析最基础的格式，对标R中的seurat对象或sce对象

这个方法可以具体参考： https://www.jianshu.com/p/c438d545f696 他写得更加详细一点

（2）Seurat官方设置 loom格式 也是可行方式之一：
参见： https://www.jianshu.com/p/147da295fc34
将Seurat对象转为loom:

在python环境下读入loom，成为adata:

总而言之，经个人使用和实际操作来说， 第二个存储成为h5ad的方法时最好用的、也最友好 ，除了只能保存一个assay之外，就都没问题，但是如果许多保存多个assay成为anndata中的多个layer，这就的用手动的方法了。