Scanpy 是一个基于 Python 分析单细胞数据的软件包，内容包括预处理，可视化，聚类，拟时序分析和差异表达分析等。
https://genomebiology.biomedcentral.com/articles/10.1186/s13059-017-1382-0
https://scanpy.readthedocs.io/en/stable/tutorials.html#integrating-datasets
有人可能会说：单细胞分析使用Seurat，monocle等R包会更加方便。但是实际分析中，当单细胞数据过多时，Seurat和monocle会产生内存不足的问题，尤其是monocle2，单细胞数据大于50k的时候，往往会产生内存不足，即使是在服务器上运行。而基于python的单细胞转录分析包scanpy，能很好的解决内存不足的问题，根据官网说明，scanpy可以解决超过1,000,000个单细胞的数据集，这是基于R分析包不敢想象的。(Seurat:https://satijalab.org/seurat/articles/pbmc3k_tutorial.html)

单细胞转录数据分析之Scanpy：https://www.jianshu.com/p/e22a947e6c60
单细胞转录组之Scanpy - 轨迹推断/拟时序分析：https://www.jianshu.com/p/0b2ca0e0b544
单细胞转录组之Scanpy - 样本整合分析：https://www.jianshu.com/p/beef8a8be360
单细胞空间转录分析之Scanpy：
单细胞空间转录分析之Scanpy-结合单细胞转录组：

Fig1. 网址主页

这儿将系统全面演示了如何使用Scanpy来处理单细胞转录组学数据。

scanpy相关python 包安装（安装好python3之后，终端运行）：

pip install scanpy==1.6.1 -i https://pypi.doubanio.com/simple/

数据下载，这儿我们使用了Seurat官网使用的pbmc3k数据进行测试（我是很早之前就下载过了）：

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
tar -xzf pbmc3k_filtered_gene_bc_matrices.tar.gz

运行python3，导入相关包：

import scanpy as sc
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
sc.settings.verbosity = 3 # verbosity 的取值表示测试结果显示的详细程度，数字越大越详细
sc.logging.print_versions() # 输出版本号
sc.settings.set_figure_params(dpi=80) # set_figure_params 设置图片的分辨率/大小以及其他样式
import os #在服务器运行，习惯性会设置一个输出路径，用于保存pdf图片
os.getcwd()  #查看当前路径
os.chdir('./filtered_gene_bc_matrices/scanpy') #修改路径
os.getcwd()
results_file = 'pbmc3k.h5ad' ##置结果文件保存路径

读取并查看数据：

# 导入 10X 数据
data=sc.read_10x_mtx('./filtered_gene_bc_matrices/hg19/',var_names='gene_symbols',   cache=True)
data.var_names_make_unique()  # 索引去重，若上一步中使用 `var_names='gene_ids'` 则这一步非必须进行
#data.X 存储 count matrix，数据类型为稀疏矩阵 scipy.sparse.csr.csr_matrix
#data.obs 存储关于 obervations(cells) 的 metadata，数据类型为 dataframe
#data.var 存储关于 variables(genes) 的 metadata，数据类型为 dataframe
#AnnData.uns 存储后续附加的其他非结构化信息
#data.obs_names 和 adata.var_names index
#细胞名和基因名可分别通过 adata.obs_names 和 adata.var_names 查看。 AnnData 对象可以像 dataframe 一样进行切片操作，例如，data_subset = data[:, list_of_gene_names]

# 当然scanpy可以直接读取10Xgenomics的.h5格式数据
data=sc.read_10x_h5("./pbmc3K.h5",genome=None,gex_only=True)
adata.var_names_make_unique()`

数据预处理（基因细胞过滤）：

可视化所有细胞中计数最多的基因。我们可以看到前20个基因中有14个基因是属于核糖体基因，说明核糖体基因在这此数据集中表达丰度很高，有时候我们会去除核糖体占比过高的细胞，认为这些细胞质量低，这儿我们先不考虑核糖体基因，考虑线粒体相关基因。

sc.pl.highest_expr_genes(data, n_top=20)
plt.savefig("Highest_expr_genes.pdf")
#基础过滤：去除表达基因200以下的细胞；去除在3个细胞以下表达的基因。
sc.pp.filter_cells(data, min_genes=200)
sc.pp.filter_genes(data, min_cells=3)`

高表达基因top 20

计算线粒体基因占所有基因的比例

mito_genes=data.var_names.str.startswith('MT-')
data.obs['percent_mito']=np.sum(data[:,mito_genes].X,axis=1).A1/np.sum(data.X,axis=1).A1
data.obs['n_counts']=data.X.sum(axis=1).A1
sc.pl.violin(data, ['n_genes', 'n_counts', 'percent_mito'],jitter=0.4, multi_panel=True)
plt.savefig("QC_violin.pdf")
sc.pl.scatter(data, x='n_counts', y='percent_mito')
plt.savefig("QC_dot.pdf")
sc.pl.scatter(data,x='n_counts',y='n_genes')
plt.savefig("QC_dot1.pdf")
##过滤线粒体基因比例 > 5% 和基因总数 >2500 的细胞。
data = data[data.obs.n_genes < 2500, :]
data = data[data.obs.percent_mito < 0.05, :]

数据预处理（标准化，挑选HVG基因）

sc.pp.normalize_total(data, target_sum=1e4) ##标准化
sc.pp.log1p(data)
data.raw = data
sc.pp.highly_variable_genes(data, min_mean=0.0125, max_mean=3, min_disp=0.5)
sc.pl.highly_variable_genes(data)
plt.savefig("highly_variable_genes.pdf")
data = data[:, data.var.highly_variable] 
sc.pp.regress_out(data, ['n_counts', 'percent_mito']) #校正细胞基因计数和线粒体基因比例的影响。
sc.pp.scale(data, max_value=10)

HVG基因的挑选

PCA降维，聚类，Umap可视化

sc.tl.pca(data, svd_solver='arpack')# svd_solver 指定奇异值分解 SVD 的方法
sc.pl.pca_variance_ratio(data, log=True)
plt.savefig("PCA.pdf")
data.write(results_file)
sc.pp.neighbors(data, n_neighbors=10, n_pcs=40)
sc.tl.leiden(data) #使用leiden进行聚类，注意安装对应的python包，conda install -c conda-forge leidenalg ，当然也可使用其他的聚类算法，如louvain，sc.tl.louvain(data),sc.pl.umap(adata, color=['louvain'])，比较了一下，聚类结果差异不大
sc.tl.umap(data)
sc.pl.umap(data, color=['leiden'])
plt.savefig("Umap.pdf")
data.write(results_file)

PC得分

预测簇 marker genes

sc.tl.rank_genes_groups(data, 'leiden', method='wilcoxon')
sc.pl.rank_genes_groups(data, n_genes=30, sharey=False)
plt.savefig("dif_gene.pdf")
data.write(results_file)

cluster 相关差异表达基因

umap图

data = sc.read(results_file)
pd.DataFrame(data.uns['rank_genes_groups']['names']).head(5)
data.uns['rank_genes_groups'].keys()
#dict_keys(['logfoldchanges', 'names', 'params', 'pvals', 'pvals_adj', 'scores'])
result = data.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)
res = pd.DataFrame(    {group + '_' + key: result[key][group]    for group in groups for key in ['names', 'pvals','logfoldchanges','pvals_adj','scores']})
res.to_csv("dif.csv") #基因差异情况输出到本地保存

