(未测试)单细胞数据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
请关注“恒诺新知”微信公众号,感谢“R语言“,”数据那些事儿“,”老俊俊的生信笔记“,”冷🈚️思“,“珞珈R”,“生信星球”的支持!