写在前面
学习一个软件最好的方法就是啃它的官方文档。本着自己学习、分享他人的态度,分享官方文档的中文教程。软件可能随时更新,建议配合官方文档一起阅读。
::: 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_names
为 obs
和 var
轴提供索引。
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.obs
和 adata.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
的更多说明:
- “array-like” metadata 可以源自 Pandas DataFrame、scipy sparse matrix 或 numpy dense array。
- 使用 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多平台发布