與上次翻車實驗的不同：

1. 過濾了更多的細(xì)胞和基因。在處理數(shù)據(jù)時，低質(zhì)量的細(xì)胞一定要清除掉，告誡大家寧缺毋濫。。。否則后續(xù)分析真的是一堆的噪點。
2. 跳過了scanpy中的降噪步驟。scanpy的降噪算法做的真不咋地，一個好圖降噪之后更加混亂了。這個故事告訴我們一個道理：不要迷信作者說的話！??！

對翻車實驗感興趣的話：http://www.itdecent.cn/p/e688646a451b

什么是軌跡推斷/偽時間分析？

軌跡推斷對應(yīng)英文是trajectory inference，偽時間分析對應(yīng)英文是pseudotime analysis。其實它們做的是同一件事。

在單細(xì)胞轉(zhuǎn)錄組測序完成的時候，無論是細(xì)胞在發(fā)育、分化還是凋亡的過程中，它們的生命活動已經(jīng)在這一時刻靜止了，因此我們可以將這些測序數(shù)據(jù)看作是這一個時刻細(xì)胞呈現(xiàn)的表達(dá)狀態(tài)的瞬時截圖。

而軌跡推斷/偽時間分析則構(gòu)建了一個細(xì)胞發(fā)育和分化的模型。這一個瞬時截圖內(nèi)的細(xì)胞在各種各樣不同的發(fā)育狀態(tài)，有的發(fā)育早，有的發(fā)育晚，有的分化了，有的未分化，甚至有的處于中間態(tài)。而決定細(xì)胞分化往往就是基因表達(dá)的變化。

軌跡推斷利用各類細(xì)胞基因表達(dá)的連續(xù)性建立了一個“分化軌跡”，再將這些細(xì)胞投射到這個“分化軌跡”上。如此一來，這個靜態(tài)的截圖就呈現(xiàn)出了動態(tài)的過程。

PS：細(xì)胞類型分類真是一場艱苦卓絕的戰(zhàn)斗?。?！一入文獻(xiàn)深似海?。。。?/p>

0. 數(shù)據(jù)預(yù)處理

先加載需要的包：

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 anndata==0.6.19 numpy==1.14.5 scipy==1.1.0 pandas==0.23.4 scikit-learn==0.19.2 statsmodels==0.9.0 python-igraph==0.7.1 louvain==0.6.1 

讀取數(shù)據(jù)（數(shù)據(jù)來源詳見http://www.itdecent.cn/p/b190efae4d31）
這個數(shù)據(jù)是利用scanpy進(jìn)行聚類后的對象。

adata = sc.read_h5ad("/data1/zll/project/deepBase3/HCA/bone_marrow/scanpy/3_29_PC16_filterMore/umap_tsne_3_29.h5ad")

預(yù)處理數(shù)據(jù)，計算距離并可視化。

sc.tl.draw_graph(adata)
sc.pl.draw_graph(adata, color='louvain', legend_loc='on data',title = "")

output_4_0.png

1. Partition-based Graph Abstraction（PAGA）

這是一種基于空間劃分的抽提細(xì)胞分化“骨架”的一種算法，用于顯示細(xì)胞的分化軌跡，得到哪些cluster之間的關(guān)系更接近。

sc.tl.paga(adata, groups='louvain')
sc.pl.paga(adata, color=['louvain'],title = "")

output_13_1.png

AVP是造血干細(xì)胞表達(dá)的一個基因，給這個基因上色，我們可以看到基本上都富集在了cluster13的位置。那么有很大可能這個cluster就是造血干細(xì)胞。

sc.pl.paga(adata, color=['AVP'])

output_15_1.png

2.注釋細(xì)胞類型

上面的每個cluster都是用數(shù)字編號代替，為了能更清楚地知道這些細(xì)胞之間的關(guān)系，將細(xì)胞類型的信息注釋上去。細(xì)胞類型的確定主要利用了marker基因以及現(xiàn)有的文獻(xiàn)和資料，是一個非常復(fù)雜的過程?？傊谂c文獻(xiàn)的一番斗爭之后，找到了以下這些骨髓的細(xì)胞類型（僅供參考）：

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

adata.obs['louvain'].cat.categories

返回共有24個cluster：

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

