AnnData 教程 1:开始使用 anndata

写在前面

学习一个软件最好的方法就是啃它的官方文档。本着自己学习、分享他人的态度,分享官方文档的中文教程。软件可能随时更新,建议配合官方文档一起阅读。


::: block-1

目录

  • 1 初始化 AnnData
  • 2 添加对齐的 metadata
  • 3 观测/变量-水平矩阵
  • 4 非结构化 metadata
  • 5 层
  • 6 转换为 DataFrames
  • 7 将结果写入磁盘
  • 8 总结介绍
  • 9 视图和副本
  • 10 部分读取大数据
    :::

官网教程:https://anndata.readthedocs.io/en/latest/tutorials/notebooks/getting-started.html

在本教程中,我们介绍中心对象 AnnData(“Annotated Data”)的基本属性。

AnnData 专为类似矩阵的数据而设计。这意味着我们有 n 个观测值(observations),每个观测值可以表示为 d 维向量,其中每个维度对应于一个变量或特征。此 n×d 矩阵的行和列在索引的意义上是特殊的。

例如,在 scRNA-seq 数据中,每行对应于具有一个 barcode 的细胞,每列对应于具有 gene id 的基因。此外,对于每个细胞和每个基因,我们可能有额外的 metadata,如(1)每个细胞的供体信息,或(2)每个基因的替代基因符号。最后,我们可能有其他 unstructured metadata,如调色板,用于绘图。在不深入研究每一种基于 Python 的奇特数据结构的情况下,我们认为时至今日,还不存在其他真正的替代方案:

  • 处理稀疏性
  • 处理非结构化数据
  • 处理观测-水平和特征-水平的 metadata
  • 用户友好
import numpy as np
import pandas as pd
import anndata as ad
from scipy.sparse import csr_matrix
print(ad.__version__)
## 0.8.0

1 初始化 AnnData

让我们从构建一个包含一些稀疏计数信息的基本 AnnData 对象开始,可能代表 gene expression counts。

counts = csr_matrix(np.random.poisson(1, size=(100, 2000)), dtype=np.float32)
adata = ad.AnnData(counts)
adata
## AnnData object with n_obs × n_vars = 100 × 2000

我们看到 AnnData 提供了具有数据摘要统计的表示形式。我们传递的初始数据可以使用 adata.X 作为稀疏矩阵进行访问。

adata.X
## <100x2000 sparse matrix of type '<class 'numpy.float32'>'
##	with 126543 stored elements in Compressed Sparse Row format>

现在,我们使用 .obs_names.var_namesobsvar 轴提供索引。

adata.obs_names = [f"Cell_{i:d}" for i in range(adata.n_obs)]
adata.var_names = [f"Gene_{i:d}" for i in range(adata.n_vars)]
print(adata.obs_names[:10])
## Index(['Cell_0', 'Cell_1', 'Cell_2', 'Cell_3', 'Cell_4', 'Cell_5', 'Cell_6','Cell_7', 'Cell_8', 'Cell_9'], dtype='object')

AnnData 取子集

这些索引值可用于对 AnnData 进行取子集,从而提供 AnnData 对象的视图。我们可以想象这对于将 AnnData 选取特定的细胞类型或感兴趣的基因模块的子集很有用。AnnData 取子集的规则与 Pandas DataFrame 的规则非常相似。您可以使用 obs_names/var_names、boolean masks, or cell index integer。

adata[["Cell_1", "Cell_10"], ["Gene_5", "Gene_1900"]]
## View of AnnData object with n_obs × n_vars = 2 × 2

2 添加对齐的 metadata

观测/变量水平

所以我们有了对象的核心,现在我们想给 observation 和 variable levels 添加 metadata。这对于 AnnData 来说非常简单,adata.obsadata.var 都是 Pandas DataFrames。

ct = np.random.choice(["B", "T", "Monocyte"], size=(adata.n_obs,))
adata.obs["cell_type"] = pd.Categorical(ct)  # Categoricals are preferred for efficiency
adata.obs

