Scanpy 是一个基于 Python 单细胞数据分析软件包，内容包括预处理，可视化，聚类，拟时序分析和差异表达分析等。在单细胞数据过多时，使用R进行一些单细胞分析比如monocle等即使使用服务器会出现内存不足的情况，而Scanpy则能很好的解决这个问题。
官网：https://scanpy-tutorials.readthedocs.io/en/latest/pbmc3k.html
0. 数据准备
确保python3已装好
下载scanpy
pip install scanpy
下载练习数据集，还是熟悉的pbmc3k
wget https://cf.10xgenomics.com/samples/cell/pbmc3k/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 数据
adata=sc.read_10x_mtx('./filtered_gene_bc_matrices/hg19/',var_names='gene_symbols', cache=True)
adata.var_names_make_unique() # 索引去重，若上一步中使用 `var_names='gene_ids'` 则这一步非必须进行
#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 一样进行切片操作，例如，data_subset = adata[:, list_of_gene_names]
# 当然scanpy可以直接读取10Xgenomics的.h5格式数据
adata=sc.read_10x_h5("./pbmc3K.h5",genome=None,gex_only=True)
adata.var_names_make_unique()
1. 数据预处理（基因细胞过滤）：
# 可视化所有细胞中计数最多的基因
sc.pl.highest_expr_genes(adata, n_top=20)
plt.savefig("Highest_expr_genes.pdf")
可以看到前20个基因中有14个基因是属于核糖体基因RP-，说明核糖体基因在这此数据集中表达丰度很高，有时候我们会去除核糖体占比过高的细胞，认为这些细胞质量低，这儿我们先不考虑核糖体基因和线粒体相关基因。
#基础过滤：去除表达基因200以下的细胞；去除在3个细胞以下表达的基因。
sc.pp.filter_cells(data, min_genes=200)
sc.pp.filter_genes(data, min_cells=3)
计算线粒体基因占所有基因的比例
adata.var['mt'] = adata.var_names.str.startswith('MT-') # annotate the group of mitochondrial genes as '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)
plt.savefig("QC_violin.pdf")
sc.pl.scatter(adata, x='total_counts', y='pct_counts_mt')
sc.pl.scatter(adata, x='total_counts', y='n_genes_by_counts')
过滤
##过滤线粒体基因比例 > 5% 和基因总数 >2500 的细胞。
adata = adata[adata.obs.n_genes_by_counts < 2500, :]
adata = adata[adata.obs.pct_counts_mt < 5, :]
数据预处理（标准化，挑选HVG基因）
这三步就相当于Seurat中的NormalizeData()和FindVariableFeatures()
sc.pp.normalize_total(adata, target_sum=1e4) ##标准化
sc.pp.log1p(adata) #归一化
sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5) #鉴定高变基因
sc.pl.highly_variable_genes(adata)
Scale数据，相当于Seurat中的ScaleData()
adata.raw = adata #备份
adata = adata[:, adata.var.highly_variable] #将保守基因去除，留下差异表达的基因用于后续分析
sc.pp.regress_out(adata, ['total_counts', 'pct_counts_mt']) #占内存。校正细胞基因计数和线粒体基因比例的影响。
sc.pp.scale(adata, max_value=10) #将数据放缩到方差为1
2. PCA降维，聚类，Umap可视化
主成分分析是一种将数据降维的分析方法，是考察多个变量间相关性一种多元统计方法，研究如何通过少数几个主成分来揭示多个变量间的内部结构，即从原始变量中导出少数几个主成分，使它们尽可能多地保留原始变量的信息，且彼此间互不相关.通常数学上的处理就是将原来P个指标作线性组合，作为新的综合指标。
sc.tl.pca(adata, svd_solver='arpack')# svd_solver 指定奇异值分解 SVD 的方法
sc.pl.pca(adata, color='CST3')
#碎石图观测主成分的质量，用于选择后续应该使用多少个PC，用于计算细胞间的相邻距离。
sc.pl.pca_variance_ratio(adata, log=True)
adata.write(results_file) #保存数据
adata
3. 聚类分析
非线性降维
# n_neighbors指的是每个点的邻近点的数量。neighbors的个数越多，聚类数会越少。
sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40)
# tl.paga(adata)
# pl.paga(adata, plot=False) # remove `plot=False` if you want to see the coarse-grained graph
# tl.umap(adata, init_pos='paga')
# 将距离嵌入图中
sc.tl.umap(adata)
sc.pl.umap(adata, color=['CST3', 'NKG7', 'PPBP'], use_raw=False)
聚类
这里用的是leiden聚类算法。
Seurat的函数FindClusters()有一个参数algorithm，默认是1。也就是说
Louvain是Seurat做聚类时的默认算法。而设置为4时，就是Leiden算法，是Scanpy做聚类的默认算法。Louvain算法有三个缺点：1. 社区划分的精度有局限性；2. 分组内细胞分布密度的大小会影响亚群的鉴定；3. 被鉴定为同一个分群的细胞群内，存在两个没有连线的小分群。 Leiden算法主要针对上述的第3个缺点，对louvain算法进行优化。参考：Leiden算法对Louvain算法的优化。
sc.tl.leiden(adata)
sc.pl.umap(adata, color=['leiden', 'CST3', 'NKG7'])
adata.write(results_file)
得到8个cluster
4. 寻找marker基因用于细胞群注释
这里用的是wilcoxon（推荐），也可以用t.test，logreg等。
#sc.tl.rank_genes_groups(adata, 'leiden', method='t-test')
#sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)
sc.settings.verbosity = 2 # reduce the verbosity
sc.tl.rank_genes_groups(adata, 'leiden', method='wilcoxon')
sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False) #绘制每个群top25的marker基因
adata.write(results_file)
##定义marker基因list用于后续画图
marker_genes = ['IL7R', 'CD79A', 'MS4A1', 'CD8A', 'CD8B', 'LYZ', 'CD14',
'LGALS3', 'S100A8', 'GNLY', 'NKG7', 'KLRB1',
'FCGR3A', 'MS4A7', 'FCER1A', 'CST3', 'PPBP']
adata = sc.read(results_file) #读入计算的marker基因文件
#数据框展示0-7群的top 10的marker基因
pd.DataFrame(adata.uns['rank_genes_groups']['names']).head(5)
#Get a table with the scores and groups.
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，寻找差异基因
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=8) #小提琴图展示
同一基因在不同组间的小提琴图
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.obs['leiden_anno'] = adata.obs['leiden']
adata.rename_categories('leiden_anno', new_cluster_names)
细胞群的marker基因展示
sc.pl.dotplot(adata, marker_genes, groupby='leiden_anno');
sc.pl.stacked_violin(adata, marker_genes, groupby='leiden_anno', rotation=90);
adata.write('./pbmc3k.h5ad')
