這篇文章學(xué)習(xí)scanpy官網(wǎng)的第二部分，這部分介紹了如何使用scanpy進(jìn)行軌跡推斷。
官網(wǎng)地址：https://scanpy-tutorials.readthedocs.io/en/latest/paga-paul15.html

#load數(shù)據(jù)，這部分學(xué)習(xí)的數(shù)據(jù)不用下載，貌似是scanpy自帶的
>>> 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()
scanpy==1.4.6 anndata==0.7.1 umap==0.4.1 numpy==1.18.2 scipy==1.4.1 pandas==1.0.3 scikit-learn==0.22.2.post1 statsmodels==0.11.1 python-igraph==0.8.0
>>> results_file = './write/paul15.h5ad'
>>> adata = sc.datasets.paul15()
>>> adata
AnnData object with n_obs × n_vars = 2730 × 3451 
    obs: 'paul15_clusters'
    uns: 'iroot'

數(shù)據(jù)處理與可視化

這里的數(shù)據(jù)處理官網(wǎng)用了scnapy里的一種自帶的處理過程，你也可以使用上一篇文章里的數(shù)據(jù)預(yù)處理方法。關(guān)于這個zheng17的方法的具體代碼，可以看單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析|| scanpy教程：PAGA軌跡推斷。這里我就不贅述了。

>>> sc.pp.recipe_zheng17(adata)
running recipe zheng17
normalizing counts per cell
    finished (0:00:00)
extracting highly variable genes
    finished (0:00:00)
normalizing counts per cell
    finished (0:00:00)
    finished (0:00:00)

跑PCA（降維）：

>>> sc.tl.pca(adata, svd_solver='arpack')
computing PCA with n_comps = 50
    finished (0:00:00)
#計算neighbor graph
>>> sc.pp.neighbors(adata, n_neighbors=4, n_pcs=20)
computing neighbors
    using 'X_pca' with n_pcs = 20
    finished: added to `.uns['neighbors']`
    'distances', distances for each pair of neighbors
    'connectivities', weighted adjacency matrix (0:00:02)
>>> sc.tl.draw_graph(adata)
#出圖
>>> sc.pl.draw_graph(adata, color='paul15_clusters', legend_loc='on data')

降低噪音（Denoising the graph）

上圖看起來是不是很亂？
為了讓上圖看起來有序一點，我們試著用另一種方法進(jìn)行降維：diffusion map。關(guān)于diffusion map降維的介紹，可以參考我之前看視頻做的筆記Single cell RNA-seq data analysis with R視頻學(xué)習(xí)筆記（八）。

>>> sc.tl.diffmap(adata)
computing Diffusion Maps using n_comps=15(=n_dcs)
computing transitions
    finished (0:00:00)
    eigenvalues of transition matrix
    [1.         1.         0.9989645  0.9967852  0.9944013  0.98928535
     0.9882636  0.98712575 0.98383176 0.98297554 0.9789326  0.97689945
     0.9744091  0.9727858  0.9661876 ]
    finished: added
    'X_diffmap', diffmap coordinates (adata.obsm)
    'diffmap_evals', eigenvalues of transition matrix (adata.uns) (0:00:00)
>>> sc.pp.neighbors(adata, n_neighbors=10, use_rep='X_diffmap')
computing neighbors
    finished: added to `.uns['neighbors']`
    'distances', distances for each pair of neighbors
    'connectivities', weighted adjacency matrix (0:00:00)
>>> sc.tl.draw_graph(adata)
drawing single-cell graph using layout 'fa'
WARNING: Package 'fa2' is not installed, falling back to layout 'fr'.To use the faster and better ForceAtlas2 layout, install package 'fa2' (`pip install fa2`).
    finished: added
    'X_draw_graph_fr', graph_drawing coordinates (adata.obsm) (0:00:12)
>>> sc.pl.draw_graph(adata, color='paul15_clusters', legend_loc='on data')

看起來依然很亂。但官網(wǎng)給出的解釋是：有些分化過程的分支被過度繪制了。

聚類和PAGA

這里用louvain來進(jìn)行聚類（起始這里不太理解的是，上一步實際上已經(jīng)聚類了，而且還標(biāo)記了細(xì)胞類型，但官網(wǎng)這里仍然進(jìn)行了聚類）。
PAGA可以生成粗粒度的可視化圖像 （coarse‐grained visualizations），從而可以簡化單細(xì)胞數(shù)據(jù)的解釋，尤其是在測序細(xì)胞量大或整合了大量細(xì)胞的情況下。（參考：https://zhuanlan.zhihu.com/p/108918012）

