參考官網(wǎng) https://scanpy-tutorials.readthedocs.io/en/latest/paga-paul15.html 學(xué)習(xí)
重建髓樣和紅系分化 數(shù)據(jù)為 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]
數(shù)據(jù)是包里用一下命運引用，不用自己再下載哦

adata = sc.datasets.paul15()

[4]

adata

以比默認(rèn)的“ float32”更高的精度進(jìn)行工作，以確保在不同的計算平臺上獲得完全相同的結(jié)果。
[5]

adata.X = adata.X.astype('float64')  # this is not required and results will be comparable without it

預(yù)處理和可視化

可參考之前講解的預(yù)處理方法也可以參考此應(yīng)用簡單的預(yù)處理方法。
[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')

看起來很亂
需要降噪

選擇性：對圖形進(jìn)行降噪

為了使圖去噪，在擴散圖空間（而不是在PCA空間）中表示它。計算幾個擴散分量內(nèi)的距離就等于對圖進(jìn)行去噪-只考慮第一個一些分量。這與使用PCA去噪數(shù)據(jù)矩陣非常相似。該方法已在幾篇論文中使用，例如參見Schiebinger（2017）或Tabaka（2018）。這也與MAGIC Dijk背后的原理有關(guān)（2018）。

對于PAGA，聚類或偽時間估計，這都不是必需的步驟。可以繼續(xù)使用非去噪圖。在許多情況下（也在此處），這將為您帶來非常不錯的效果。
[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，現(xiàn)在 使用sc.tl.louvain為了復(fù)現(xiàn)論文結(jié)果。
[13]

sc.tl.louvain(adata, resolution=1.0)

使用基因marker 注釋細(xì)胞群

對于簡單的粗粒度可視化，計算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'])

實際注釋細(xì)胞 -注意Cma1是肥大細(xì)胞標(biāo)記，僅出現(xiàn)在祖細(xì)胞/干細(xì)胞簇8的一小部分細(xì)胞中，參見下面的單細(xì)胞分解圖。
[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')

現(xiàn)在可以在有意義的布局中以單細(xì)胞分辨率查看所有基因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路徑重建基因變化

選擇一個根細(xì)胞群進(jìn)行偽時間擴散
[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

使用完整的原始數(shù)據(jù)進(jìn)行可視化
[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()