单细胞分析的 Python 包 Scanpy（图文详解）

一、安装

如果没有conda 基础，参考： Conda 安装使用图文详解（2021版）

pip install scanpy
conda install -y -c conda-forge leidenalg

二、使用

1、准备工作

# 载入包
import numpy as np
import pandas as pd
import scanpy as sc

# 设置
sc.settings.verbosity = 3             # 设置日志等级: errors (0), warnings (1), info (2), hints (3)
sc.logging.print_header()
sc.settings.set_figure_params(dpi=80, facecolor='white')

# 用于存储分析结果文件的路径
results_file = 'write/pbmc3k.h5ad'

# 载入文件
adata = sc.read_10x_mtx(
    './filtered_gene_bc_matrices/hg19/',  # mtx 文件目录
    var_names='gene_symbols',             # 使用 gene_symbols 作为变量名
    cache=True)                           # 载入包
import numpy as np
import pandas as pd
import scanpy as sc

# 设置
sc.settings.verbosity = 3             # 设置日志等级: errors (0), warnings (1), info (2), hints (3)
sc.logging.print_header()
sc.settings.set_figure_params(dpi=80, facecolor='white')

# 用于存储分析结果文件的路径
results_file = 'write/pbmc3k.h5ad'

# 载入文件
adata = sc.read_10x_mtx(
    './filtered_gene_bc_matrices/hg19/',  # mtx 文件目录
    var_names='gene_symbols',             # 使用 gene_symbols 作为变量名
    cache=True)                           # 写入缓存，可以更快的读取文件

2、预处理

显示在所有细胞中在每个单细胞中产生最高计数分数的基因

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

过滤低质量细胞样本

过滤在少于三个细胞中表达，或一个细胞中表达少于200个基因的细胞样本

sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)

过滤包含线粒体基因和表达基因过多的细胞

线粒体基因的转录本比单个转录物分子大，并且不太可能通过细胞膜逃逸。因此，检测出高比例的线粒体基因，表明细胞质量差（Islam et al. 2014; Ilicic et al. 2016）。

表达基因过多可能是由于一个油滴包裹多个细胞，从而检测出比正常检测要多的基因数，因此要过滤这些细胞。

adata.var['mt'] = adata.var_names.str.startswith('MT-')  # 将线粒体基因标记为 mt
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=False, inplace=True)
sc.pl.violin(adata, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'],
             jitter=0.4, multi_panel=adata.var['mt'] = adata.var_names.str.startswith('MT-')  # 将线粒体基因标记为 mt
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=False, inplace=True)
sc.pl.violin(adata, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'],
             jitter=0.4, multi_panel=True)

生成的三张小提琴图代表：表达基因的数量，每个细胞包含的表达量，线粒体基因表达量的百分比。

过滤

