今天給大家分享一個用于擬時分析的軟件Palantir，當然，目前擬時分析的軟件已經(jīng)有很多了，至于哪個適合于你，就要看你的目的和結(jié)果是否符合你的預(yù)期了。
Palantir是2020年3月份發(fā)表于nature biotechnology，文章名稱是Characterization of cell fate probabilities in single-cell data with Palantir,大家可以看一下，這里我們來分享一下這個軟件的用法。

介紹

Palantir是一種沿分化軌跡對齊細胞的算法。 Palantir將分化建模為一個隨機過程，在該過程中，干細胞通過低維表型歧管通過一系列步驟分化為最終分化的細胞。 Palantir(是一個python模塊)有效地捕獲了細胞狀態(tài)的連續(xù)性和細胞命運決定的隨機性。
首先

import palantir
import scanpy as sc
import numpy as np
import os
# Plotting 
import matplotlib
import matplotlib.pyplot as plt
# Reset random seed
np.random.seed(5)

第二步，載入數(shù)據(jù)

###Loading data
ad = sc.read(dir + 'marrow_sample_scseq_counts.h5ad')
####這里是根據(jù)scanpy讀取的h5ad文件，如果是矩陣文件，如下用法：
adata=sc.read_csv（dir + 'marrow_sample_scseq_counts.csv'）
####可不可以直接用python讀取R語言的rds或者Rdata呢？期待大家開發(fā)

第三步：
Data processing

sc.pp.normalize_per_cell(ad)    ####這個就是scanpy的均一化步驟，里面有參數(shù)可以調(diào)整。

我們建議對數(shù)據(jù)進行l(wèi)og轉(zhuǎn)換（對于單細胞這種稀疏矩陣，大范圍符合負二項分布的數(shù)據(jù)最好進行l(wèi)og轉(zhuǎn)化）。 請注意，一些數(shù)據(jù)集在線性尺度上顯示出更好的信號（非稀疏矩陣更合適），而其他數(shù)據(jù)集在對數(shù)尺度上顯示出更強的信號。這也就是為什么單細胞數(shù)據(jù)矩陣均一化以后為什么還要log的原因。
下面的函數(shù)使用0.1的偽計數(shù)代替1。
第四步：

palantir.preprocess.log_transform(ad)###也就是為了增強信號

Highly variable gene selection

sc.pp.highly_variable_genes(ad, n_top_genes=1500, flavor='cell_ranger')
###這個要和根據(jù)自己的數(shù)據(jù)類型進行修改,參數(shù)的含義如下
flavor
        Choose the flavor for computing normalized dispersion. In their default
        workflows, Seurat passes the cutoffs whereas Cell Ranger passes
        `n_top_genes`。

PCA
PCA是Palantir數(shù)據(jù)處理的第一步。 為了克服在單細胞RNA-seq數(shù)據(jù)中普遍存在的廣泛缺失，必須采用這種表示法。
我們建議不要使用固定數(shù)量的PC，而應(yīng)使用能夠解釋高度可變的基因選擇后數(shù)據(jù)中85％的方差的主成分。

pca_projections, _ = palantir.utils.run_pca(ad, use_hvg=False)

Diffusion maps
這里稍微解釋一下，在對細胞進行軌跡分析的過程中，一般在挑選擬時分析所用的基因過程中，有三種方式：
1、采用類似于Seurat的方式，選擇高變基因，但是高變基因方差較大，在細胞中表達大部分出現(xiàn)斷層，不利于擬時過程的連續(xù)性（但也有效果好的，比如沒有分化關(guān)系的用這個基因做擬時效果好一點）。
2、采用類似于PAGA的方式進行擬時分析，就是現(xiàn)在用到的Diffusion maps，會檢測單細胞數(shù)據(jù)中連續(xù)變化的基因，這樣做擬時分析更有連續(xù)性，但是對于沒有分化關(guān)系的細胞類型，這個方法不太適用。
3、人工挑選基因，上面選擇的基因是否跟發(fā)育等研究目的相關(guān)，不太清楚，如果有很深的生物學(xué)背景，人工選擇更好一點。

# Run diffusion maps
dm_res = palantir.utils.run_diffusion_maps(pca_projections, n_components=5)

使用以下函數(shù)基于特征間隙估計數(shù)據(jù)的低維嵌入

ms_data = palantir.utils.determine_multiscale_space(dm_res)

如果在上述步驟中手動指定特征向量的數(shù)量，請確保指定的參數(shù)> 2

接下來的就是可視化

