scanpy | 基本单细胞数据分析流程

作者: 尘世中一个迷途小书僮 | 来源:发表于2023-07-21 15:38 被阅读0次

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 as pp.calculate_qc_metrics, pp.filter_cells, pp.filter_genes, ...;
tools: scanpy.tl, such as tl.pca, tl.tsne, tl.umap, ...;
plots: scanpy.pl, such as pl.scatter, pl.highest_expr_genes, pl.umap, ...

https://scanpy.readthedocs.io/en/stable/api.html

接下来，我们先使用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

https://anndata.readthedocs.io/en/latest/

scanpy | 基本单细胞数据分析流程

scanpy | Preprocessing and clustering 3k PBMCs

环境配置

Preprocessing

Computing the neighborhood graph

Embedding the neighborhood graph

Clustering the neighborhood graph

Finding marker genes

总结

相关文章

网友评论

延伸阅读

深度阅读

栏目导航

热点阅读