## cell_type
## Cell_0	B
## Cell_1	B
## Cell_2	Monocyte
## Cell_3	Monocyte
## Cell_4	Monocyte
## ...	...
## Cell_95	B
## Cell_96	B
## Cell_97	B
## Cell_98	Monocyte
## Cell_99	Monocyte
## 100 rows × 1 columns

我们现在还可以看到 AnnData 表示已经更新:

adata
## AnnData object with n_obs × n_vars = 100 × 2000
##    obs: 'cell_type'

使用 metadata 取子集

我们还可以使用这些随机生成的细胞类型对 AnnData 进行取子集:

bdata = adata[adata.obs.cell_type == "B"]
bdata
## View of AnnData object with n_obs × n_vars = 36 × 2000
##    obs: 'cell_type'

3 观测/变量-水平矩阵

我们也可能在任一 level 都有 metadata,这些 metadata 具有许多维度,例如数据的 UMAP embedding。对于这种类型的 metadata,AnnData 具有 .obsm/.varm 属性。我们使用 keys 来识别我们插入的不同矩阵。.obsm/.varm 的限制是 .obsm 矩阵的长度必须等于 .n_obs(observations的数目),而 .varm 矩阵的长度必须等于 .n_vars。它们可以各自独立地具有不同数量的维度。

让我们从一个随机生成的矩阵开始,我们可以将其解释为我们想要存储的数据的 UMAP embedding,以及一些随机的 gene-level metadata:

adata.obsm["X_umap"] = np.random.normal(0, 1, size=(adata.n_obs, 2))
adata.varm["gene_stuff"] = np.random.normal(0, 1, size=(adata.n_vars, 5))
adata.obsm
## AxisArrays with keys: X_umap

同样,可以看到 AnnData 表示已经更新:

adata
## AnnData object with n_obs × n_vars = 100 × 2000
##    obs: 'cell_type'
##    obsm: 'X_umap'
##    varm: 'gene_stuff'

关于 .obsm/.varm 的更多说明:

  1. “array-like” metadata 可以源自 Pandas DataFrame、scipy sparse matrix 或 numpy dense array。
  2. 使用 scanpy 时,它们的 values(columns) 不容易绘制,取而代之的是 .obs 中的项目很容易绘制,例如 UMAP 图。

4 非结构化 metadata

AnnData 有 .uns,它允许任何 unstructured metadata。这可以是任何东西,例如包含一些对我们的数据分析有用的一般信息的 list 或 dictionary。

adata.uns["random"] = [1, 2, 3]
adata.uns
## OverloadedDict, wrapping:
## 	OrderedDict([('random', [1, 2, 3])])
## With overloaded keys:
##	['neighbors'].

5 层

最后,我们可能有不同形式的原始核心数据,可能一种是 normalized 的,另一种不是。这些可以存储在 AnnData 的不同 layers 中。例如,让我们对原始数据进行 log transform 并将其存储在一个 layer 中:

adata.layers["log_transformed"] = np.log1p(adata.X)
adata
## AnnData object with n_obs × n_vars = 100 × 2000
##     obs: 'cell_type'
##     uns: 'random'
##     obsm: 'X_umap'
##     varm: 'gene_stuff'
##     layers: 'log_transformed'

6 转换为 DataFrames

我们还可以访问 AnnData 从其中一个 layer 返回一个 DataFrame:

