什么是拟时序分析？拟时序（pseudotime）分析，又称细胞轨迹（cell trajectory）分析，通过拟时分析可以推断出发育过程细胞的分化轨迹或细胞亚型的演化过程。
我们可以理解为在一堆细胞中包含各种各样不同的发育状态的细胞，有的发育早，有的发育晚，有的分化了，有的未分化，有的处于中间态。利用算法基于基因表达推断每个细胞的相对分化时间，从而确定分化轨迹。
monocle是 进行拟时序分析常用的包，这是基于R完成的。但是之前也说了，monocle对于内存消耗很大，很容易出现内存不足的问题，scanpy则不会出现这个问题，而且scanpy内嵌轨迹推断函数，可以无缝衔接之前的单细胞分析。
scanpy作者使用了小鼠造血髓样数据进行了轨迹分析，我们这儿为了方便，我们直接使用pbmc3k数据进行测试。
注：pbmc这套数据集因为本身就是基本分化完全的细胞，分化轨迹没有啥实际生物学意义，这儿只是做测试。
单细胞转录数据分析之Scanpy：https://www.jianshu.com/p/e22a947e6c60
单细胞转录组之Scanpy - 轨迹推断/拟时序分析：https://www.jianshu.com/p/0b2ca0e0b544
单细胞转录组之Scanpy - 样本整合分析：https://www.jianshu.com/p/beef8a8be360
单细胞空间转录分析之Scanpy：
单细胞空间转录分析之Scanpy-结合单细胞转录组：
Scanpy trajectory inference
导入相关包：
import numpy as np
import pandas as pd
import matplotlib.pyplot as pt
from matplotlib import rcParams
import scanpy as sc
sc.settings.verbosity = 3
sc.logging.print_versions()
import os
os.getcwd() ##查看当前路径
os.chdir(./scanpy/pseudo') ##修改路径
os.getcwd()
results_file = 'pbmc3k_pseudo.h5ad'
sc.settings.set_figure_params(dpi=80, frameon=False, figsize=(3, 3), facecolor='white')
读取数据(这儿我们使用了前面跑完scanpy流程输出的pbmc3k.h5ad)
#读取数据,scanpy进行聚类后的对象
data = sc.read_h5ad("./scanpy/pbmc3k.h5ad")
预处理数据，计算距离并可视化
sc.tl.draw_graph(data)
sc.pl.draw_graph(data, color='leiden', legend_loc='on data',title = "")
pt.savefig("draw_graph.pdf")
draw_graph
汗，作者自己认为这儿做的这个图类别很乱，因此作者进行了优化，就是去噪。
对图形进行降噪
##去噪,扩散图空间表示
sc.tl.diffmap(data)
sc.pp.neighbors(data, n_neighbors=10, use_rep='X_diffmap')
sc.tl.draw_graph(data)
sc.pl.draw_graph(data, color='leiden', legend_loc='on data') #legend_loc='right margin'
pt.savefig("draw_graph_diffmap.pdf")
draw_graph +去噪
汗，作者认为仍旧有点乱（我心也乱了）
因此，因此，因此
作者又提供了一种方法：Clustering and PAGA
PAGA（Partition-based Graph Abstraction）是一种基于空间划分的抽提细胞分化“骨架”的一种算法，用于显示细胞的分化轨迹，评估cluster之间的关系紧密程度。
在这儿，作者又使用sc.tl.louvain来对细胞进行聚类，想重现使用的数据的结果，好吧，我也试试louvain聚类，发现在pbmc中聚类结果基本一致
sc.tl.louvain(data) #可以使用resolution调节聚类的簇的数据，如resolution=1.0
sc.tl.paga(data, groups='louvain')
sc.pl.paga(data, color=['louvain', 'MS4A1', 'NKG7', 'PPBP']) ##随便挑选了几个基因
pt.savefig("paga_celltype.pdf")
papg
当然我们根据已知marker基因识别细胞类别，可以将细胞类型的信息注释上去
data.obs['louvain'].cat.categories
data.obs['louvain_anno'] = data.obs['louvain']
data.obs['louvain_anno'].cat.categories = ['CD4 T', 'CD14 Monocytes','B', 'CD8 T','NK', 'FCGR3A Monocytes','Dendritic', 'Megakaryocytes']
sc.tl.paga(data, groups='louvain_anno')
sc.pl.paga(data, threshold=0.03, show=False) #细胞与细胞之间的关系，距离越近表示关系越接近
pt.savefig("paga_celltype1.pdf")
paga_celltype
#利用PAGA重新计算细胞之间的距离
sc.tl.draw_graph(data, init_pos='paga')
sc.pl.draw_graph(data, color=['louvain_anno','MS4A1', 'NKG7', 'PPBP'], legend_loc='on data')
pt.savefig("paga_celltype_graph.pdf")
paga_celltype
查看颜色，可以自行定义颜色
pt.figure(figsize=(8, 2))
for i in range(28):
pt.scatter(i, 1, c=sc.pl.palettes.zeileis_28[i], s=200)
pt.show()
pt.savefig("colour.pdf")
colour
替换颜色
zeileis_colors = np.array(sc.pl.palettes.zeileis_28)
new_colors = np.array(data.uns['louvain_anno_colors'])
new_colors[[0]] = zeileis_colors[[12]] # CD4 T colors / green
new_colors[[1]] = zeileis_colors[[5]] # CD14 Monocytes colors / red
new_colors[[2]] = zeileis_colors[[17]] # B colors / yellow
new_colors[[3]] = zeileis_colors[[2]] # CD8 T / grey
new_colors[[4]] = zeileis_colors[[18]] # NK / turquoise
new_colors[[5]] = zeileis_colors[[6]] # FCGR3A Monocytes / light blue
new_colors[[6]] = zeileis_colors[[0]] # Dendritic / dark blue
new_colors[[7]] = zeileis_colors[[25]] # Megakaryocytes / grey
#new_colors[[10, 17, 5, 3, 15, 6, 18, 13, 7, 12]] = zeileis_colors[[5, 5, 5, 5, 11, 11, 10, 9, 21, 21]] # CD14 Monocytes colors / red
data.uns['louvain_anno_colors'] = new_colors
sc.pl.paga_compare(
data, threshold=0.03, title='', right_margin=0.2, size=10, edge_width_scale=0.5,
legend_fontsize=12, fontsize=12, frameon=False, edges=True)
pt.savefig("paga_compare.pdf")
paga_compare
定义分化起点，计算每个细胞的拟时间，画拟时间分布（这儿我是随便取的B细胞作为root，请根据自身数据细胞类别选取）
data.uns['iroot'] = np.flatnonzero(data.obs['louvain_anno'] == 'B')[0] ##假设分化起点为B cells,当然自己分析的时候需要根据数据实际情况选择分化起点
sc.tl.dpt(data)
sc.pl.draw_graph(data, color=['louvain_anno', 'dpt_pseudotime'], legend_loc='on data',title = ['','pseudotime'], frameon=True)
pt.savefig("paga_peudotime.pdf")
sc.pl.draw_graph(data, color=['louvain', 'dpt_pseudotime'], legend_loc='on data',title = ['','pseudotime'], frameon=True)
pt.savefig("paga_peudotime1.pdf")
data.write(results_file)
paga_peudotime
paga_peudotime
针对给定的一组基因，沿PAGA路径重建基因变化
adata =sc.read(results_file)
gene_names = ['IL7R', 'CD14', 'LYZ', 'MS4A1', 'CD8A', 'CD8B', 'FCGR3A', 'MS4A7','GNLY', 'NKG7', 'FCER1A', 'CST3'] # 选取了一列marker 基因，要根据实际情况选取
data_raw = sc.read_10x_mtx('./filtered_gene_bc_matrices/hg19/', var_names='gene_symbols', cache=True)
data_raw.var_names_make_unique()
sc.pp.filter_cells(data_raw, min_genes=200) # 去除表达基因200以下的细胞
sc.pp.filter_genes(data_raw, min_cells=3) # 去除在3个细胞以下表达的基因
mito_genes=data_raw.var_names.str.startswith('MT-')
data_raw.obs['percent_mito']=np.sum(data_raw[:,mito_genes].X,axis=1).A1/np.sum(data_raw.X,axis=1).A1
data_raw.obs['n_counts']=data_raw.X.sum(axis=1).A1
data_raw = data_raw[data_raw.obs.n_genes < 2500, :]
data_raw = data_raw[data_raw.obs.percent_mito < 0.05, :]
sc.pp.normalize_total(data_raw, target_sum=1e4)
sc.pp.log1p(data_raw)
sc.pp.scale(data_raw)
adata.raw = data_raw #使用完整的原始数据进行可视化
paths = [('NK', [2,0,3]), ('FCGR3A Monocytes', [2,6,1])]
adata.obs['distance'] = adata.obs['dpt_pseudotime']
adata.obs['clusters'] = adata.obs['louvain'] # just a cosmetic change
adata.uns['clusters_colors'] = adata.uns['louvain_anno_colors']
_, axs = pt.subplots(ncols=2, figsize=(6, 2.5), gridspec_kw={'wspace': 0.05, 'left': 0.12})
pt.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('./paga_path_{}.csv'.format(descr))
pt.savefig('paga_path_pbmc.pdf')
pt.show()
paga_path_pbmc
