- 學(xué)習(xí)資料來(lái)源：
- scanpy主頁(yè)：https://scanpy.readthedocs.io/en/stable/
- 官網(wǎng)：https://scanpy-tutorials.readthedocs.io/en/latest/paga-paul15.html【注意教程有兩個(gè)版本播玖，這里是latest版本的學(xué)習(xí)筆記】
數(shù)據(jù)說(shuō)明：
使用來(lái)自的Paul et al. (2015)數(shù)據(jù)重建髓系myeloid和紅細(xì)胞系erythroid的分化
本次會(huì)用到一個(gè)新的軟件枣抱，安裝如下：
conda activate scanpy
conda install -c anaconda pandas -y
01 數(shù)據(jù)導(dǎo)入
首先去：將數(shù)據(jù)下載下來(lái)剩蟀，這里直接封裝到了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')
# 保存數(shù)據(jù)路徑
outdir = './Trajectory/'
adata = sc.datasets.paul15()
adata
# 讓我們用比默認(rèn)的float32更高的精度來(lái)確保在不同的計(jì)算平臺(tái)上得到完全相同的結(jié)果
# this is not required and results will be comparable without it
adata.X = adata.X.astype('float64')
adata
總共有：2730個(gè)細(xì)胞陨界，3451個(gè)基因
AnnData object with n_obs × n_vars = 2730 × 3451
obs: 'paul15_clusters'
uns: 'iroot'
02 預(yù)處理和可視化
這里用到了scanpy.pp.recipe_zheng17這個(gè)函數(shù)诞仓，主要是將數(shù)據(jù)預(yù)處理的幾個(gè)步驟包裝成一個(gè)函數(shù)，處理方式來(lái)自文章：
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)
進(jìn)行預(yù)處理：
# 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 去噪：可選部分
為了去除圖的噪聲，我們?cè)跀U(kuò)散映射空間(而不是 PCA 空間)中表示它衔彻。計(jì)算幾個(gè)擴(kuò)散成分（diffusion components）內(nèi)的距離相當(dāng)于去噪圖形-我們只取幾個(gè)第一個(gè)光譜成分（the first spectral components）薇宠。這與使用PCA去噪數(shù)據(jù)矩陣非常相似。
該方法已在幾篇論文中使用艰额，如Schiebinger et al. (2017) or Tabaka et al. (2018)澄港。這也與MAGIC Dijk et al. (2018)等人背后的原則有關(guān)。
Note：這不是一個(gè)必要的步驟柄沮，如 PAGA回梧、聚類(lèi)、偽時(shí)間估計(jì)等分析都不是一個(gè)必要的步驟祖搓。你最好使用一個(gè)無(wú)噪聲的圖形狱意。在許多情況下(也在這里) ，這將給你帶來(lái)非常體面的結(jié)果棕硫。
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()
這看起來(lái)仍然很混亂髓涯，但以一種不同的方式:許多分支被過(guò)度繪制：
1658510114518.png
04 聚類(lèi) and PAGA
Note：一般我們使用sc.tl.leiden，這里我們使用sc.tl.louvain是為了再現(xiàn)論文結(jié)果
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注釋細(xì)胞群，marker如下蚓再，想要復(fù)制的可以去官網(wǎng)復(fù)制：
1658578907866.png
對(duì)于簡(jiǎn)單的粗粒度可視化滑肉，計(jì)算PAGA圖，這是一個(gè)粗粒度的簡(jiǎn)化(抽象)圖摘仅。粗粒度圖中的非顯著邊被閾值化靶庙。
sc.tl.paga(adata, groups='louvain')
sc.pl.paga(adata, color=['louvain', 'Hba-a2', 'Elane', 'Irf8'])
pl.savefig(outdir + "03-paul15_PAGA.png")
louvain路徑可視化，以及三個(gè)基因在軌跡上的可視化：
1658579586591.png
在可視化三個(gè)基因看看：
sc.pl.paga(adata, color=['louvain', 'Itga2b', 'Prss34', 'Cma1'])
pl.savefig(outdir + "03-paul15_PAGA-1.png")
1658579691897.png
實(shí)際上注釋細(xì)胞簇-注意Cma1是肥大細(xì)胞標(biāo)記娃属，只出現(xiàn)在祖細(xì)胞/干細(xì)胞簇8中的一小部分細(xì)胞中六荒，見(jiàn)下面的單細(xì)胞分辨圖：
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')
注釋?zhuān)?/p>
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重新計(jì)算embedding
對(duì)于UMAP來(lái)說(shuō)，下面的內(nèi)容也是可能的：
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")
結(jié)果：
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)
結(jié)果：
1658580953602.png
06 對(duì)于一組給定的基因矾端，重組基因沿著PAGA路徑變化
對(duì)于diffusion pseudotime選擇根細(xì)胞
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")
可視化結(jié)果：
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()
三種細(xì)胞的發(fā)育軌跡：紅細(xì)胞掏击，中性粒細(xì)胞，單核細(xì)胞
1658583363488.png
這個(gè)教程沒(méi)有太多的解釋與說(shuō)明秩铆，應(yīng)該不是最初版的軌跡教程~