將編號所對應(yīng)的細(xì)胞類型注釋上去：

adata.obs['louvain_anno'] = adata.obs['louvain']

# annote them with names
adata.obs['louvain_anno'].cat.categories = ['0/nT', '1/B', '2nT', '3/CD8+T',
                                            '4/Neu', '5/nT', '6/Neu', '7/NK',
                                            '8/Eos', '9/CD8+T', '10/Eos', '11/pre-B',
                                            '12/NRBC','13/HSC', '14/Mon', '15/plasma', 
                                            '16/preDC', '17/DC', '18/plt', '19/pro-B', 
                                            '20/strom','21/Eos','22/','23/pro-B']

3. 計算 PAGA并可視化

sc.tl.paga(adata, groups='louvain_anno')
sc.pl.paga(adata, threshold=0.03)

output_23_1.png

于是就得到了這樣像“骨架”一樣的結(jié)果。這個圖顯示的就是細(xì)胞與細(xì)胞之間的關(guān)系，距離越近表示關(guān)系越接近。

4. 利用PAGA重新計算細(xì)胞之間的距離

還記得我們第0步計算的距離嗎？現(xiàn)在我們要將細(xì)胞在PAGA這個“骨架”上重現(xiàn)出來，就利用PAGA的計算結(jié)果，把細(xì)胞放上去。

sc.tl.draw_graph(adata, init_pos='paga')
sc.pl.draw_graph(adata, color=['louvain_anno'], legend_loc='right margin')

output_27_0.png

5.（可選）美化圖片

Choose the colors of the clusters a bit more consistently.

pl.figure(figsize=(8, 2))
for i in range(28):
    pl.scatter(i, 1, c=sc.pl.palettes.zeileis_26[i], s=200)
pl.show()

output_30_0.png

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

new_colors[[13]] = zeileis_colors[[17]]  # Stem colors / yellow
new_colors[[12]] = zeileis_colors[[5]]  # Ery colors / red
new_colors[[4,6]] = zeileis_colors[[12,12]] # Neutrophil / green
new_colors[[8,10]] = zeileis_colors[[11,11]] # Eosinophil / light red
new_colors[[14]] = zeileis_colors[[19]] # monocyte / light green
new_colors[[16,17]] = zeileis_colors[[18,26]]  # DC / yellow
new_colors[[0,2,3,5,9]] = zeileis_colors[[7,7,6,7,6]]  # naive T & CD8+T / light blue
new_colors[[7]] = zeileis_colors[[0]]  # NK cell / dark blue
new_colors[[1,11,19,23]] = zeileis_colors[[23,22,9,9]]  # B cell / pink

new_colors[[22]] = zeileis_colors[[25]]  # Not known / grey
new_colors[[15]] = zeileis_colors[[3]] #plasma cell / blue grey
new_colors[[18]] = zeileis_colors[[4]] #platelet / pink 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)

--> added 'pos', the PAGA positions (adata.uns['paga'])
saving figure to file ./figures/paga_compare.pdf

output_35_1.png

6. 計算偽時間值

首先需要定義一個偽時間值為0的細(xì)胞群體，一般來說就是干細(xì)胞或者祖細(xì)胞?？傊褪亲钤嫉募?xì)胞類型。

# the most primitive cell is refered as 0 persudotime.
# Group 13 is the nearest cell population to Hematopoietic stem cell.
adata.uns['iroot'] = np.flatnonzero(adata.obs['louvain_anno']  == '13/HSC')[0]
sc.tl.dpt(adata)

sc.pl.draw_graph(adata, color=['louvain_anno', 'dpt_pseudotime'],
                 legend_loc='right margin',title = ['','pseudotime'])

output_38_0.png

偽時間圖的分化軌跡無法用色彩體現(xiàn)出來，可能是偽時間參數(shù)的問題。

上面的圖注被遮住了，所以單獨做一個圖。

#plot again to see full legends info
sc.pl.draw_graph(adata, color=['louvain_anno'],
                 legend_loc='right margin',title = ['']) 

output_39_0.png

保存結(jié)果文件：

adata.write("/data1/zll/project/deepBase3/HCA/bone_marrow/scanpy/3_31_traj/trajectory_3_31.h5ad")