Scanpy 是一個(gè)基于 Python 分析單細(xì)胞數(shù)據(jù)的軟件包争群，內(nèi)容包括預(yù)處理峰弹，可視化距境，聚類，擬時(shí)序分析和差異表達(dá)分析等垮卓。
https://genomebiology.biomedcentral.com/articles/10.1186/s13059-017-1382-0
https://scanpy.readthedocs.io/en/stable/tutorials.html#integrating-datasets
有人可能會(huì)說(shuō)：?jiǎn)渭?xì)胞分析使用Seurat垫桂，monocle等R包會(huì)更加方便。但是實(shí)際分析中粟按，當(dāng)單細(xì)胞數(shù)據(jù)過(guò)多時(shí)诬滩，Seurat和monocle會(huì)產(chǎn)生內(nèi)存不足的問(wèn)題，尤其是monocle2灭将，單細(xì)胞數(shù)據(jù)大于50k的時(shí)候疼鸟，往往會(huì)產(chǎn)生內(nèi)存不足，即使是在服務(wù)器上運(yùn)行庙曙。而基于python的單細(xì)胞轉(zhuǎn)錄分析包scanpy空镜，能很好的解決內(nèi)存不足的問(wèn)題，根據(jù)官網(wǎng)說(shuō)明捌朴，scanpy可以解決超過(guò)1,000,000個(gè)單細(xì)胞的數(shù)據(jù)集吴攒，這是基于R分析包不敢想象的。(Seurat:https://satijalab.org/seurat/articles/pbmc3k_tutorial.html)

單細(xì)胞轉(zhuǎn)錄數(shù)據(jù)分析之Scanpy：http://www.reibang.com/p/e22a947e6c60
單細(xì)胞轉(zhuǎn)錄組之Scanpy - 軌跡推斷/擬時(shí)序分析：http://www.reibang.com/p/0b2ca0e0b544
單細(xì)胞轉(zhuǎn)錄組之Scanpy - 樣本整合分析：http://www.reibang.com/p/beef8a8be360
單細(xì)胞空間轉(zhuǎn)錄分析之Scanpy：
單細(xì)胞空間轉(zhuǎn)錄分析之Scanpy-結(jié)合單細(xì)胞轉(zhuǎn)錄組：

Fig1. 網(wǎng)址主頁(yè)

這兒將系統(tǒng)全面演示了如何使用Scanpy來(lái)處理單細(xì)胞轉(zhuǎn)錄組學(xué)數(shù)據(jù)砂蔽。

scanpy相關(guān)python 包安裝（安裝好python3之后洼怔，終端運(yùn)行）：

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

數(shù)據(jù)下載，這兒我們使用了Seurat官網(wǎng)使用的pbmc3k數(shù)據(jù)進(jìn)行測(cè)試（我是很早之前就下載過(guò)了）：

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

運(yùn)行python3左驾，導(dǎo)入相關(guān)包：

import scanpy as sc
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
sc.settings.verbosity = 3 # verbosity 的取值表示測(cè)試結(jié)果顯示的詳細(xì)程度镣隶，數(shù)字越大越詳細(xì)
sc.logging.print_versions() # 輸出版本號(hào)
sc.settings.set_figure_params(dpi=80) # set_figure_params 設(shè)置圖片的分辨率/大小以及其他樣式
import os #在服務(wù)器運(yùn)行，習(xí)慣性會(huì)設(shè)置一個(gè)輸出路徑诡右，用于保存pdf圖片
os.getcwd()  #查看當(dāng)前路徑
os.chdir('./filtered_gene_bc_matrices/scanpy') #修改路徑
os.getcwd()
results_file = 'pbmc3k.h5ad' ##置結(jié)果文件保存路徑

讀取并查看數(shù)據(jù)：

# 導(dǎo)入 10X 數(shù)據(jù)
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'` 則這一步非必須進(jìn)行
#data.X 存儲(chǔ) count matrix，數(shù)據(jù)類型為稀疏矩陣 scipy.sparse.csr.csr_matrix
#data.obs 存儲(chǔ)關(guān)于 obervations(cells) 的 metadata帆吻，數(shù)據(jù)類型為 dataframe
#data.var 存儲(chǔ)關(guān)于 variables(genes) 的 metadata域那，數(shù)據(jù)類型為 dataframe
#AnnData.uns 存儲(chǔ)后續(xù)附加的其他非結(jié)構(gòu)化信息
#data.obs_names 和 adata.var_names index
#細(xì)胞名和基因名可分別通過(guò) adata.obs_names 和 adata.var_names 查看。 AnnData 對(duì)象可以像 dataframe 一樣進(jìn)行切片操作桅锄，例如琉雳，data_subset = data[:, list_of_gene_names]

# 當(dāng)然scanpy可以直接讀取10Xgenomics的.h5格式數(shù)據(jù)
data=sc.read_10x_h5("./pbmc3K.h5",genome=None,gex_only=True)
adata.var_names_make_unique()`

數(shù)據(jù)預(yù)處理（基因細(xì)胞過(guò)濾）：

可視化所有細(xì)胞中計(jì)數(shù)最多的基因样眠。我們可以看到前20個(gè)基因中有14個(gè)基因是屬于核糖體基因友瘤，說(shuō)明核糖體基因在這此數(shù)據(jù)集中表達(dá)豐度很高，有時(shí)候我們會(huì)去除核糖體占比過(guò)高的細(xì)胞檐束，認(rèn)為這些細(xì)胞質(zhì)量低辫秧，這兒我們先不考慮核糖體基因，考慮線粒體相關(guān)基因被丧。

sc.pl.highest_expr_genes(data, n_top=20)
plt.savefig("Highest_expr_genes.pdf")
#基礎(chǔ)過(guò)濾：去除表達(dá)基因200以下的細(xì)胞盟戏；去除在3個(gè)細(xì)胞以下表達(dá)的基因绪妹。
sc.pp.filter_cells(data, min_genes=200)
sc.pp.filter_genes(data, min_cells=3)`

高表達(dá)基因top 20

計(jì)算線粒體基因占所有基因的比例

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")
##過(guò)濾線粒體基因比例 > 5% 和基因總數(shù) >2500 的細(xì)胞。
data = data[data.obs.n_genes < 2500, :]
data = data[data.obs.percent_mito < 0.05, :]

數(shù)據(jù)預(yù)處理（標(biāo)準(zhǔn)化柿究，挑選HVG基因）

sc.pp.normalize_total(data, target_sum=1e4) ##標(biāo)準(zhǔn)化
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']) #校正細(xì)胞基因計(jì)數(shù)和線粒體基因比例的影響邮旷。
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進(jìn)行聚類婶肩，注意安裝對(duì)應(yīng)的python包，conda install -c conda-forge leidenalg 貌夕，當(dāng)然也可使用其他的聚類算法律歼，如louvain，sc.tl.louvain(data),sc.pl.umap(adata, color=['louvain'])啡专，比較了一下险毁，聚類結(jié)果差異不大
sc.tl.umap(data)
sc.pl.umap(data, color=['leiden'])
plt.savefig("Umap.pdf")
data.write(results_file)

PC得分

預(yù)測(cè)簇 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 相關(guān)差異表達(dá)基因

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

當(dāng)然任意兩個(gè)簇之間也可以比較差異:

sc.tl.rank_genes_groups(data, 'leiden', groups=['0'], reference='1', method='wilcoxon') #cluster0 與 cluster1之間進(jìn)行比較
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")

任意基因表達(dá)可視化(小提琴圖、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")

基因表達(dá)小提琴圖

基因表達(dá)Umap圖

從之前的Seurat分析中们童，我們可以了解pbmc包含的一些細(xì)胞類別的已知marker gene：

pbmc已知的細(xì)胞類別及marker genes

我們可以清晰的看到B細(xì)胞相關(guān)的MS4A1在簇2中表達(dá)畔况，因此我們可以定于這簇為B細(xì)胞，同樣的NKG7在簇4表達(dá)最高慧库，簇4為NK細(xì)胞问窃。

確定細(xì)胞類別，打上細(xì)胞類別標(biāo)簽完沪，繪制已知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

用已知的細(xì)胞類別標(biāo)簽替換原始分簇域庇。

在不同細(xì)胞類別中marker gene表達(dá)點(diǎn)圖

已知的marker genes在類別中的表達(dá)情況，我們可以明顯的看出這些marker genes在不同簇中特異表達(dá)覆积。

在不同細(xì)胞類別中marker gene表達(dá)小提琴圖