adata.to_df(layer="log_transformed")
## Gene_0	Gene_1	Gene_2	Gene_3	Gene_4	Gene_5	Gene_6	Gene_7	Gene_8	Gene_9	...	Gene_1990	Gene_1991	Gene_1992	Gene_1993	Gene_1994	Gene_1995	Gene_1996	Gene_1997	Gene_1998	Gene_1999
## Cell_0	0.693147	0.000000	0.693147	1.098612	0.693147	0.693147	0.693147	0.000000	0.000000	0.693147	...	0.000000	1.098612	0.000000	0.693147	0.693147	0.693147	1.098612	0.693147	0.693147	0.693147
## Cell_1	0.000000	1.098612	0.000000	1.098612	0.693147	0.693147	0.693147	0.000000	0.693147	0.693147	...	0.000000	0.693147	0.000000	0.693147	0.000000	1.098612	1.386294	0.000000	0.000000	1.098612
## Cell_2	1.098612	0.000000	0.000000	1.098612	0.693147	0.693147	1.098612	0.000000	0.000000	0.693147	...	1.098612	0.693147	1.098612	1.098612	0.693147	0.000000	1.098612	0.693147	0.693147	0.693147
## Cell_3	0.000000	1.098612	0.693147	1.386294	0.693147	1.098612	0.693147	1.098612	0.000000	1.098612	...	0.000000	0.000000	0.693147	1.386294	1.098612	0.693147	0.693147	0.000000	0.000000	0.693147
## Cell_4	0.693147	1.098612	0.693147	1.386294	0.693147	1.386294	0.693147	1.386294	1.098612	0.000000	...	1.098612	0.000000	1.386294	1.098612	0.693147	0.000000	0.000000	0.693147	0.693147	1.098612
## ...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
## Cell_95	1.609438	1.386294	0.693147	0.693147	1.609438	0.693147	0.000000	1.609438	0.693147	1.386294	...	0.000000	1.386294	0.693147	0.693147	0.000000	0.000000	1.386294	0.693147	1.098612	0.693147
## Cell_96	0.693147	0.693147	0.000000	0.000000	1.098612	0.000000	1.098612	0.693147	0.000000	0.000000	...	0.000000	1.098612	1.386294	1.098612	0.000000	0.693147	1.609438	0.693147	0.693147	0.000000
## Cell_97	1.098612	0.693147	0.000000	0.000000	0.000000	0.000000	0.000000	1.098612	0.693147	0.000000	...	0.693147	0.693147	0.000000	0.693147	0.000000	1.098612	0.000000	0.000000	1.098612	0.693147
## Cell_98	0.693147	0.693147	0.693147	0.693147	1.386294	0.000000	0.693147	1.098612	0.000000	0.693147	...	0.000000	0.693147	0.000000	0.000000	0.693147	0.693147	0.693147	0.000000	0.693147	1.098612
## Cell_99	0.000000	0.693147	0.693147	1.098612	1.791759	0.693147	0.693147	0.693147	0.693147	0.693147	...	0.693147	0.693147	0.000000	1.609438	1.098612	0.000000	0.693147	0.000000	0.693147	0.000000
100 rows × 2000 columns

我们看到 .obs_names/.var_names 被用于创建这个 Pandas 对象。

7 将结果写入磁盘

AnnData 带有自己的基于 HDF5 的持久文件格式:h5ad。如果类别数量较少的字符串列还不是类别,AnnData 将自动转换为类别。

adata.write('my_results.h5ad', compression="gzip")
!h5ls 'my_results.h5ad'
## X                        Group
## layers                   Group
## obs                      Group
## obsm                     Group
## obsp                     Group
## uns                      Group
## var                      Group
## varm                     Group
## varp                     Group

8 总结介绍

AnnData 已成为 Python 中单细胞分析的标准,这是有充分理由的——它使用简单,并通过其基于 key 的存储促进了更多可重复的分析。转换为流行的基于 R 的格式进行单细胞分析甚至变得更加容易。

9 视图和副本

为了好玩,让我们看看另一个 metadata 的使用场景。想象一下,观察结果来自一项多年研究中表征 10 个读数的仪器,该研究的样本取自不同地点的不同受试者。我们通常会以某种格式获取该信息,然后将其存储在 DataFrame 中:

obs_meta = pd.DataFrame({
        'time_yr': np.random.choice([0, 2, 4, 8], adata.n_obs),
        'subject_id': np.random.choice(['subject 1', 'subject 2', 'subject 4', 'subject 8'], adata.n_obs),
        'instrument_type': np.random.choice(['type a', 'type b'], adata.n_obs),
        'site': np.random.choice(['site x', 'site y'], adata.n_obs),
    },
    index=adata.obs.index,    # these are the same IDs of observations as above!
)