sc.pl.scatter(adata, x='total_counts', y='pct_counts_mt')
sc.pl.scatter(adata, x='total_counts', y=sc.pl.scatter(adata, x='total_counts', y='pct_counts_mt')
sc.pl.scatter(adata, x='total_counts', y='n_genes_by_counts')

过滤线粒体基因表达过多或总数过多的细胞，也就是红框标识的样本。

# 获取线粒体基因占比在 5% 以下的细胞样本
adata = adata[adata.obs.pct_counts_mt < 5, :]
# 获取表达基因数在 2500 以下的细胞样本
adata = adata[adata.obs.n_genes_by_counts < # 获取线粒体基因占比在 5% 以下的细胞样本
adata = adata[adata.obs.pct_counts_mt < 5, :]
# 获取表达基因数在 2500 以下的细胞样本
adata = adata[adata.obs.n_genes_by_counts < 2500, :]

3、检测特异性基因

归一化

# 归一化，使得不同细胞样本间可比
sc.pp.normalize_total(adata, target_sum=# 归一化，使得不同细胞样本间可比
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)

存储数据

将 AnnData 对象的 .raw 属性设置为归一化和对数化的原始基因表达，以便以后用于基因表达的差异测试和可视化。这只是冻结了 AnnData 对象的状态。

adata.raw = adata

识别特异性基因

# 计算
sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)
# 绘制特异性基因散点图
sc.pl.highly_variable_genes(adata)

获取只有特异性基因的数据集

# 获取只有特异性基因的数据集
adata = adata[:, adata.var.highly_variable]
# 回归每个细胞的总计数和表达的线粒体基因的百分比的影响。
sc.pp.regress_out(adata, ['total_counts', 'pct_counts_mt'])
# 将每个基因缩放到单位方差。阈值超过标准偏差 10。
sc.pp.scale(adata, max_value=# 获取只有特异性基因的数据集
adata = adata[:, adata.var.highly_variable]
# 回归每个细胞的总计数和表达的线粒体基因的百分比的影响。
sc.pp.regress_out(adata, ['total_counts', 'pct_counts_mt'])
# 将每个基因缩放到单位方差。阈值超过标准偏差 10。
sc.pp.scale(adata, max_value=10)

4、主成分分析（Principal component analysis）

通过运行主成分分析 (PCA) 来降低数据的维数，可以对数据进行去噪并揭示不同分群的主因素。

# 绘制 PCA 图
sc.tl.pca(adata, svd_solver='arpack')
sc.pl.pca(adata, color=# 绘制 PCA 图
sc.tl.pca(adata, svd_solver='arpack')
sc.pl.pca(adata, color='CST3')

检查单个 PC 对数据总方差的贡献，这可以提供给我们应该考虑多少个 PC 以计算细胞的邻域关系的信息，例如用于后续的聚类函数 sc.tl.louvain() 或 tSNE 聚类 sc.tl.tsne()。

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

5、领域图，聚类图（Neighborhood graph）

使用数据矩阵的 PCA 表示来计算细胞的邻域图。为了重现 Seurat 的结果，我们采用以下值。

建议使用 UMAP ，它可能比 tSNE 更忠实于流形的全局连通性，因此能更好地保留轨迹。

sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40)
sc.tl.umap(adata)
# 如果设置了 adata 的 .raw 属性时，下图显示了“raw”（标准化、对数化但未校正）基因表达矩阵。
sc.pl.umap(adata, color=['CST3', 'NKG7', sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40)
sc.tl.umap(adata)
# 如果设置了 adata 的 .raw 属性时，下图显示了“raw”（标准化、对数化但未校正）基因表达矩阵。
sc.pl.umap(adata, color=['CST3', 'NKG7', 'PPBP'])

为了绘制缩放矫正的基因表达聚类图，需要使用 use_raw=False 参数。

sc.pl.umap(adata, color=['CST3', 'NKG7', 'PPBP'], use_raw=sc.pl.umap(adata, color=['CST3', 'NKG7', 'PPBP'], use_raw=False)

目前还没有计算出各个细胞类群，下面进行聚类

Leiden 图聚类

# 计算
sc.tl.leiden(adata)
# 绘制
sc.pl.umap(adata, color=[# 计算
sc.tl.leiden(adata)
# 绘制
sc.pl.umap(adata, color=['leiden'])

6、检索标记基因

先计算每个 leiden 分群中高度差异基因的排名，取排名前 25 的基因。

默认情况下，使用 AnnData 的 .raw 属性。

T-test

最简单和最快的方法是 t 检验。

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

Wilcoxon rank-sum

Wilcoxon rank-sum (Mann-Whitney-U) 检验的结果非常相似，还可以使用其他的差异分析包，如 MAST、limma、DESeq2 和 diffxpy。

sc.tl.rank_genes_groups(adata, 'leiden', method='wilcoxon')
sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)
sc.tl.rank_genes_groups(adata, 'leiden', method='wilcoxon')
sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)
# 保存这次的数据结果
adata.write(results_file)

逻辑回归

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

使用逻辑回归对基因进行排名 Natranos et al. (2018)，这里使用多变量方法，而传统的差异测试是单变量 Clark et al. (2014)

除了仅由 t 检验发现的 IL7R 和由其他两种方法发现的 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', marker_genes = ['IL7R', 'CD79A', 'MS4A1', 'CD8A', 'CD8B', 'LYZ', 'CD14',
                'LGALS3', 'S100A8', 'GNLY', 'NKG7', 'KLRB1',
                'FCGR3A', 'MS4A7', 'FCER1A', 'CST3', 'PPBP']

载入数据

# 使用 Wilcoxon Rank-Sum 测试结果重新加载已保存的对象
adata = sc.read(results_file)

获取聚类分组和分数

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(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)

Group 1 与 Group 2 比较，进行差异分析

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

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

Group 0 与其余组的比较进行差异分析

adata = sc.read(results_file)
sc.pl.rank_genes_groups_violin(adata, groups='0', n_genes=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=sc.pl.violin(adata, ['CST3', 'NKG7', 'PPBP'], groupby='leiden')

根据已知的细胞标记，注释细胞类型

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=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')

可视化每个类群的标记基因

气泡图显示：

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

小提琴图显示

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

7、保存数据

保存压缩文件

如果只想将其用于可视化的人共享此文件，减少文件大小的一种简单方法是删除缩放和校正的数据矩阵。

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

保存为 h5ad 数据

adata.raw.to_adata().write(adata.raw.to_adata().write('./write/pbmc3k_withoutX.h5ad')

读取使用 adata = sc.read_h5ad(‘./write/pbmc3k_withoutX.h5ad’)

导出数据子集

# 导出聚类数据
adata.obs[['n_counts', 'louvain_groups']].to_csv('./write/pbmc3k_corrected_louvain_groups.csv')
# 导出PCA数据
adata.obsm.to_df()[['X_pca1', 'X_pca2']].to_csv(# 导出聚类数据
adata.obs[['n_counts', 'louvain_groups']].to_csv('./write/pbmc3k_corrected_louvain_groups.csv')
# 导出PCA数据
adata.obsm.to_df()[['X_pca1', 'X_pca2']].to_csv('./write/pbmc3k_corrected_X_pca.csv')

8、番外

我之前在处理较多数据量的时候，会有些地方不一样，具体每个数据集的处理也会有比较大的自由度，比如：

在检测线粒体基因时，这里在质控时，已经把线粒体基因直接剔除。

在做 UMAP 时，可以看到一些类群间的联系和轨迹。

做 TSNE时，可以看到类群间比较干净利索，整体比较“饱满”。

其他 Scanpy 的使用教程：
scanpy 单细胞分析包图文详解 01 | 深入理解 AnnData 数据结构