>>> sc.tl.louvain(adata, resolution=1.0)
>>> sc.tl.paga(adata, groups='louvain')

>>> sc.pl.paga(adata, color=['louvain', 'Hba-a2', 'Elane', 'Irf8'])

這一步我的電腦報錯了，顯示AttributeError: module 'matplotlib.cbook' has no attribute 'is_numlike'。如果你出現(xiàn)了相同的報錯，可以嘗試pip unstall matplotlib下載，然后安裝低版本的pip install matplotlib ==2.2.3

>>> sc.pl.paga(adata, color=['louvain', 'Itga2b', 'Prss34', 'Cma1'])

下面對cluster進(jìn)行注釋：

>>>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'],
      dtype='object')
>>> adata.obs['louvain_anno'] = adata.obs['louvain']
#下面在對cluster進(jìn)行注釋的時候，你需要提前知道哪一個cluster表達(dá)干細(xì)胞的基因，或者表達(dá)lineage特異基因。這里我就先隨便標(biāo)注了
>>> adata.obs['louvain_anno'].cat.categories = ['1/Stem', '2', '3', '4', '5', '6', '7', '8', '9', '10/Ery', '11', '12','13', '14', '15', '16', '17', '18', '19/Neu', '20/Mk', '21', '22/Baso', '23', '24/Mo']
#注釋
>>> sc.tl.paga(adata, groups='louvain_anno')
running PAGA
    finished: added
    'paga/connectivities', connectivities adjacency (adata.uns)
    'paga/connectivities_tree', connectivities subtree (adata.uns) (0:00:00)
#出圖
>>> sc.pl.paga(adata, threshold=0.03)

利用PAGA初始化重新計算embedding

>>> sc.tl.draw_graph(adata, init_pos='paga')
#畫marker基因
>>> sc.pl.draw_graph(adata, color=['louvain_anno', 'Itga2b', 'Prss34', 'Cma1'], legend_loc='on data')

把上圖改顏色：

>>> 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()

這樣新的顏色，每一種都有編號

>>> zeileis_colors = np.array(sc.pl.palettes.zeileis_28)
>>> new_colors = np.array(adata.uns['louvain_anno_colors'])

#把擬時間上每一個點重新分配顏色
>>> new_colors[[0]] = zeileis_colors[[12]]  
>>> new_colors[[2, 4, 12, 15, 11, 3, 7, 18, 10]] = zeileis_colors[[2, 3, 9, 10, 10, 11, 11, 5, 5]] 
>>> new_colors[[8, 20]] = zeileis_colors[[16, 17]] 
>>> new_colors[[17]] = zeileis_colors[[25]]  
>>> new_colors[[16]] = zeileis_colors[[18]]  
>>> new_colors[[19, 5, 6, 9]] = zeileis_colors[[8, 7, 6, 0]]  
>>> new_colors[[1, 13, 14, 21]] = zeileis_colors[[27, 27, 27, 27]]  
>>> new_colors[[22, 23]] = zeileis_colors[[1, 1]] 

>>> adata.uns['louvain_anno_colors'] = new_colors
>>> 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)

利用已知的基因集重構(gòu)PAGA Paths上的基因變化

首先確定擬時間上的root:

>>> adata.uns['iroot'] = np.flatnonzero(adata.obs['louvain_anno']  == '1/Stem')[0]
>>> sc.tl.dpt(adata)
computing Diffusion Pseudotime using n_dcs=10
    finished: added
    'dpt_pseudotime', the pseudotime (adata.obs) (0:00:00)

#給定一組已知的marker基因
#Select some of the marker gene names.
>>> gene_names = ['Gata2', 'Gata1', 'Klf1', 'Epor', 'Hba-a2',  # erythroid
              'Elane', 'Cebpe', 'Gfi1',                    # neutrophil
              'Irf8', 'Csf1r', 'Ctsg']                     # monocyte
#Use the full raw data for visualization.
>>> 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')

#你可以定義每一條lineage是通過哪一條途徑發(fā)育的
>>> paths = [('erythrocytes', [1, 3, 5, 13, 12, 16, 4, 8]),
         ('neutrophils', [1, 20, 2, 14, 15, 22]),
         ('monocytes', [1, 20, 6, 7, 10])]
>>> 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()