在這個例子中，我們采用tSNE投影，using diffusion components to visualize the data?，F(xiàn)在，我們建議使用力導(dǎo)向布局來可視化軌跡。 力導(dǎo)向的布局可以由用于確定擴散圖的相同自適應(yīng)內(nèi)核來計算。 仍然可以使用函數(shù)tsne=palantir.utils.run_tsne（ms_data）計算擴散成分上的tSNE
（當然這里我們也可以使用UMAP ）
我們的軟件包Harmony提供了用于計算力導(dǎo)向布局的界面(大家可以看一下)。

import harmony
fdl = harmony.plot.force_directed_layout(dm_res['kernel'], ad.obs_names)
fig, ax = palantir.plot.plot_tsne(fdl)

11.png

為了與先前的教程和手稿保持一致，我們將為此數(shù)據(jù)集使用預(yù)先計算的tSNE。

import pandas as pd
tsne = pd.read_pickle(palantir_dir + 'data/sample_tsne.p')
fig, ax = palantir.plot.plot_tsne(tsne)

23.png

MAGIC 歸因

MAGIC是Peerer實驗室開發(fā)的一種插補技術(shù)，用于單細胞數(shù)據(jù)插補。 Palantir使用MAGIC估算數(shù)據(jù)以可視化并確定基因表達趨勢。

imp_df = palantir.utils.run_magic_imputation(ad, dm_res)

使用plot_gene_expression函數(shù)，可以在tSNE圖上顯示基因表達。 genes參數(shù)是字符串的可迭代的基因，它們是列名表達的子集。 下面的函數(shù)繪制了HSC基因CD34，骨髓基因MPO和類紅細胞前體基因GATA1和樹突狀細胞基因IRF8的表達。

palantir.plot.plot_gene_expression(imp_df, tsne, ['CD34', 'MPO', 'GATA1', 'IRF8'])

圖片.png

The same functions can be used to plot gene expression on force directed layouts.

palantir.plot.plot_gene_expression(imp_df, fdl, ['CD34', 'MPO', 'GATA1', 'IRF8'])

Diffusion maps visualization

The computed diffusion components can be visualized with the following snippet.

palantir.plot.plot_diffusion_components(tsne, dm_res)

圖片.png

Running Palantir

可以通過指定一個近似的早期細胞來運行Palantir。
Palantir也可以自動確定終端狀態(tài)。 在此數(shù)據(jù)集中，我們知道了終端狀態(tài)，我們將使用terminal_states參數(shù)進行設(shè)置

terminal_states = pd.Series(['DC', 'Mono', 'Ery'], 
                           index=['Run5_131097901611291', 'Run5_134936662236454', 'Run4_200562869397916'])
start_cell = 'Run5_164698952452459'
pr_res = palantir.core.run_palantir(ms_data, start_cell, num_waypoints=500, terminal_states=terminal_states.index)

也就是說，最好有一定的先驗知識。
Palantir generates the following results
1、Pseudotime: Pseudo time ordering of each cell
2、Terminal state probabilities: Matrix of cells X terminal states. Each entry represents the probability of the corresponding cell reaching the respective terminal state
3、Entropy: A quantiative measure of the differentiation potential of each cell computed as the entropy of the multinomial terminal state probabilities（有點意思）
The terminal states in this dataset are renamed to reflect the known biology below

pr_res.branch_probs.columns = terminal_states[pr_res.branch_probs.columns]

Visualizing Palantir results

Palantir results can be visualized on the tSNE map using the plot_palantir_results function

palantir.plot.plot_palantir_results(pr_res, tsne)

可以使用plot_terminal_state_probs函數(shù)可視化各個單元格的終端狀態(tài)概率分布

cells = ['Run5_164698952452459', 'Run5_170327461775790', 'Run4_121896095574750', ]
palantir.plot.plot_terminal_state_probs(pr_res, cells) 

基因表達趨勢

Palantir使用通用加性模型（GAM）確定沿不同譜系的基因表達趨勢。 可以使用以下代碼段確定標記趨勢。 這將計算所有譜系的趨勢。 可以使用lineages參數(shù)使用血統(tǒng)的子集。

genes = ['CD34', 'MPO', 'GATA1', 'IRF8']
gene_trends = palantir.presults.compute_gene_trends( pr_res, imp_df.loc[:, genes])
palantir.plot.plot_gene_trends(gene_trends)

圖片.png

palantir.plot.plot_gene_trend_heatmaps(gene_trends)

圖片.png

Clustering

這個聚類跟我們常用的Seurat不是一回事
Cells can be clustered and visualized using Phenograph with the following snippet.

clusters = palantir.utils.determine_cell_clusters(pca_projections)

圖片.png

類似于monocle的state。
當然，這里我們可以人為指定cluster，與Seurat進行對接

palantir.plot.plot_gene_trend_clusters(trends, gene_clusters)

圖片.png

這個方法雖然用的不多，但是算法還是很好的，大家可以多借鑒借鑒。