（未测试）单细胞数据scanpy 学习笔记

Scanpy 是一个基于 Python 分析单细胞数据的软件包，内容包括预处理，可视化，聚类，拟时序分析和差异表达分析等。本文翻译自 scanpy 的官方教程 Preprocessing and clustering 3k PBMCs^[1]，用 scanpy 重现 Seurat 聚类教程^[2] 中的绝大部分内容。

0. scanpy 安装

直接用 pip ，安装成功。

pip install scanpy[louvain]

Docker

docker pull fastgenomics/scanpy:1.4-p368-v1-stretch-slim

1. 载入数据

# 下载PBMC 数据集
## 其实就是 Seurat 那个示例数据，之前下过就不用重复下了
!mkdir data
!wget http://cf.10xgenomics.com/samples/cell-exp/1.1.0/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz -O data/pbmc3k_filtered_gene_bc_matrices.tar.gz
!cd data; tar -xzf pbmc3k_filtered_gene_bc_matrices.tar.gz

import numpy as np
import pandas as pd
import scanpy as sc

/Users/baozhiwei/anaconda3/lib/python3.7/site-packages/anndata/core/anndata.py:17: FutureWarning: pandas.core.index is deprecated and will be removed in a future version.  The public classes are available in the top-level namespace.
  from pandas.core.index import RangeIndex

# verbosity 的取值表示测试结果显示的详细程度，数字越大越详细
## errors (0), warnings (1), info (2), hints (3)
sc.settings.verbosity = 3
# 输出版本号
sc.logging.print_versions()
# set_figure_params 设置图片的分辨率/大小以及其他样式
sc.settings.set_figure_params(dpi=80)

scanpy==1.4.5.post3 anndata==0.6.22.post1 umap==0.3.10 numpy==1.18.1 scipy==1.4.1 pandas==1.0.1 scikit-learn==0.22.1 statsmodels==0.11.0 python-igraph==0.7.1

# 设置结果文件保存路径
results_file = './pbmc3k.h5ad'

# 导入 10X 数据
adata = sc.read_10x_mtx(
    './data/filtered_gene_bc_matrices/hg19/',  # 包含有 `.mtx` 文件的目录
    var_names='gene_symbols',                  # 用 gene symbols 作为变量名 (variables-axis index)
    cache=True)                                # 使用缓存文件加快读取

... writing an h5ad cache file to speedup reading next time

scanpy 的结果文件是 AnnData 对象，意为 annotated data。AnnData 用了一些广义的单词来描述「细胞」和「基因」：将细胞称为观察值 observations ，将基因称为变量 variables 。AnnData 包括四个可以存储信息的区域：

•adata.X 存储 count matrix，数据类型为稀疏矩阵 scipy.sparse.csr.csr_matrix

•adata.obs 存储关于 obervations(cells) 的 metadata，数据类型为 dataframe

•adata.var 存储关于 variables(genes) 的 metadata，数据类型为 dataframe

•AnnData.uns 存储后续附加的其他非结构化信息

•adata.obs_names 和 adata.var_names index

细胞名和基因名可分别通过 adata.obs_names 和 adata.var_names 查看。 AnnData 对象可以像 dataframe 一样进行切片操作，例如，adata_subset = adata[:, list_of_gene_names]。

# 索引去重，若上一步中使用 `var_names='gene_ids'` 则这一步非必须进行
# 其实也可以直接用 pandas 判断索引是否重复
# adata.var_names.is_unique
adata.var_names_make_unique()

adata

AnnData object with n_obs × n_vars = 2700 × 32738 
    var: 'gene_ids'

2. 数据预处理

可视化所有细胞中计数最多的基因。

sc.pl.highest_expr_genes(adata, n_top=20, )

normalizing counts per cell
    finished (0:00:00)

数据初筛

# 保留至少在三个细胞中表达的基因，保留至少包含 200 个基因的细胞
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)

filtered out 19024 genes that are detected in less than 3 cells

标注线粒体基因并计算每个细胞中的线粒体基因比例。

mito_genes = adata.var_names.str.startswith('MT-')
# 计算每个细胞中的线粒体基因比例
# 使用`.A1` 将 numpy.matrix 转为一维数组 ndarray
adata.obs['percent_mito'] = np.sum(
    adata[:, mito_genes].X, axis=1).A1 / np.sum(adata.X, axis=1).A1
# 把每个细胞的基因计数添加到 adata 中
adata.obs['n_counts'] = adata.X.sum(axis=1).A1

方法二

还可以使用 pp.calculate_qc_metrics 来计算每个细胞中线粒体基因的数量，以及其他一系列 QC 数据。关于更多 preprocessing 的内容可参见文档^[3]。

# 标记线粒体基因
adata.var['mito'] = mito_genes