Top5 marker genes

Top5 marker genes

当然任意两个簇之间也可以比较差异:

sc.tl.rank_genes_groups(data, 'leiden', groups=['0'], reference='1', method='wilcoxon') #cluster0 与 cluster1之间进行比较
sc.pl.rank_genes_groups(data, groups=['0'], n_genes=20)
plt.savefig("gene_exp_C0_C1_score.pdf")
sc.pl.rank_genes_groups_violin(data, groups='0', n_genes=8)
plt.savefig("gene_exp_C0_C1_vio.pdf")
data = sc.read(results_file)
sc.pl.rank_genes_groups_violin(data, groups='0', n_genes=8)
plt.savefig("gene_exp_C0_all_vio.pdf")

任意基因表达可视化(小提琴图、umap)

sc.pl.violin(data, ['MS4A1', 'NKG7', 'PPBP'], groupby='leiden')
plt.savefig("marker_exp_Vio.pdf")
sc.pl.umap(data, color=['MS4A1', 'NKG7', 'PPBP'])
plt.savefig("marker_exp_umap.pdf")

基因表达小提琴图

基因表达Umap图

从之前的Seurat分析中，我们可以了解pbmc包含的一些细胞类别的已知marker gene：

pbmc已知的细胞类别及marker genes

我们可以清晰的看到B细胞相关的MS4A1在簇2中表达，因此我们可以定于这簇为B细胞，同样的NKG7在簇4表达最高，簇4为NK细胞。

确定细胞类别，打上细胞类别标签，绘制已知marker gene的小提琴图

marker_genes = ['IL7R', 'CD14', 'LYZ',  'MS4A1', 'CD8A', 'CD8B', 'FCGR3A', 'MS4A7','GNLY', 'NKG7', 'FCER1A', 'CST3', 'PPBP']
new_cluster_names = ['CD4 T', 'CD14 Monocytes','B', 'CD8 T','NK', 'FCGR3A Monocytes','Dendritic', 'Megakaryocytes']
data.rename_categories('leiden', new_cluster_names)
sc.pl.umap(data, color='leiden', legend_loc='on data', title='', frameon=True)
plt.savefig("umap_celltype.pdf")
ax = sc.pl.dotplot(data, marker_genes, groupby='leiden')
plt.savefig("gene_Dot_celltype.pdf")
ax = sc.pl.stacked_violin(data, marker_genes, groupby='leiden', rotation=90)
plt.savefig("gene_Vio_celltype.pdf")

Umap -celltype

用已知的细胞类别标签替换原始分簇。

在不同细胞类别中marker gene表达点图

已知的marker genes在类别中的表达情况，我们可以明显的看出这些marker genes在不同簇中特异表达。

在不同细胞类别中marker gene表达小提琴图