用VSCode Jupyter 学习Scanpy——轨迹推断
2020-12-20 本文已影响0人
米妮爱分享
参考官网 https://scanpy-tutorials.readthedocs.io/en/latest/paga-paul15.html 学习
重建髓样和红系分化 数据为 Paul et al. (2015).
引入包
[1]
import numpy as np
import pandas as pd
import matplotlib.pyplot as pl
from matplotlib import rcParams
import scanpy as sc
[2]
sc.settings.verbosity = 3 # verbosity: errors (0), warnings (1), info (2), hints (3)
sc.logging.print_header()
results_file = './write/paul15.h5ad'
sc.settings.set_figure_params(dpi=80, facecolor='white')
[3]
数据是包里用一下命运引用，不用自己再下载哦
adata = sc.datasets.paul15()
[4]
adata
以比默认的“ float32”更高的精度进行工作，以确保在不同的计算平台上获得完全相同的结果。
[5]
adata.X = adata.X.astype('float64') # this is not required and results will be comparable without it
预处理和可视化
可参考之前讲解的预处理方法也可以参考此应用简单的预处理方法。
[6]
sc.pp.recipe_zheng17(adata)
[7]
sc.tl.pca(adata, svd_solver='arpack')
[8]
sc.pp.neighbors(adata, n_neighbors=4, n_pcs=20)
sc.tl.draw_graph(adata)
[9]
sc.pl.draw_graph(adata, color='paul15_clusters', legend_loc='on data')
看起来很乱
需要降噪
选择性：对图形进行降噪
为了使图去噪，在扩散图空间（而不是在PCA空间）中表示它。计算几个扩散分量内的距离就等于对图进行去噪-只考虑第一个一些分量。这与使用PCA去噪数据矩阵非常相似。该方法已在几篇论文中使用，例如参见Schiebinger（2017）或Tabaka（2018）。这也与MAGIC Dijk背后的原理有关（2018）。
对于PAGA，聚类或伪时间估计，这都不是必需的步骤。可以继续使用非去噪图。在许多情况下（也在此处），这将为您带来非常不错的效果。
[10]
sc.tl.diffmap(adata)
sc.pp.neighbors(adata, n_neighbors=10, use_rep='X_diffmap')
[11]
sc.tl.draw_graph(adata)
[12]
sc.pl.draw_graph(adata, color='paul15_clusters', legend_loc='on data')
这看起来仍然很杂乱，但是方式却不同：许多分支都被过度绘制。
分群 和 PAGA
请注意，在这里使用了sc.tl.leiden，现在 使用sc.tl.louvain为了复现论文结果。
[13]
sc.tl.louvain(adata, resolution=1.0)
使用基因marker 注释细胞群
| cell type | marker |
|---|---|
| HSCs | Procr |
| Erythroids | Gata1, Klf1, Epor, Gypa, Hba-a2, Hba-a1, Spi1 |
| Neutrophils | Elane, Cebpe, Ctsg, Mpo, Gfi1 |
| Monocytes | Irf8, Csf1r, Ctsg, Mpo |
| Megakaryocytes | Itga2b (encodes protein CD41), Pbx1, Sdpr, Vwf |
| Basophils | Mcpt8, Prss34 |
| B cells | Cd19, Vpreb2, Cd79a |
| Mast cells | Cma1, Gzmb, CD117/C-Kit |
| Mast cells & Basophils | Ms4a2, Fcer1a, Cpa3, CD203c (human) |
对于简单的粗粒度可视化，计算PAGA图，粗粒度和简化（抽象）的图形。粗粒度图中的非重要边缘将被阈值化。
[14]
sc.tl.paga(adata, groups='louvain')
[15]
sc.pl.paga(adata, color=['louvain', 'Hba-a2', 'Elane', 'Irf8'])
[16]
sc.pl.paga(adata, color=['louvain', 'Itga2b', 'Prss34', 'Cma1'])
实际注释细胞 -注意Cma1是肥大细胞标记，仅出现在祖细胞/干细胞簇8的一小部分细胞中，参见下面的单细胞分解图。
[17]
adata.obs['louvain'].cat.categories
[18]
adata.obs['louvain_anno'] = adata.obs['louvain']
[19]
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']
对PAGA使用带注释的群集
[20]
sc.tl.paga(adata, groups='louvain_anno')
[21]
sc.pl.paga(adata, threshold=0.03, show=False)
使用PAGA初始化重计算嵌入
对于UMAP而言，以下是可能的
[22]
sc.tl.draw_graph(adata, init_pos='paga')
现在可以在有意义的布局中以单细胞分辨率查看所有基因marker。
[23]
sc.pl.draw_graph(adata, color=['louvain_anno', 'Itga2b', 'Prss34', 'Cma1'], legend_loc='on data')
选择与群集更加一致的颜色
[24]
pl.figure(figsize=(8, 2))
for i in range(28):
pl.scatter(i, 1, c=sc.pl.palettes.zeileis_28[i], s=200)
pl.show()
[25]
zeileis_colors = np.array(sc.pl.palettes.zeileis_28)
new_colors = np.array(adata.uns['louvain_anno_colors'])
[26]
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
[27]
adata.uns['louvain_anno_colors'] = new_colors
并在某些群集名称中添加一些空格。这里显示的布局不同于原paper的布局，可以在此处看到。但是这些差异只是表面上的。从随机PCA和float32精度提高到float64精度，不得不更改布局。
[28]
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)
针对给定的一组基因，沿PAGA路径重建基因变化
选择一个根细胞群进行伪时间扩散
[29]
adata.uns['iroot'] = np.flatnonzero(adata.obs['louvain_anno'] == '16/Stem')[0]
[30]
sc.tl.dpt(adata)
选择一些基因marker名称
[31]
gene_names = ['Gata2', 'Gata1', 'Klf1', 'Epor', 'Hba-a2', # erythroid
'Elane', 'Cebpe', 'Gfi1', # neutrophil
'Irf8', 'Csf1r', 'Ctsg'] # monocyte
使用完整的原始数据进行可视化
[31]
adata_raw = sc.datasets.paul15()
sc.pp.log1p(adata_raw)
sc.pp.scale(adata_raw)
adata.raw = adata_raw
[33]
sc.pl.draw_graph(adata, color=['louvain_anno', 'dpt_pseudotime'], legend_loc='on data')
[34]
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])]
[35]
adata.obs['distance'] = adata.obs['dpt_pseudotime']
[36]
adata.obs['clusters'] = adata.obs['louvain_anno'] # just a cosmetic change
[37]
adata.uns['clusters_colors'] = adata.uns['louvain_anno_colors']
[38]
!mkdir write
[39]
_, 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()