calculate_qc_metrics 函数会返回由两个 dataframe 组成的元组，一个是细胞的 QC 矩阵，另一个是基因的 QC 矩阵。

qc = sc.pp.calculate_qc_metrics(adata,qc_vars=['mito'])
cell_qc_dataframe = qc[0]
gene_qc_dataframe = qc[1]

cell_qc_dataframe.iloc[:4, :4]

gene_qc_dataframe.iloc[:4, :4]

用小提琴图可视化

sc.pl.violin(adata, ['n_genes', 'n_counts', 'percent_mito'],
             jitter=0.4, multi_panel=True)

可视化特征之间的关系。

sc.pl.scatter(adata, x='n_counts', y='percent_mito')
sc.pl.scatter(adata, x='n_counts', y='n_genes')

adata

AnnData object with n_obs × n_vars = 2700 × 13714 
    obs: 'n_genes', 'percent_mito', 'n_counts'
    var: 'gene_ids', 'n_cells'

过滤线粒体基因比例 > 5% 和基因总数 >2500 的细胞。

adata = adata[adata.obs.n_genes < 2500, :]
adata = adata[adata.obs.percent_mito < 0.05, :]

adata

AnnData object with n_obs × n_vars = 2638 × 13714 
    obs: 'n_genes', 'percent_mito', 'n_counts'
    var: 'gene_ids', 'n_cells'

标准化数据

sc.pp.normalize_total(adata, target_sum=1e4)

normalizing counts per cell
    finished (0:00:01)

进行自然对数转换

sc.pp.log1p(adata)

可以将 AnnData 对象的 .raw 属性设置为经归一化和对数化的原始基因表达值，供之后的可视化分析使用。

adata.raw = adata

选择差异基因

sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)

extracting highly variable genes
    finished (0:00:00)
--> added
    'highly_variable', boolean vector (adata.var)
    'means', float vector (adata.var)
    'dispersions', float vector (adata.var)
    'dispersions_norm', float vector (adata.var)

这一步在 adata.var 中添加了四列内容（highly_variable，means，dispersions，dispersions_norm）。

sc.pl.highly_variable_genes(adata)

取出高度差异的基因。

adata = adata[:, adata.var.highly_variable]

校正细胞基因计数和线粒体基因比例的影响。

sc.pp.regress_out(adata, ['n_counts', 'percent_mito'])

regressing out ['n_counts', 'percent_mito']
    sparse input is densified and may lead to high memory use
    finished (0:00:04)

数据缩放

sc.pp.scale(adata, max_value=10)

3. PCA

sc.tl.pca(adata, svd_solver='arpack')
# svd_solver 指定奇异值分解 SVD 的方法

computing PCA with n_comps = 50
    on highly variable genes
    finished (0:00:00)

绘制碎石图，确定数据维度。

sc.pl.pca_variance_ratio(adata, log=True)

保存结果。

adata.write(results_file)

adata

AnnData object with n_obs × n_vars = 2638 × 1838
    obs: 'n_genes', 'percent_mito', 'n_counts'
    var: 'gene_ids', 'n_cells', 'highly_variable', 'means', 'dispersions', 'dispersions_norm'
    uns: 'log1p', 'pca'
    obsm: 'X_pca'
    varm: 'PCs'

4. 细胞聚类

为了重现 Seurat 的结果，我们使用下面的参数。

sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40)

computing neighbors
    using 'X_pca' with n_pcs = 40
    finished: added to `.uns['neighbors']`
    'distances', distances for each pair of neighbors
    'connectivities', weighted adjacency matrix (0:00:01)

sc.tl.leiden(adata)

running Leiden clustering
    finished: found 8 clusters and added
    'leiden', the cluster labels (adata.obs, categorical) (0:00:00)

非线性降维 UMAP

sc.tl.umap(adata)

computing UMAP
    finished: added
    'X_umap', UMAP coordinates (adata.obsm) (0:00:04)

绘制聚类图

sc.pl.umap(adata, color=['leiden', 'CST3', 'NKG7'])

保存结果。

adata.write(results_file)

5. 找 marker gene

接下来计算每个 cluster 中高度差异基因的排名。默认，若之前已初始化 AnnData 的 .raw 属性，则使用会该属性。

t 检验

sc.tl.rank_genes_groups(adata, 'leiden', method='t-test')
sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)

ranking genes
    finished: added to `.uns['rank_genes_groups']`
    'names', sorted np.recarray to be indexed by group ids
    'scores', sorted np.recarray to be indexed by group ids
    'logfoldchanges', sorted np.recarray to be indexed by group ids
    'pvals', sorted np.recarray to be indexed by group ids
    'pvals_adj', sorted np.recarray to be indexed by group ids (0:00:00)

sc.settings.verbosity = 2  # reduce the verbosity

