scanpy | Preprocessing and clustering 3k PBMCs
本文记录使用scanpy处理3k PBMCs scRNA-seq数据的流程。
环境配置
创建一个虚拟环境以方便管理相关的库。
conda create --name pysc python=3.9
conda activate pysc
conda install -c anaconda ipykernel
python -m ipykernel install --user --name pysc
pip3 install scanpy
pip3 install pandas
pip3 install loompy
本文使用PBMCs 3k数据可以在该网址下载（http://cf.10xgenomics.com/samples/cell-exp/1.1.0/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz）.
下载后，将文件解压至当前的data目录下.
# !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
# !mkdir output
import numpy as np
import pandas as pd
import scanpy as sc
sc.settings.verbosity = 3 # verbosity: errors (0), warnings (1), info (2), hints (3)
sc.logging.print_header()
sc.settings.set_figure_params(dpi=80, facecolor='white')
e:\miniconda3\envs\pysc\lib\site-packages\umap\distances.py:1063: NumbaDeprecationWarning: �[1mThe 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.�[0m
@numba.jit()
e:\miniconda3\envs\pysc\lib\site-packages\umap\distances.py:1071: NumbaDeprecationWarning: �[1mThe 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.�[0m
@numba.jit()
e:\miniconda3\envs\pysc\lib\site-packages\umap\distances.py:1086: NumbaDeprecationWarning: �[1mThe 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.�[0m
@numba.jit()
e:\miniconda3\envs\pysc\lib\site-packages\tqdm\auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
from .autonotebook import tqdm as notebook_tqdm
scanpy==1.9.3 anndata==0.9.1 umap==0.5.3 numpy==1.24.4 scipy==1.11.1 pandas==2.0.3 scikit-learn==1.3.0 statsmodels==0.14.0 python-igraph==0.10.6 pynndescent==0.5.10
e:\miniconda3\envs\pysc\lib\site-packages\umap\umap_.py:660: NumbaDeprecationWarning: �[1mThe 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.�[0m
@numba.jit()
results_file='output/pbmc3k.h5ad' # the file that will store the analysis results
adata = sc.read_10x_mtx('data/filtered_gene_bc_matrices/hg19/', # the directory with the `.mtx` file
var_names='gene_symbols', # use gene symbols for the variable names (variables-axis index)
cache=True # write a cache file for faster subsequent reading
)
... reading from cache file cache\data-filtered_gene_bc_matrices-hg19-matrix.h5ad
# remove duplicated symbols
adata.var_names_make_unique() # this is unnecessary if using `var_names='gene_ids'` in `sc.read_10x_mtx`
adata
AnnData object with n_obs × n_vars = 2700 × 32738
var: 'gene_ids'
我们读入的数据为AnnData格式。
AnnData是python中处理带注释的数据的一种格式，读入的数据并不直接读入内存当中，而是创建与磁盘数据的链接来进行处理。对于单细胞测序数据而言，一般与细胞相关数据存储于.obs，而与feature相关数据存储于.var中
anndata_schema
Preprocessing
scanpy包提供的api有几个主要的modules，包括：
- preprocessing:
scanpy.pp, such aspp.calculate_qc_metrics,pp.filter_cells,pp.filter_genes, ...; - tools:
scanpy.tl, such astl.pca,tl.tsne,tl.umap, ...; - plots:
scanpy.pl, such aspl.scatter,pl.highest_expr_genes,pl.umap, ...
接下来，我们先使用scanpy包进行数据预处理，展示高表达基因在各个细胞中的counts。
sc.pl.highest_expr_genes(adata, n_top=20)
normalizing counts per cell
finished (0:00:00)
执行简单的过滤操作。保留至少有200个基因表达的细胞，至少有3个细胞表达的基因。
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
然后，再根据基因的counts和线粒体基因表达进行进一步过滤。
首先，计算线粒体基因比例.
adata.var 存的是feature-level相关的信息，adata.obs 存的是cell-level的信息
adata.var['mt'] = adata.var_names.str.startswith("MT-")
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=False, inplace=True)
adata.obs['pct_counts_mt']
AAACATACAACCAC-1 3.017776
AAACATTGAGCTAC-1 3.793596
AAACATTGATCAGC-1 0.889736
AAACCGTGCTTCCG-1 1.743085
AAACCGTGTATGCG-1 1.224490
...
TTTCGAACTCTCAT-1 2.110436
TTTCTACTGAGGCA-1 0.929422
TTTCTACTTCCTCG-1 2.197150
TTTGCATGAGAGGC-1 2.054795
TTTGCATGCCTCAC-1 0.806452
Name: pct_counts_mt, Length: 2700, dtype: float32
sc.pl.violin(adata, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'], jitter=0.4, multi_panel=True)
e:\miniconda3\envs\pysc\lib\site-packages\seaborn\axisgrid.py:118: UserWarning: The figure layout has changed to tight
self._figure.tight_layout(*args, **kwargs)
output_12_1.png
sc.pl.scatter(adata, x='total_counts', y='pct_counts_mt')
sc.pl.scatter(adata, x='total_counts', y='n_genes_by_counts')
output_13_0.png
output_13_1.png
过滤基因数和线粒体基因占比过高的细胞
通过切片的方法过滤
adata = adata[adata.obs.n_genes_by_counts < 2500, :]
adata = adata[adata.obs.pct_counts_mt < 5, :]
adata
View of AnnData object with n_obs × n_vars = 2638 × 13714
obs: 'n_genes', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt'
var: 'gene_ids', 'n_cells', 'mt', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts'
接下来，我们对read counts进行文库大小校正，并进行log转换。
随后，根据normalized expression鉴定highly-variable genes以供后续PCA等分析使用。
# make a copy of raw count
adata.layers['counts'] = adata.X.copy()
# Total-count normalization
sc.pp.normalize_total(adata, target_sum=1e4)
# log transformation
sc.pp.log1p(adata)
# store normalized data
adata.layers['data'] = adata.X.copy()
normalizing counts per cell
finished (0:00:00)
e:\miniconda3\envs\pysc\lib\site-packages\scanpy\preprocessing\_normalization.py:170: UserWarning: Received a view of an AnnData. Making a copy.
view_to_actual(adata)
# Identify highly-variable genes
sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)
sc.pl.highly_variable_genes(adata)
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)
output_18_1.png
鉴定的基因存储在adata.var.highly_variable中，后续可被PCA、clustering和UMAP/tSNE等函数识别，所以不需要再对数据进行过滤。
adata.var.highly_variable
AL627309.1 False
AP006222.2 False
RP11-206L10.2 False
RP11-206L10.9 False
LINC00115 False
...
AC145212.1 False
AL592183.1 False
AL354822.1 False
PNRC2-1 False
SRSF10-1 False
Name: highly_variable, Length: 13714, dtype: bool
这里将adata.raw设置为校正后的数据，以供后续差异分析和可视化使用。相当于是把默认的表达矩阵替换为校正过的。可以通过.raw.to_adata()再替换为原始数据。
adata.raw = adata
print(adata.raw)
Raw AnnData with n_obs × n_vars = 2638 × 13714
var: 'gene_ids', 'n_cells', 'mt', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts', 'highly_variable', 'means', 'dispersions', 'dispersions_norm'
# filtering hvg
adata = adata[:, adata.var.highly_variable]
# regress out effects of total counts per cell and the percentage of mt genes expressed
sc.pp.regress_out(adata, ['total_counts', 'pct_counts_mt'])
# cale each gene to unit variance. clip values exceeding standard deviation 10.
sc.pp.scale(adata, max_value=10)
# store scaled data
adata.layers['scaled'] = adata.X.copy()
regressing out ['total_counts', 'pct_counts_mt']
sparse input is densified and may lead to high memory use
finished (0:00:03)
# PCA
sc.tl.pca(adata, svd_solver='arpack')
# elbow plot
sc.pl.pca_variance_ratio(adata, log=True)
computing PCA
on highly variable genes
with n_comps=50
finished (0:00:00)
output_24_1.png
adata
AnnData object with n_obs × n_vars = 2638 × 1838
obs: 'n_genes', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt'
var: 'gene_ids', 'n_cells', 'mt', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts', 'highly_variable', 'means', 'dispersions', 'dispersions_norm', 'mean', 'std'
uns: 'log1p', 'hvg', 'pca'
obsm: 'X_pca'
varm: 'PCs'
Computing the neighborhood graph
接下来，我们基于PCA计算细胞的邻接图
sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40)
computing neighbors
using 'X_pca' with n_pcs = 40
finished: added to `.uns['neighbors']`
`.obsp['distances']`, distances for each pair of neighbors
`.obsp['connectivities']`, weighted adjacency matrix (0:00:02)
Embedding the neighborhood graph
使用UMAP可视化
sc.tl.umap(adata)
sc.pl.umap(adata, color=['CST3','NKG7', 'PPBP'])
computing UMAP
finished: added
'X_umap', UMAP coordinates (adata.obsm) (0:00:03)
output_29_1.png
pl.umap默认使用adata.raw中的值作图，由于上面我们将其替换为了normalized后的值。如果想用scaled values作图，可以使用参数use_raw=False
sc.pl.umap(adata, color=['CST3','NKG7', 'PPBP'], use_raw=False)
output_31_0.png
Clustering the neighborhood graph
使用Leiden graph-clustering进行聚类
sc.tl.leiden(adata, resolution=1)
sc.pl.umap(adata, color=['leiden'])
running Leiden clustering
finished: found 8 clusters and added
'leiden', the cluster labels (adata.obs, categorical) (0:00:00)
e:\miniconda3\envs\pysc\lib\site-packages\scanpy\plotting\_tools\scatterplots.py:392: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored
cax = scatter(
output_33_2.png
# save the result
adata.write(results_file)
Finding marker genes
Wilcoxon rank-sum test to identify markers
# compare each group with the rest of cell
sc.tl.rank_genes_groups(adata, 'leiden', method='wilcoxon')
sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False) # sharey: Controls if the y-axis of each panels should be shared.
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:01)
output_36_1.png
adata.write(results_file)
# marker gene
marker_genes = ['IL7R', 'CD79A', 'MS4A1', 'CD8A', 'CD8B', 'LYZ', 'CD14',
'LGALS3', 'S100A8', 'GNLY', 'NKG7', 'KLRB1',
'FCGR3A', 'MS4A7', 'FCER1A', 'CST3', 'PPBP']
# top 5 ranked genes per cluster
pd.DataFrame(adata.uns['rank_genes_groups']['names']).head(5)
image.png
result = adata.uns['rank_genes_groups']
groups = result['names'].dtype.names
pd.DataFrame(
{group + '_' + key[:1] : result[key][group] # create new dictionary
for group in groups for key in ['names', 'pvals']}
).head(5)
image.png
# compare between groups
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: 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)
output_40_1.png
# 0 vs. 1
sc.pl.rank_genes_groups_violin(adata, groups='0', n_genes=8)
e:\miniconda3\envs\pysc\lib\site-packages\seaborn\categorical.py:166: FutureWarning: Setting a gradient palette using color= is deprecated and will be removed in version 0.13. Set `palette='dark:black'` for same effect.
warnings.warn(msg, FutureWarning)
output_41_1.png
# Reload the object with the computed differential expression (i.e. DE via a comparison with the rest of the groups):
adata = sc.read(results_file)
sc.pl.rank_genes_groups_violin(adata, groups='0', n_genes=8)
e:\miniconda3\envs\pysc\lib\site-packages\seaborn\categorical.py:166: FutureWarning: Setting a gradient palette using color= is deprecated and will be removed in version 0.13. Set `palette='dark:black'` for same effect.
warnings.warn(msg, FutureWarning)
output_42_1.png
# plot expression across clusters
sc.pl.violin(adata, ['CST3', 'NKG7', 'PPBP'], groupby='leiden')
output_43_0.png
# rename cluster based on cell types
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
e:\miniconda3\envs\pysc\lib\site-packages\scanpy\plotting\_tools\scatterplots.py:392: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap' will be ignored
cax = scatter(
output_44_2.png
# dotplot for marker genes
sc.pl.dotplot(adata, marker_genes, groupby='leiden')
e:\miniconda3\envs\pysc\lib\site-packages\scanpy\plotting\_dotplot.py:749: UserWarning: No data for colormapping provided via 'c'. Parameters 'cmap', 'norm' will be ignored
dot_ax.scatter(x, y, **kwds)
output_45_1.png
# violin plot
sc.pl.stacked_violin(adata, marker_genes, groupby='leiden')
output_46_0.png
adata
AnnData object with n_obs × n_vars = 2638 × 1838
obs: 'n_genes', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt', 'leiden'
var: 'gene_ids', 'n_cells', 'mt', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts', 'highly_variable', 'means', 'dispersions', 'dispersions_norm', 'mean', 'std'
uns: 'hvg', 'leiden', 'leiden_colors', 'log1p', 'neighbors', 'pca', 'rank_genes_groups', 'umap'
obsm: 'X_pca', 'X_umap'
varm: 'PCs'
obsp: 'connectivities', 'distances'
总结
scanpy中scRNA-seq分析流程包括：
# read in data
adata = sc.read_10x_mtx('data/filtered_gene_bc_matrices/hg19/', cache=True)
# normalization
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
# highly variable genes
sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)
# scale
sc.pp.scale(adata, max_value=10)
# PCA
sc.tl.pca(adata, svd_solver='arpack')
# find neighbors
sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40)
# clustering
sc.tl.leiden(adata)
# UMAP
sc.tl.umap(adata)
sc.pl.umap(adata, color=['leiden'])
# find markers
sc.tl.rank_genes_groups(adata, 'leiden', method='wilcoxon')
由于后续需要在Geneformer中使用该数据，这里我将行名替换为ensembl id，并在细胞注释.obs中增加一列cell_type指定细胞类型和organ_major指定组织类型。同时，替换默认表达矩阵为原始counts。随后，以loom格式存储。
# make dataset suitable for Geneformer
# convert row attributes to Ensembl IDs
adata.var['gene_name'] = adata.var_names
adata.var_names = adata.var['gene_ids']
adata.var.rename(columns={"gene_ids": "ensembl_id"}, inplace=True)
adata.var
image.png
# add cell_type and organ_major columns for cell
adata.obs['cell_type'] = adata.obs['leiden']
adata.obs['organ_major'] = 'immune'
adata.obs.rename(columns={"total_counts": "n_counts"}, inplace=True)
adata.obs
image.png
# replace count matrix with raw counts
adata.X = adata.layers['counts']
adata.X.toarray()[10:15,0:3]
array([[0., 1., 0.],
[0., 0., 0.],
[0., 0., 0.],
[0., 0., 0.],
[0., 0., 0.]], dtype=float32)
# save in h5ad
adata.write(results_file)
# save in loom
adata.write_loom("output/pbmc3k.loom")
The loom file will lack these fields:
{'X_umap', 'PCs', 'X_pca'}
Use write_obsm_varm=True to export multi-dimensional annotations
e:\miniconda3\envs\pysc\lib\site-packages\loompy\bus_file.py:68: NumbaDeprecationWarning: �[1mThe 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.�[0m
def twobit_to_dna(twobit: int, size: int) -> str:
e:\miniconda3\envs\pysc\lib\site-packages\loompy\bus_file.py:85: NumbaDeprecationWarning: �[1mThe 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.�[0m
def dna_to_twobit(dna: str) -> int:
e:\miniconda3\envs\pysc\lib\site-packages\loompy\bus_file.py:102: NumbaDeprecationWarning: �[1mThe 'nopython' keyword argument was not supplied to the 'numba.jit' decorator. The implicit default value for this argument is currently False, but it will be changed to True in Numba 0.59.0. See https://numba.readthedocs.io/en/stable/reference/deprecation.html#deprecation-of-object-mode-fall-back-behaviour-when-using-jit for details.�[0m
def twobit_1hamming(twobit: int, size: int) -> List[int]:
Ref:
https://scanpy-tutorials.readthedocs.io/en/latest/pbmc3k.html#Clustering-3K-PBMCs