obs_meta
## time_yr	subject_id	instrument_type	site
## Cell_0	0	subject 4	type b	site x
## Cell_1	4	subject 8	type a	site y
## Cell_2	8	subject 2	type a	site y
## Cell_3	0	subject 2	type a	site y
## Cell_4	4	subject 8	type a	site y
## ...	...	...	...	...
## Cell_95	8	subject 2	type a	site y
## Cell_96	0	subject 1	type b	site y
## Cell_97	2	subject 1	type a	site y
## Cell_98	2	subject 1	type a	site x
## Cell_99	4	subject 1	type a	site y
## 100 rows × 4 columns

这就是我们将读出数据与 metadata 连接起来的方式。当然,下面对 X 的调用的第一个参数也可以只是一个 DataFrame。

adata = ad.AnnData(adata.X, obs=obs_meta, var=adata.var)

现在我们再次拥有一个单个数据容器,可以跟踪所有内容。

print(adata)
## AnnData object with n_obs × n_vars = 100 × 2000
##    obs: 'time_yr', 'subject_id', 'instrument_type', 'site'

对联合数据矩阵进行子集化对于关注变量或观察的子集或为机器学习模型定义训练测试拆分可能很重要。

访问两个 variables 的前 5 行:

adata[:5, ['Gene_1', 'Gene_3']]
## View of AnnData object with n_obs × n_vars = 5 × 2
##    obs: 'time_yr', 'subject_id', 'instrument_type', 'site'

这是一个视图!如果我们想要一个将数据保存在内存中的 AnnData,让我们调用 .copy()

adata_subset = adata[:5, ['Gene_1', 'Gene_3']].copy()

对于视图,我们还可以设置列的前 3 个元素。

print(adata[:3, 'Gene_1'].X.toarray().tolist())
adata[:3, 'Gene_1'].X = [0, 0, 0]
print(adata[:3, 'Gene_1'].X.toarray().tolist())
## [[0.0], [2.0], [0.0]]
## [[0.0], [0.0], [0.0]]

如果您尝试访问 AnnData 视图的部分内容,则会自动复制内容并生成数据存储对象。

adata_subset = adata[:3, ['Gene_1', 'Gene_2']]
adata_subset
## View of AnnData object with n_obs × n_vars = 3 × 2
##    obs: 'time_yr', 'subject_id', 'instrument_type', 'site'

adata_subset.obs['foo'] = range(3)

现在 adata_subset 存储实际数据,不再只是对 adata 的引用。

adata_subset
## AnnData object with n_obs × n_vars = 3 × 2
##    obs: 'time_yr', 'subject_id', 'instrument_type', 'site', 'foo'

显然,您可以使用所有 pandas 对序列或布尔索引进行切片。

adata[adata.obs.time_yr.isin([2, 4])].obs.head()
##    time_yr	subject_id	instrument_type	site
## Cell_1	4	subject 8	type a	site y
## Cell_4	4	subject 8	type a	site y
## Cell_5	2	subject 4	type b	site y
## Cell_7	4	subject 8	type b	site x
## Cell_9	4	subject 1	type a	site y

10 部分读取大数据

如果单个 .h5ad 非常大,您可以使用 backed 模式将其部分读入内存:

adata = ad.read('my_results.h5ad', backed='r')
adata.isbacked
## True

如果你这样做,你需要记住 AnnData 对象与用于读取的文件有一个打开的连接:

adata.filename
## PosixPath('my_results.h5ad')

由于我们在只读模式下使用它,我们不能损坏任何东西。要继续本教程,我们仍然需要显式关闭它:

adata.file.close()

像往常一样,您应该使用 with 语句来避免悬空打开的文件(即将推出的功能)。

操作磁盘上的对象是可能的,但对于稀疏数据是实验性的。因此,我们将其排除在本教程之外。


结束

本文由mdnice多平台发布

  • 0
    点赞
  • 4
    收藏
    觉得还不错? 一键收藏
  • 1
    评论
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值