• 学习资料来源：
  • scanpy主页：https://scanpy.readthedocs.io/en/stable/
  • 官网：https://scanpy-tutorials.readthedocs.io/en/latest/paga-paul15.html【注意教程有两个版本，这里是latest版本的学习笔记】

数据说明：

使用来自的Paul et al. (2015)数据重建髓系myeloid和红细胞系erythroid的分化

本次会用到一个新的软件，安装如下：

conda activate scanpy
conda install -c anaconda pandas -y

01 数据导入

首先去：将数据下载下来，这里直接封装到了scanpy包中：sc.datasets.paul15()

import numpy as np
import pandas as pd
import matplotlib.pyplot as pl
from matplotlib import rcParams
import scanpy as sc

# verbosity: errors (0), warnings (1), info (2), hints (3)
sc.settings.verbosity = 3 
sc.logging.print_versions()
# low dpi (dots per inch) yields small inline figures
#sc.settings.set_figure_params(dpi=80, frameon=False, figsize=(3, 3), facecolor='white')  

# 保存数据路径
outdir = './Trajectory/'
adata = sc.datasets.paul15()
adata

# 让我们用比默认的float32更高的精度来确保在不同的计算平台上得到完全相同的结果
# this is not required and results will be comparable without it
adata.X = adata.X.astype('float64')
adata

总共有：2730个细胞，3451个基因

AnnData object with n_obs × n_vars = 2730 × 3451
    obs: 'paul15_clusters'
    uns: 'iroot'

02 预处理和可视化

这里用到了scanpy.pp.recipe_zheng17这个函数，主要是将数据预处理的几个步骤包装成一个函数，处理方式来自文章：

Zheng et al. (2017), Massively parallel digital transcriptional profiling of single cells, Nature Communications

包装步骤包括：

# only consider genes with more than 1 count
sc.pp.filter_genes(adata, min_counts=1)         
# normalize with total UMI count per cell
sc.pp.normalize_per_cell(adata, key_n_counts='n_counts_all')
# select highly-variable genes
filter_result = sc.pp.filter_genes_dispersion(adata.X, flavor='cell_ranger', n_top_genes=n_top_genes, log=False)
# subset the genes
adata = adata[:, filter_result.gene_subset]
# renormalize after filtering
sc.pp.normalize_per_cell(adata) 
# log transform: adata.X = log(adata.X + 1)
if log: sc.pp.log1p(adata)
    # scale to unit variance and shift to zero mean
    sc.pp.scale(adata)                              

进行预处理：

# Apply a simple preprocessing recipe.
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)
sc.pl.draw_graph(adata, color='paul15_clusters', legend_loc='on data',size=8)
pl.savefig(outdir + "01-paul15_clusters.png")
pl.close()

初始轨迹图，This looks pretty messy

1658508738306.png

03 去噪：可选部分

为了去除图的噪声，我们在扩散映射空间(而不是 PCA 空间)中表示它。计算几个扩散成分（diffusion components）内的距离相当于去噪图形-我们只取几个第一个光谱成分（the first spectral components）。这与使用PCA去噪数据矩阵非常相似。

该方法已在几篇论文中使用，如Schiebinger et al. (2017) or Tabaka et al. (2018)。这也与MAGIC Dijk et al. (2018)等人背后的原则有关。

Note：这不是一个必要的步骤，如 PAGA、聚类、伪时间估计等分析都不是一个必要的步骤。你最好使用一个无噪声的图形。在许多情况下(也在这里) ，这将给你带来非常体面的结果。

sc.tl.diffmap(adata)
sc.pp.neighbors(adata, n_neighbors=10, use_rep='X_diffmap')
sc.tl.draw_graph(adata)
sc.pl.draw_graph(adata, color='paul15_clusters', legend_loc='on data')
pl.savefig("./paul15/paul15_test2.png")
pl.close()

这看起来仍然很混乱，但以一种不同的方式:许多分支被过度绘制：

1658510114518.png

04 聚类 and PAGA

Note：一般我们使用sc.tl.leiden，这里我们使用sc.tl.louvain是为了再现论文结果

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

running Louvain clustering
    using the "louvain" package of Traag (2017)
    finished: found 25 clusters and added
    'louvain', the cluster labels (adata.obs, categorical) (0:00:00)

接着，使用marker注释细胞群，marker如下，想要复制的可以去官网复制：

1658578907866.png

对于简单的粗粒度可视化，计算PAGA图，这是一个粗粒度的简化(抽象)图。粗粒度图中的非显著边被阈值化。

sc.tl.paga(adata, groups='louvain')
sc.pl.paga(adata, color=['louvain', 'Hba-a2', 'Elane', 'Irf8'])
pl.savefig(outdir + "03-paul15_PAGA.png")

louvain路径可视化，以及三个基因在轨迹上的可视化：

1658579586591.png

在可视化三个基因看看：

sc.pl.paga(adata, color=['louvain', 'Itga2b', 'Prss34', 'Cma1'])
pl.savefig(outdir + "03-paul15_PAGA-1.png")

1658579691897.png

实际上注释细胞簇-注意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']

sc.tl.paga(adata, groups='louvain_anno')
sc.pl.paga(adata, threshold=0.03, show=False)
pl.savefig(outdir + "03-paul15_PAGA-anno.png")

1658580217642.png

05 使用PAGA-initialization重新计算embedding

对于UMAP来说，下面的内容也是可能的：

sc.tl.draw_graph(adata, init_pos='paga')

# Now we can see all marker genes also at single-cell resolution in a meaningful layout.
sc.pl.draw_graph(adata, color=['louvain_anno', 'Itga2b', 'Prss34', 'Cma1'], legend_loc='on data')
pl.savefig(outdir + "03-paul15_PAGA-anno1.png")

结果：

1658580486265.png

修改颜色：

# Choose the colors of the clusters
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[[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


# And add some white space to some cluster names. 
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)

结果：

1658580953602.png

06 对于一组给定的基因，重组基因沿着PAGA路径变化

对于diffusion pseudotime选择根细胞

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

# 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')
pl.savefig(outdir + "04-paul15_visualization.png")

可视化结果：

1658582588244.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']

# just a cosmetic change
adata.obs['clusters'] = adata.obs['louvain_anno']  
adata.uns['clusters_colors'] = adata.uns['louvain_anno_colors']

_, 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(outdir+'./paga_path_{}.csv'.format(descr))
pl.savefig(outdir+'./paga_path_paul15.pdf')
pl.close()

三种细胞的发育轨迹：红细胞，中性粒细胞，单核细胞

1658583363488.png

这个教程没有太多的解释与说明，应该不是最初版的轨迹教程~