什么是擬時序分析？擬時序（pseudotime）分析仪搔，又稱細(xì)胞軌跡（cell trajectory）分析勃痴，通過擬時分析可以推斷出發(fā)育過程細(xì)胞的分化軌跡或細(xì)胞亞型的演化過程。
我們可以理解為在一堆細(xì)胞中包含各種各樣不同的發(fā)育狀態(tài)的細(xì)胞徙菠，有的發(fā)育早，有的發(fā)育晚郁岩，有的分化了婿奔，有的未分化，有的處于中間態(tài)问慎。利用算法基于基因表達(dá)推斷每個細(xì)胞的相對分化時間萍摊，從而確定分化軌跡。
monocle是 進(jìn)行擬時序分析常用的包如叼，這是基于R完成的冰木。但是之前也說了，monocle對于內(nèi)存消耗很大，很容易出現(xiàn)內(nèi)存不足的問題踊沸，scanpy則不會出現(xiàn)這個問題囚衔，而且scanpy內(nèi)嵌軌跡推斷函數(shù)，可以無縫銜接之前的單細(xì)胞分析雕沿。
scanpy作者使用了小鼠造血髓樣數(shù)據(jù)進(jìn)行了軌跡分析练湿，我們這兒為了方便，我們直接使用pbmc3k數(shù)據(jù)進(jìn)行測試审轮。
注：pbmc這套數(shù)據(jù)集因為本身就是基本分化完全的細(xì)胞肥哎，分化軌跡沒有啥實際生物學(xué)意義，這兒只是做測試疾渣。
單細(xì)胞轉(zhuǎn)錄數(shù)據(jù)分析之Scanpy：http://www.reibang.com/p/e22a947e6c60
單細(xì)胞轉(zhuǎn)錄組之Scanpy - 軌跡推斷/擬時序分析：http://www.reibang.com/p/0b2ca0e0b544
單細(xì)胞轉(zhuǎn)錄組之Scanpy - 樣本整合分析：http://www.reibang.com/p/beef8a8be360
單細(xì)胞空間轉(zhuǎn)錄分析之Scanpy：
單細(xì)胞空間轉(zhuǎn)錄分析之Scanpy-結(jié)合單細(xì)胞轉(zhuǎn)錄組：
導(dǎo)入相關(guān)包：
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() ##查看當(dāng)前路徑
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')
讀取數(shù)據(jù)(這兒我們使用了前面跑完scanpy流程輸出的pbmc3k.h5ad)
#讀取數(shù)據(jù),scanpy進(jìn)行聚類后的對象
data = sc.read_h5ad("./scanpy/pbmc3k.h5ad")
預(yù)處理數(shù)據(jù)篡诽，計算距離并可視化
sc.tl.draw_graph(data)
sc.pl.draw_graph(data, color='leiden', legend_loc='on data',title = "")
pt.savefig("draw_graph.pdf")
汗，作者自己認(rèn)為這兒做的這個圖類別很亂榴捡，因此作者進(jìn)行了優(yōu)化杈女，就是去噪。
對圖形進(jìn)行降噪
##去噪,擴(kuò)散圖空間表示
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")
汗吊圾，作者認(rèn)為仍舊有點(diǎn)亂（我心也亂了）
因此达椰，因此，因此
作者又提供了一種方法：Clustering and PAGA
PAGA（Partition-based Graph Abstraction）是一種基于空間劃分的抽提細(xì)胞分化“骨架”的一種算法项乒，用于顯示細(xì)胞的分化軌跡啰劲，評估cluster之間的關(guān)系緊密程度。
在這兒檀何，作者又使用sc.tl.louvain來對細(xì)胞進(jìn)行聚類蝇裤，想重現(xiàn)使用的數(shù)據(jù)的結(jié)果，好吧频鉴，我也試試louvain聚類栓辜，發(fā)現(xiàn)在pbmc中聚類結(jié)果基本一致
sc.tl.louvain(data) #可以使用resolution調(diào)節(jié)聚類的簇的數(shù)據(jù)，如resolution=1.0
sc.tl.paga(data, groups='louvain')
sc.pl.paga(data, color=['louvain', 'MS4A1', 'NKG7', 'PPBP']) ##隨便挑選了幾個基因
pt.savefig("paga_celltype.pdf")
當(dāng)然我們根據(jù)已知marker基因識別細(xì)胞類別垛孔，可以將細(xì)胞類型的信息注釋上去
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) #細(xì)胞與細(xì)胞之間的關(guān)系藕甩，距離越近表示關(guān)系越接近
pt.savefig("paga_celltype1.pdf")
#利用PAGA重新計算細(xì)胞之間的距離
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")
查看顏色，可以自行定義顏色
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")
替換顏色
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")
定義分化起點(diǎn)似炎，計算每個細(xì)胞的擬時間辛萍，畫擬時間分布（這兒我是隨便取的B細(xì)胞作為root，請根據(jù)自身數(shù)據(jù)細(xì)胞類別選认勖辍）
data.uns['iroot'] = np.flatnonzero(data.obs['louvain_anno'] == 'B')[0] ##假設(shè)分化起點(diǎn)為B cells,當(dāng)然自己分析的時候需要根據(jù)數(shù)據(jù)實際情況選擇分化起點(diǎn)
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路徑重建基因變化
adata =sc.read(results_file)
gene_names = ['IL7R', 'CD14', 'LYZ', 'MS4A1', 'CD8A', 'CD8B', 'FCGR3A', 'MS4A7','GNLY', 'NKG7', 'FCER1A', 'CST3'] # 選取了一列marker 基因悯许，要根據(jù)實際情況選取
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) # 去除表達(dá)基因200以下的細(xì)胞
sc.pp.filter_genes(data_raw, min_cells=3) # 去除在3個細(xì)胞以下表達(dá)的基因
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 #使用完整的原始數(shù)據(jù)進(jìn)行可視化
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()
