拟时分析软件Palantir

今天给大家分享一个用于拟时分析的软件Palantir，当然，目前拟时分析的软件已经有很多了，至于哪个适合于你，就要看你的目的和结果是否符合你的预期了。
Palantir是2020年3月份发表于nature biotechnology，文章名称是Characterization of cell fate probabilities in single-cell data with Palantir,大家可以看一下，这里我们来分享一下这个软件的用法。

介绍

Palantir是一种沿分化轨迹对齐细胞的算法。 Palantir将分化建模为一个随机过程，在该过程中，干细胞通过低维表型歧管通过一系列步骤分化为最终分化的细胞。 Palantir(是一个python模块)有效地捕获了细胞状态的连续性和细胞命运决定的随机性。
首先

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)

第二步，载入数据

###Loading data
ad = sc.read(dir + 'marrow_sample_scseq_counts.h5ad')
####这里是根据scanpy读取的h5ad文件，如果是矩阵文件，如下用法：
adata=sc.read_csv（dir + 'marrow_sample_scseq_counts.csv'）
####可不可以直接用python读取R语言的rds或者Rdata呢？期待大家开发

第三步：
Data processing

sc.pp.normalize_per_cell(ad)    ####这个就是scanpy的均一化步骤，里面有参数可以调整。

我们建议对数据进行log转换（对于单细胞这种稀疏矩阵，大范围符合负二项分布的数据最好进行log转化）。 请注意，一些数据集在线性尺度上显示出更好的信号（非稀疏矩阵更合适），而其他数据集在对数尺度上显示出更强的信号。这也就是为什么单细胞数据矩阵均一化以后为什么还要log的原因。
下面的函数使用0.1的伪计数代替1。
第四步：

palantir.preprocess.log_transform(ad)###也就是为了增强信号

Highly variable gene selection

sc.pp.highly_variable_genes(ad, n_top_genes=1500, flavor='cell_ranger')
###这个要和根据自己的数据类型进行修改,参数的含义如下
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数据处理的第一步。 为了克服在单细胞RNA-seq数据中普遍存在的广泛缺失，必须采用这种表示法。
我们建议不要使用固定数量的PC，而应使用能够解释高度可变的基因选择后数据中85％的方差的主成分。

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

Diffusion maps
这里稍微解释一下，在对细胞进行轨迹分析的过程中，一般在挑选拟时分析所用的基因过程中，有三种方式：
1、采用类似于Seurat的方式，选择高变基因，但是高变基因方差较大，在细胞中表达大部分出现断层，不利于拟时过程的连续性（但也有效果好的，比如没有分化关系的用这个基因做拟时效果好一点）。
2、采用类似于PAGA的方式进行拟时分析，就是现在用到的Diffusion maps，会检测单细胞数据中连续变化的基因，这样做拟时分析更有连续性，但是对于没有分化关系的细胞类型，这个方法不太适用。
3、人工挑选基因，上面选择的基因是否跟发育等研究目的相关，不太清楚，如果有很深的生物学背景，人工选择更好一点。

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

使用以下函数基于特征间隙估计数据的低维嵌入

ms_data = palantir.utils.determine_multiscale_space(dm_res)

如果在上述步骤中手动指定特征向量的数量，请确保指定的参数> 2

接下来的就是可视化

在这个例子中，我们采用tSNE投影，using diffusion components to visualize the data。现在，我们建议使用力导向布局来可视化轨迹。 力导向的布局可以由用于确定扩散图的相同自适应内核来计算。 仍然可以使用函数tsne=palantir.utils.run_tsne（ms_data）计算扩散成分上的tSNE
（当然这里我们也可以使用UMAP ）
我们的软件包Harmony提供了用于计算力导向布局的界面(大家可以看一下)。

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

11.png

为了与先前的教程和手稿保持一致，我们将为此数据集使用预先计算的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实验室开发的一种插补技术，用于单细胞数据插补。 Palantir使用MAGIC估算数据以可视化并确定基因表达趋势。

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

使用plot_gene_expression函数，可以在tSNE图上显示基因表达。 genes参数是字符串的可迭代的基因，它们是列名表达的子集。 下面的函数绘制了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也可以自动确定终端状态。 在此数据集中，我们知道了终端状态，我们将使用terminal_states参数进行设置

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函数可视化各个单元格的终端状态概率分布

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

基因表达趋势

Palantir使用通用加性模型（GAM）确定沿不同谱系的基因表达趋势。 可以使用以下代码段确定标记趋势。 这将计算所有谱系的趋势。 可以使用lineages参数使用血统的子集。

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

这个方法虽然用的不多，但是算法还是很好的，大家可以多借鉴借鉴。