3. 分析框架和工具介绍
3.1 单细胞分析框架和强大系统
如前所述，获得计数矩阵后，探索性数据分析阶段开始。由于数据的大小和复杂性，需要专门的工具。虽然在单细胞分析的早期，人们习惯于使用自定义脚本来分析数据，但现在已经存在专门用于此目的的框架。三个最受欢迎的选项是基于R的Bioconductor和Seurat生态系统以及基于Python的scverse生态系统。它们不仅在所使用的编程语言方面有所不同，而且在底层数据结构和可用的专门分析工具方面也有所不同。
Bioconductor是一个开发、支持和共享免费开源软件的项目，重点是对包括单细胞在内的许多不同生物测定的数据进行严格且可重复的分析。大量的开发人员和优质的用户体验以及带有用户友好小插图的丰富文档是Bioconductor的最大优势。Seurat是一款备受推崇的R软件包，专为分析单细胞数据而设计。它为分析的所有步骤提供工具，包括多模式和空间数据。然而，对于极大的数据集（超过 50 万个细胞），这两种基于R的选择都会遇到可扩展性问题，这促使基于Python的社区开发scverse生态系统。scverse 是一个致力于生命科学基础工具的组织和生态系统，最初重点关注单细胞。可扩展性、可扩展性以及与现有Python数据和机器学习工具的强大互操作性是scverse生态系统的一些优势。
在以下部分中，将更详细地介绍scverse生态系统，并解释最重要的概念，重点关注最重要的数据结构。
3.2 使用AnnData存储单模态数据
如前所述，在比对和基因注释之后，基因组数据通常被总结为特征矩阵。该矩阵的形式为number_observations x number_variables——其中scRNA-seq观察结果是细胞条形码，变量是注释基因。在分析过程中，该矩阵的观察结果和变量用计算得出的测量值（例如质量控制指标或潜在空间嵌入）和先验知识（例如源供体或替代基因标识符）进行注释。在scverse生态系统中，AnnData用于将数据矩阵与这些注释相关联。为了实现快速且内存高效的转换，AnnData还支持稀疏矩阵和部分读取。
虽然AnnData与R生态系统的数据结构（例如Bioconductor的SummarizedExperiment或Seurat的object）大致相似，但R包使用转置特征矩阵。
AnnData对象的核心是在X中存储稀疏或密集矩阵（scRNA-Seq 中的计数矩阵）。该矩阵的维度为obs_names x var_names，其中obs（观察值）对应于细胞的条形码var（变量）对应于基因标识符。该矩阵X被Pandas DataFramesobs和var包含，它们分别保存细胞和基因的注释。此外，AnnData保存具有相应维度的观测值 (obsm) 或变量 (varm) 的整个计算矩阵。将细胞与细胞或基因与基因相关联的图形结构通常保存在obsp和varp中。任何不适合其他slot的非结构化数据都将保存在uns中。还可以在layers中存储更多的X值。例如，在countslayer中存储原始的、未标准化的计数数据，在未命名的默认层中存储标准化数据。
AnnData主要针对单模态（例如scRNA-Seq）数据而设计。然而，AnnData的扩展允许高效存储和访问多模态数据。
3.2.1 安装
AnnData可以在PyPI和Conda上安装使用：
pip install anndata
conda install -c conda-forge anndata
3.2.2 初始化一个AnnData对象
首先创建一个具有稀疏计数信息的简单AnnData对象，可以代表基因表达计数。首先，我们导入所需的包。
import numpy as np
import pandas as pd
import anndata as ad
from scipy.sparse import csr_matrix
下一步，我们使用随机泊松分布数据初始化AnnData对象。将AnnData对象命名为adata。
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对象有100个观测值和2000个变量。这相当于100个细胞和2000个基因。我们传递的初始数据可以使用adata.X作为稀疏矩阵进行访问。
adata.X
输出结果：
<100x2000 sparse matrix of type '<class 'numpy.float32'>'
with 126413 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')
3.2.3 添加匹配的元数据metadata
3.2.3.1 观察obs或变量var水平
我们的AnnData对象的核心部分已经完整了。下一步，我们在观察和变量级别添加元数据。此类注释存储在AnnData对象的.obs和.var槽中，分别用于细胞和基因注释。
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 T
Cell_2 B
Cell_3 T
Cell_4 Monocyte
... ...
Cell_95 T
Cell_96 T
Cell_97 B
Cell_98 T
Cell_99 Monocyte
100 rows × 1 columns
如果我们现在再次检查AnnData对象，我们会注意到它已经对obs中的cell_type信息进行了更新。
adata
输出结果：
AnnData object with n_obs × n_vars = 100 × 2000
obs: 'cell_type'
3.2.3.2 使用元数据进行子集化
我们还可以使用随机生成的细胞类型对AnnData对象进行子集化。
bdata = adata[adata.obs.cell_type == "B"]
bdata
输出结果：
View of AnnData object with n_obs × n_vars = 26 × 2000
obs: 'cell_type'
3.2.4 观察/变量水平矩阵
我们还可以拥有具有多个维度的元数据，例如数据的UMAP嵌入。对于此类元数据，AnnData具有.obsm/.varm属性。我们使用keys来标识我们插入的不同矩阵。.obsm/.varm的要求是.obsm矩阵的长度必须等于观测值的数量，因为.n_obs和.varm矩阵的长度必须等于.n_vars。它们各自可以独立地具有不同数量的维度。
让我们从一个随机生成的矩阵开始，我们可以将其解释为我们想要存储的数据的UMAP嵌入，以及一些随机的基因级元数据。
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.“数组样”的元数据可以来源于Pandas DataFrame、scipy稀疏矩阵或numpy密集数组。
- 2.使用scanpy时，它们的值（列）不容易绘制，而
.obs中的内容很容易绘制在图上(例如UMAP)。
3.2.5 非结构化元数据
如上所述，AnnData有.uns，它允许任何非结构化元数据。这可以是任何东西，例如包含一些对数据分析有用的一般信息的列表或字典。最好仅将此slot用于无法有效存储在其他slots中的数据。
adata.uns["random"] = [1, 2, 3]
adata.uns
输出结果：
OverloadedDict, wrapping:
OrderedDict([('random', [1, 2, 3])])
With overloaded keys:
['neighbors'].
3.2.6 Layers层
最后，我们可能有不同形式的原始核心数据，可能一种是标准化的，也可能不是标准化的。这些可以存储在AnnData的不同层中。例如，让我们对原始数据进行log转换并将其存储在一个层中。
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'
我们的原始矩阵X没有修改并且仍然可以访问。我们可以通过将原始X与新层进行比较来验证这一点。
(adata.X != adata.layers["log_transformed"]).nnz == 0
输出结果：
False
3.2.7 转换为DataFrames
我们可以从其中一层layer获取Pandas 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 1.386294 1.386294 0.693147 0.000000 0.693147 1.098612 0.000000 1.098612 ... 0.693147 0.693147 0.000000 0.693147 1.098612 0.000000 0.000000 1.098612 0.693147 0.000000
Cell_1 1.098612 1.098612 0.000000 0.000000 0.000000 0.693147 0.000000 0.693147 1.386294 0.000000 ... 1.098612 1.386294 1.098612 0.000000 0.693147 0.693147 1.386294 0.000000 0.000000 0.693147
Cell_2 0.693147 0.693147 1.098612 0.693147 1.386294 1.386294 0.693147 0.693147 0.000000 1.098612 ... 0.000000 0.693147 0.693147 0.000000 1.386294 0.693147 0.000000 0.000000 0.693147 0.693147
Cell_3 0.000000 0.000000 0.000000 0.693147 0.693147 0.000000 1.386294 1.098612 0.693147 0.693147 ... 0.693147 0.693147 0.000000 0.693147 0.693147 0.693147 0.693147 0.000000 1.386294 0.693147
Cell_4 0.000000 0.693147 0.000000 1.609438 1.098612 0.693147 0.000000 0.000000 0.693147 0.000000 ... 1.098612 1.098612 0.000000 1.386294 1.791759 1.098612 1.609438 0.693147 0.000000 0.693147
... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ... ...
Cell_95 0.693147 1.098612 0.693147 1.386294 1.098612 0.693147 0.000000 0.000000 0.000000 0.693147 ... 0.000000 0.693147 0.000000 0.000000 0.000000 0.693147 0.000000 0.693147 1.386294 0.693147
Cell_96 0.693147 0.000000 0.000000 1.098612 0.000000 0.693147 0.693147 1.098612 0.000000 0.000000 ... 1.386294 0.693147 1.098612 0.693147 0.000000 0.693147 1.098612 1.386294 0.693147 0.000000
Cell_97 0.693147 0.000000 0.000000 0.000000 0.000000 0.000000 1.098612 0.000000 0.000000 0.693147 ... 0.000000 0.693147 1.386294 0.000000 0.000000 0.000000 0.693147 0.693147 1.609438 0.693147
Cell_98 0.693147 0.693147 0.000000 1.098612 0.000000 0.693147 0.693147 1.098612 0.000000 0.000000 ... 0.000000 0.000000 0.693147 0.693147 0.000000 0.693147 1.386294 0.693147 1.098612 1.386294
Cell_99 1.098612 1.791759 0.000000 0.000000 0.000000 0.000000 1.098612 1.609438 0.693147 0.693147 ... 0.693147 0.000000 0.000000 0.693147 0.693147 0.693147 0.693147 0.000000 0.000000 1.386294
100 rows × 2000 columns
3.2.8 AnnData对象的读取和写入
AnnData对象可以在磁盘上以分层数组存储（如 HDF5 或 Zarr）。AnnData带有自己的基于HDF5的持久文件格式：h5ad。现在，我们将以h5ad格式保存AnnData对象。
adata.write("my_results.h5ad", compression="gzip")
!h5ls 'my_results.h5ad'
输出结果：
X Group
layers Group
obs Group
obsm Group
uns Group
var Group
varm Group
以及将其读回。
adata_new = ad.read_h5ad("my_results.h5ad")
adata_new
输出结果：
AnnData object with n_obs × n_vars = 100 × 2000
obs: 'cell_type'
uns: 'random'
obsm: 'X_umap'
varm: 'gene_stuff'
layers: 'log_transformed'
3.3 scanpy分析单模态数据
现在我们了解了单模态单细胞分析的基本数据结构，但是我们如何实际分析存储的数据呢？scverse生态系统中存在多种用于分析特定组学数据的工具。例如，scanpy提供用于一般RNA-Seq重点分析的工具，squidpy重点关注空间转录组学，而scirpy提供用于分析T细胞受体 (TCR) 和B细胞受体 (BCR) 数据的工具。尽管存在许多针对各种数据模式的横向扩展，但它们通常使用scanpy的一些预处理和可视化功能。
更具体地说，scanpy是一个建立在AnnData之上的Python包，用于促进单细胞基因表达数据的分析。可通过scanpy访问预处理、嵌入、可视化、聚类、差异基因表达测试、伪时间和轨迹推断以及基因调控网络模拟的多种方法。基于Python数据科学和机器学习库，使scanpy能够扩展到数百万个细胞。
3.3.1 安装
Scanpy可以在PyPI和Conda安装使用：
pip install scanpy
conda install -c conda-forge scanpy
3.3.2 Scanpy的API设计
scanpy的框架设计是将属于同一步骤的函数分组到相应的模块中。例如，所有预处理功能都可以在scanpy.preprocessing模块中使用，所有非预处理的数据矩阵转换都可以在scanpy.tools中使用，所有可视化都可以在scanpy.plot中使用。这些模块通常在导入scanpy之后访问，例如import scanpy as sc以及相应的缩写sc.pp（用于预处理）、sc.tl（用于工具）和sc.pl（用于绘图）。所有读或写数据的模块都是直接访问的。此外，用于各种数据集的模块可用作sc.datasets。所有函数及其相应的参数和潜在的示例图都记录在scanpy API文档中。
3.3.3 Scanpy示范
在下面的内容中，我们将很快演示使用scanpy进行分析的工作流程。
选择的数据集是健康捐献者的2700个外周血单核细胞的数据集，这些细胞在Illumina NextSeq 500上进行了测序。
作为第一步，我们导入scanpy。 我们设置Matplotlib绘图作为所有scanpy绘图的默认值，并且最后打印scanpy的header。包含当前环境中所有相关Python包的版本，包括Scanpy和AnnData。
import scanpy as sc
sc.settings.set_figure_params(dpi=80, facecolor="white")
sc.logging.print_header()
输出结果：
scanpy==1.9.1 anndata==0.7.8 umap==0.5.2 numpy==1.20.3 scipy==1.7.2 pandas==1.4.3 scikit-learn==1.0.1 statsmodels==0.13.1 python-igraph==0.9.10 pynndescent==0.5.5
接下来，我们使用scanpy作为AnnData对象获取3k PBMC数据集。数据集会自动下载并保存，也可以使用scanpy设置指定的路径。默认情况下它是在文件夹data。
adata = sc.datasets.pbmc3k()
adata
输出结果：
100%|██████████| 5.58M/5.58M [00:01<00:00, 3.25MB/s]
AnnData object with n_obs × n_vars = 2700 × 32738
var: 'gene_ids'
返回的AnnData对象有2700个细胞和32738个基因，var中还包含基因ID。
adata.var
输出结果：
gene_ids
index
MIR1302-10 ENSG00000243485
FAM138A ENSG00000237613
OR4F5 ENSG00000186092
RP11-34P13.7 ENSG00000238009
RP11-34P13.8 ENSG00000239945
... ...
AC145205.1 ENSG00000215635
BAGE5 ENSG00000268590
CU459201.1 ENSG00000251180
AC002321.2 ENSG00000215616
AC002321.1 ENSG00000215611
32738 rows × 1 columns
如上所述，scanpy的所有分析函数都可以通过sc.[pp, tl, pl]访问。作为第一步，我们使用scanpy来显示在所有细胞中每个单细胞中产生最高计数分数的基因。我们只需调用sc.pl.highest_expr_genes函数，传递AnnData对象，该对象在几乎所有情况下都是任何scanpy函数的第一个参数，并指定我们希望显示前20个表达基因。
sc.pl.highest_expr_genes(adata, n_top=20)
显然，MALAT1是表达最多的基因，经常在Poly-A捕获的scRNA-Seq数据中检测到，与实验方案无关。该基因已被证明与细胞健康呈负相关。特别是死亡/垂死的细胞具有更高的MALAT1表达。
现在，我们使用scanpy的预处理模块过滤检测到的基因少于200个的细胞以及在少于3个细胞中发现的基因，以获得粗略的质量阈值。
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)
单细胞RNA-Seq分析的一个常见步骤是通过主成分分析 (PCA) 等进行降维，以揭示主要的差异变量。这也会对数据进行去噪。scanpy的preprocessing或者tools函数提供PCA方法。
sc.tl.pca(adata, svd_solver="arpack")
相应的绘图函数允许我们将基因传递给颜色参数。相应的值会自动从AnnData对象中提取。
sc.pl.pca(adata, color="CST3")
任何高级嵌入和下游计算的基本步骤是使用数据矩阵的PCA表示来计算邻域图。它会自动用于需要它的其他工具，例如UMAP的计算。
sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40)
我们现在使用计算邻域图将单元格嵌入到UMAP中，UMAP是在scanpy中实现的许多高级降维算法之一。
sc.tl.umap(adata)
sc.pl.umap(adata, color=["CST3", "NKG7", "PPBP"])
scanpy的文档提供了教程，需要的可以认真阅读。
3.4 后记
其实除了处理单模态数据外，AnnData还可以存储和操作多模态数据，如CITE-seq等同时测量RNA和表面蛋白数据，但这一部分是比较进阶的内容，目前先只介绍简单的普通单细胞RNA-seq分析的内容，多模态的数据，之后再详细分享给大家。