Wilcoxon 秩和检验

与 t 检验的结果非常相似。我们建议在文章中使用后者，参见 Sonison＆Robinson（2018）^[4]。

sc.tl.rank_genes_groups(adata, 'leiden', method='wilcoxon')
sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)

ranking genes
    finished (0:00:02)

保存结果

adata.write(results_file)

逻辑回归

或者，也可以根据 Natranos et al. (2018)^[5] 的建议使用逻辑回归对基因进行排名。本质区别在于，这里使用多元变量，而传统的差异检验是单变量。更多细节可参考 Clark et al. (2014)^[6] 。

sc.tl.rank_genes_groups(adata, 'leiden', method='logreg')
sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)

ranking genes
    finished (0:00:03)

除 IL7R 仅在 t 检验的结果中发现，以及仅在其他两种检验方法中发现的 FCER1A以外，其他标记基因均可通过所有检验方法得到。

Louvain Group	Markers	Cell Type
0	IL7R	CD4 T cells
1	CD14, LYZ	CD14+ Monocytes
2	MS4A1	B cells
3	CD8A	CD8 T cells
4	GNLY, NKG7	NK cells
5	FCGR3A, MS4A7	FCGR3A+ Monocytes
6	FCER1A, CST3	Dendritic Cells
7	PPBP	Megakaryocytes

定义一个标记基因的列表，方便之后的分析。

marker_genes = ['IL7R', 'CD79A', 'MS4A1', 'CD8A', 'CD8B', 'LYZ', 'CD14',
                'LGALS3', 'S100A8', 'GNLY', 'NKG7', 'KLRB1',
                'FCGR3A', 'MS4A7', 'FCER1A', 'CST3', 'PPBP']

重新加载保存有 Wilcoxon 秩和检验结果的对象。

显示每个 cluster 中排名前五的基因。

pd.DataFrame(adata.uns['rank_genes_groups']['names']).head(5)

显示每个基因在每一 cluster 中的对应 p 值。

result = adata.uns['rank_genes_groups']
groups = result['names'].dtype.names
pd.DataFrame(
    {group + '_' + key[:1]: result[key][group]
    for group in groups for key in ['names', 'pvals']}).head(5)

cluster 与 cluster 之间进行比较。

sc.tl.rank_genes_groups(adata, 'leiden', groups=['0'], reference='1', method='wilcoxon')
sc.pl.rank_genes_groups(adata, groups=['0'], n_genes=20)

ranking genes
    finished (0:00:01)

还可用 sc.pl.rank_genes_groups_violin 进一步获得信息更丰富的图。

sc.pl.rank_genes_groups_violin(adata, groups='0', n_genes=8)

重新载入对象，将单一 cluster 与其他剩余的 cluster 进行比较。

adata = sc.read(results_file)

sc.pl.rank_genes_groups_violin(adata, groups='0', n_genes=8)

如果要比较某个基因在不同组之间的差异，可使用以下代码。

sc.pl.violin(adata, ['CST3', 'NKG7', 'PPBP'], groupby='leiden')

6. 标记细胞类型

new_cluster_names = [
    'CD4 T', 'CD14 Monocytes',
    'B', 'CD8 T',
    'NK', 'FCGR3A Monocytes',
    'Dendritic', 'Megakaryocytes']
adata.rename_categories('leiden', new_cluster_names)

sc.pl.umap(adata, color='leiden', legend_loc='on data', title='', frameon=False, save='.pdf')

WARNING: saving figure to file figures/umap.pdf

可视化 marker genes。

ax = sc.pl.dotplot(adata, marker_genes, groupby='leiden')

绘制不同细胞亚群中多个基因的小提琴图。

ax = sc.pl.stacked_violin(adata, marker_genes, groupby='leiden', rotation=90)

adata

AnnData object with n_obs × n_vars = 2638 × 1838
    obs: 'n_genes', 'percent_mito', 'n_counts', 'leiden'
    var: 'gene_ids', 'n_cells', 'highly_variable', 'means', 'dispersions', 'dispersions_norm'
    uns: 'leiden', 'leiden_colors', 'neighbors', 'pca', 'rank_genes_groups', 'umap'
    obsm: 'X_pca', 'X_umap'
    varm: 'PCs'

保存数据。

adata.write(results_file, compression='gzip')

引用链接

[1]  https://scanpy-tutorials.readthedocs.io/en/latest/pbmc3k.html
[2]  http://satijalab.org/seurat/pbmc3k_tutorial.html
[3]  https://icb-scanpy.readthedocs-hosted.com/en/stable/api/index.html#module-scanpy.pp
[4]  https://doi.org/10.1038/nmeth.4612
[5]  https://doi.org/10.1101/258566
[6]  https://doi.org/10.1186/1471-2105-15-79