hi，大家好，好久不见，这次跟大家分享一个单细胞降维聚类的新的分析方法scanpy，目前大部多数文章用的单细胞分析均用的Seurat分析包，目前已经更新到了3.0版本，在原有的基础上进行了优化，最大的变化就是在T-SNE方法的基础上添加了U-map的降维分析方法，相信很多道友已经尝试过了，不知道是否达到了预期呢？今天我们就来分享一下python带来的新方法，也许，就是你想要的。
第一步，加载模块
import numpy as np
import pandas as pd
import matplotlib.pyplot as pl
from matplotlib import rcParams
import scanpy as sc
sc.settings.verbosity = 3 # verbosity: errors (0), warnings (1), info (2), hints (3)
sc.logging.print_versions()
results_file = './write/paul15.h5ad'
这几个模块除了scanpy，其他的模块均是python常用的模块，大家有兴趣的话可以查看一下我的另外一篇文章----python常用模块，这里就不多介绍了。
scanpy==1.3.7+86.g2c80c7a anndata==0.6.17+1.ga0cd0c6 numpy==1.14.6 scipy==1.1.0 pandas==0.23.4 scikit-learn==0.19.1 statsmodels==0.9.0 python-igraph==0.7.1 louvain==0.6.1
sc.settings.set_figure_params(dpi=80, frameon=False) # low dpi (dots per inch) yields small inline figures
adata = sc.datasets.paul15()
adata
会看到以下结果
Out[4]:
AnnData object with n_obs × n_vars = 2730 × 3451
obs: 'paul15_clusters'
uns: 'iroot'
让我们以比默认“float32”更高的准确度运行，以确保在不同的计算平台上获得完全相同的结果。
adata.X = adata.X.astype('float64') # this is not required and results will be comparable without it
第二步 预处理和可视化
读取到数据矩阵之后，让我们做后续的分析
sc.pp.recipe_zheng17(adata)
sc.tl.pca(adata, svd_solver='arpack')
sc.pp.neighbors(adata, n_neighbors=4, n_pcs=20)
sc.tl.draw_graph(adata)
#运行经过
computing neighbors
using 'X_pca' with n_pcs = 20
finished (0:00:02.65) --> added to `.uns['neighbors']`
'distances', distances for each pair of neighbors
'connectivities', weighted adjacency matrix
drawing single-cell graph using layout "fa"
finished (0:00:11.41) --> added
'X_draw_graph_fa', graph_drawing coordinates (adata.obsm)
sc.pl.draw_graph(adata, color='paul15_clusters', legend_loc='on data')
会得到以下的结果：
123.png
这看起来像一团糟。
第三步 （可选）对图表进行去除背景噪音
为了对图去除背景噪音，我们在扩散映射空间（而不是在PCA空间中）表示它。计算几个扩散分量内的距离相当于对图进行去噪 - 我们只需要考虑一些第一个光谱分量。 这与使用PCA对数据矩阵进行去噪非常相似。 该方法已被用于几篇论文中，参见例如 Schiebinger等人。 （2017）或Tabaka等。（2018）。 它也与MAGIC背后的原理有关
sc.tl.diffmap(adata)
sc.pp.neighbors(adata, n_neighbors=10, use_rep='X_diffmap')
#运行过程
computing Diffusion Maps using n_comps=15(=n_dcs)
eigenvalues of transition matrix
[1. 1. 0.9989278 0.99671 0.99430376 0.98939794
0.9883687 0.98731077 0.9839871 0.983007 0.97908056 0.97625476
0.9744365 0.9729161 0.9652972 ]
finished (0:00:00.07) --> added
'X_diffmap', diffmap coordinates (adata.obsm)
'diffmap_evals', eigenvalues of transition matrix (adata.uns)
computing neighbors
finished (0:00:00.46) --> added to `.uns['neighbors']`
'distances', distances for each pair of neighbors
'connectivities', weighted adjacency matrix
sc.tl.draw_graph(adata)
#运行过程
drawing single-cell graph using layout "fa"
finished (0:00:10.23) --> added
'X_draw_graph_fa', graph_drawing coordinates (adata.obsm)
sc.pl.draw_graph(adata, color='paul15_clusters', legend_loc='on data')
可以得到这样的结果
234.png
这仍然看起来很混乱，但是以不同的方式：很多分支被过度绘制。
第四步 Clustering and PAGA
sc.tl.louvain(adata, resolution=1.0)
#运行过程
running Louvain clustering
using the "louvain" package of Traag (2017)
finished (0:00:00.10) --> found 25 clusters and added
'louvain', the cluster labels (adata.obs, categorical)
使用标记基因注释cluster。
11.png
对于简单的粗粒度可视化，计算PAGA图，粗粒度和简化（抽象）图。 粗粒度图中的非重要边缘被阈值化。
sc.tl.paga(adata, groups='louvain')
#运行过程
running PAGA
finished (0:00:00.14) --> added
'paga/connectivities', connectivities adjacency (adata.uns)
'paga/connectivities_tree', connectivities subtree (adata.uns)
sc.pl.paga(adata, color=['louvain', 'Hba-a2', 'Elane', 'Irf8'])
得到如下结果：
1.png
sc.pl.paga(adata, color=['louvain', 'Itga2b', 'Prss34', 'Cma1'])
2.png
实际上注释cluster - 注意Cma1是肥大细胞标记并且仅出现在祖细胞/干细胞簇8中的一小部分细胞中，参见下面的单细胞解析图。
adata.obs['louvain'].cat.categories
#运行过程
Index(['0', '1', '2', '3', '4', '5', '6', '7', '8', '9', '10', '11', '12',
'13', '14', '15', '16', '17', '18', '19', '20', '21', '22', '23', '24'],
dtype='object')
adata.obs['louvain_anno'] = adata.obs['louvain']
adata.obs['louvain_anno'].cat.categories = ['0', '1', '2', '3', '4', '5', '6', '7', '8', '9', '10/Ery', '11', '12',
'13', '14', '15', '16/Stem', '17', '18', '19/Neu', '20/Mk', '21', '22/Baso', '23', '24/Mo']
让我们使用带注释的cluster进行PAGA。
sc.tl.paga(adata, groups='louvain_anno')
#run
running PAGA
finished (0:00:00.11) --> added
'paga/connectivities', connectivities adjacency (adata.uns)
'paga/connectivities_tree', connectivities subtree (adata.uns)
sc.pl.paga(adata, threshold=0.03)
会得到新的结果：
23.png
第五步 使用PAGA初始化重新计算嵌入
对于UMAP，以下内容也是如此
sc.tl.draw_graph(adata, init_pos='paga')
#run
drawing single-cell graph using layout "fa"
finished (0:00:09.98) --> added
'X_draw_graph_fa', graph_drawing coordinates (adata.obsm)
现在我们可以在有意义的布局中以单细胞分辨率看到所有标记基因。
sc.pl.draw_graph(adata, color=['louvain_anno', 'Itga2b', 'Prss34', 'Cma1'], legend_loc='on data')
1.png
更一致地选择簇的颜色。
pl.figure(figsize=(8, 2))
for i in range(28):
pl.scatter(i, 1, c=sc.pl.palettes.zeileis_26[i], s=200)
pl.show()
1.png
zeileis_colors = np.array(sc.pl.palettes.zeileis_26)
new_colors = np.array(adata.uns['louvain_anno_colors'])
new_colors[[16]] = zeileis_colors[[12]] # Stem colors / green
new_colors[[10, 17, 5, 3, 15, 6, 18, 13, 7, 12]] = zeileis_colors[[5, 5, 5, 5, 11, 11, 10, 9, 21, 21]] # Ery colors / red
new_colors[[20, 8]] = zeileis_colors[[17, 16]] # Mk early Ery colors / yellow
new_colors[[4, 0]] = zeileis_colors[[2, 8]] # lymph progenitors / grey
new_colors[[22]] = zeileis_colors[[18]] # Baso / turquoise
new_colors[[19, 14, 2]] = zeileis_colors[[6, 6, 6]] # Neu / light blue
new_colors[[24, 9, 1, 11]] = zeileis_colors[[0, 0, 0, 0]] # Mo / dark blue
new_colors[[21, 23]] = zeileis_colors[[25, 25]] # outliers / grey
adata.uns['louvain_anno_colors'] = new_colors
并为某些群集名称添加一些空白区域。 此处显示的布局与文中的布局不同，可在此处找到。 然而，这些差异只是装饰性的。 当我们从随机PCA和float32移动到float64精度时，我们不得不改变布局
sc.pl.paga_compare(
adata, threshold=0.03, title='', right_margin=0.2, size=10, edge_width_scale=0.5,
legend_fontsize=12, fontsize=12, frameon=False, edges=True, save=True)
2.png
第六步
针对给定基因组重建沿PAGA途径的基因变化
adata.uns['iroot'] = np.flatnonzero(adata.obs['louvain_anno'] == '16/Stem')[0]
sc.tl.dpt(adata)
Select some of the marker gene names.
gene_names = ['Gata2', 'Gata1', 'Klf1', 'Epor', 'Hba-a2', # erythroid
'Elane', 'Cebpe', 'Gfi1', # neutrophil
'Irf8', 'Csf1r', 'Ctsg'] # monocyte
数据可视化
adata_raw = sc.datasets.paul15()
sc.pp.log1p(adata_raw)
sc.pp.scale(adata_raw)
adata.raw = adata_raw
sc.pl.draw_graph(adata, color=['louvain_anno', 'dpt_pseudotime'], legend_loc='on data')
1 (2).png
paths = [('erythrocytes', [16, 12, 7, 13, 18, 6, 5, 10]),
('neutrophils', [16, 0, 4, 2, 14, 19]),
('monocytes', [16, 0, 4, 11, 1, 9, 24])]
adata.obs['distance'] = adata.obs['dpt_pseudotime']
adata.obs['clusters'] = adata.obs['louvain_anno'] # just a cosmetic change
adata.uns['clusters_colors'] = adata.uns['louvain_anno_colors']
!mkdir write
_, axs = pl.subplots(ncols=3, figsize=(6, 2.5), gridspec_kw={'wspace': 0.05, 'left': 0.12})
pl.subplots_adjust(left=0.05, right=0.98, top=0.82, bottom=0.2)
for ipath, (descr, path) in enumerate(paths):
_, data = sc.pl.paga_path(
adata, path, gene_names,
show_node_names=False,
ax=axs[ipath],
ytick_fontsize=12,
left_margin=0.15,
n_avg=50,
annotations=['distance'],
show_yticks=True if ipath==0 else False,
show_colorbar=False,
color_map='Greys',
groups_key='clusters',
color_maps_annotations={'distance': 'viridis'},
title='{} path'.format(descr),
return_data=True,
show=False)
data.to_csv('./write/paga_path_{}.csv'.format(descr))
pl.savefig('./figures/paga_path_paul15.pdf')
pl.show()
1.png
哈哈，小伙伴们，是不是已经觉得山穷水复疑无路，柳暗花明又一村呢？
至于多样本整合的方式，我们留在下次讲解~
网友评论