シングルセルデータからPalantirで細胞運命の解析をやってみた

はじめに

近年、単一細胞解析技術の進展により、細胞の分化過程や運命決定のメカニズムを解明する研究が加速しています。しかしながら、このような高次元で複雑なデータを解析するには、適切なアルゴリズムが必要不可欠です。そこで注目されているのが、Palantirというツールです。

Palantirでは、幹細胞から最終的な細胞タイプに至る道筋を実験室で実際に観察する代わりに、計算上の「疑似時間軸」として作り出します。そしてその疑似時間に基づき細胞の運命確率を計算します。ここで言う「時間」は単なる経過時間ではなく、「どの運命に向かっているか」という確率の変化を表す独自の物差しです。例えば、ある血液幹細胞が赤血球になる確率が60%、白血球になる確率が30%...といった具合に、細胞の未来を数値化します。
本記事では、Palantirの基本的な仕組みや利点、具体的な活用方法について解説します。

分化軌跡解析ツールのPalantirの説明

以下の論文に基づきPalantirの仕組みを簡単に説明します。

Setty, M., Kiseliovas, V., Levine, J. et al. Characterization of cell fate probabilities in single-cell data with Palantir. Nat Biotechnol 37, 451–460 (2019). 
https://doi.org/10.1038/s41587-019-0068-4

Palantirの原理とは、高解像度の疑似時間順序付けの細胞を生成し、各細胞状態に対して各最終状態に分化する確率を割り当てることです。これは従来の細胞運命モデルと異なります。
従来の細胞運命モデルでは、細胞分化は一連の離散的なステップを経て進行し、各段階で細胞が特定の運命に固定されると考えられていました。この枠組みでは、細胞分化は直線的であり、分岐点ごとに明確な決定が行われるという前提がありました。
しかし、単一細胞RNAシーケンス（scRNA-seq）の技術が進歩するにつれ、細胞分化や運命決定が実際には連続的かつ確率的である可能性が示唆されるようになりました。

大量の細胞を見るとそこには様々な分化の程度の細胞があり、そこから各細胞に対して疑似的な時間を生成します。そしてそうした各細胞に対して最終分化状態になる確率を割り当てることでどの程度他の運命になりうるのか（ならないのか）を定量化できるということです。

環境

実行環境は以下の通りです。

OS - 5.4.0-67-generic #75-Ubuntu SMP Fri Feb 19 18:03:38 UTC 2021 x86_64 x86_64 x86_64 GNU/Linux
CPU - Intel(R) Core(TM) i9-10920X CPU @ 3.50GHz　×23
メモリ - 131.6 GB

分析したデータ

以下のデータはPythonのsinglecell解析ライブラリであるscanpyのtutorialであるTrajectory inference for hematopoiesis in mouseで用いられているものです。
Paul et al. (2015)
このデータは、6～8週齢のメスC57BL/6マウスの造血細胞のものです。
C57BL/6マウスについては以下を参考にしてください。

解析

解析は基本的にPalantirのチュートリアルに従います。
ソースコードは以下にあります。

QC（クオリティーチェック）

QCは以下の2つを行いました。

• ミトコンドリア、リボソームで足切り、フィルタリング
• Doublet除去

これらについて以下の記事で解説しています。

具体的には以下のように行いました。
発現した遺伝子が 100 未満の細胞、発現した細胞が3未満の遺伝子をフィルタリングします。

sc.pp.filter_cells(adata, min_genes=100)
sc.pp.filter_genes(adata, min_cells=3)

Doubletを除去します。

sc.pp.scrublet(adata, batch_key="sample")

その後正規化、HVG (Highly Variable Genes)の選択、PCA、Diffusion maps、umapを行います。
PCA、umap、Diffusion mapsはともに次元削減に使われるものですが、詳しくは以下の記事が参考になります。

データの中に入っていたannotaion情報を確認します。
annotationとは単一細胞解析において、遺伝子発現データに生物学的意味を付与するプロセスです。具体的には、数千～数百万の細胞データをクラスタリングした結果に対し、既知のマーカー遺伝子（例：CD34＝造血幹細胞、INS＝膵β細胞）やデータベース参照を用いて、細胞タイプ（神経細胞・免疫細胞など）や状態（未分化・分化中・アポトーシス期）を分類・命名します。

下の図はumap上でannotation情報を表示したものです。１９個のクラスターができました。

クラスター名	細胞タイプ	説明
1Ery - 6Ery	赤血球系統（Erythroid lineage）	赤芽球（proerythroblastやerythroblast）の分化段階を表し、初期から成熟段階まで含まれる。
7MEP	巨核球-赤血球前駆細胞（MEP）	赤血球と巨核球系統への分化が可能な前駆細胞。
8Mk	巨核球（Megakaryocyte）	血小板を産生する細胞。
9GMP, 10GMP	顆粒球-単球前駆細胞（GMP）	顆粒球や単球への分化が可能な前駆細胞。
11DC	樹状細胞（Dendritic Cell, DC）	抗原提示機能を持つ免疫細胞。
12Baso, 13Baso	好塩基球（Basophil）	アレルギー反応や炎症に関与する顆粒球。
14Mo, 15Mo	単球（Monocyte）	マクロファージや樹状細胞への分化が可能な免疫細胞。
16Neu, 17Neu	好中球（Neutrophil）	細菌感染に対する主要な防御役割を果たす顆粒球。
18Eos	好酸球（Eosinophil）	寄生虫感染やアレルギー反応に関与する免疫細胞。
19Lymph	リンパ球（Lymphocyte）	B細胞、T細胞、NK細胞などのリンパ系免疫細胞を含む可能性。

その後MAGICを利用してデータを補完し、その後Palantirを実行します。
MagicとはscRNA-seqデータのノイズ除去と遺伝子発現パターンの復元を行うアルゴリズムです。
Palantir実行のために一つだけ決定しなければならないものがあります。それが初期細胞です。ここでは造血幹細胞HSCです。HSCのマーカーとしてはCD34がよく知られており、今回もCD34が最もよく発現している細胞を初期細胞とします。

start_cell = adata.obs_names[adata[:, 'Cd34'].X.argmax()]

少し立ち止まって確認してみましょう。
初期細胞のumapでの位置は以下となります。

代表的なマーカー遺伝子の発現は以下の様になっています。これはMagicにより補正された値を用いてumapでプロットしたものです。

マーカー遺伝子	細胞
CD34	造血幹細胞（HSC）
MPO	好中球
GATA1	赤血球系
IRF8	単球系前駆細胞

CD34の発現が最も大きいものを初期細胞に選んだのに、CD34発現の図の赤色の部分と位置が少しずれています。したがって初期細胞の発現量は外れ値の可能性があります。なのでCD34の発現量を縦軸に、細胞IDを横軸にとったプロットを見てみましょう。横軸はデータの並び順に依存します。

CD34の発現量の最大値はほかに比べ際立って大きいわけではなく外れ値というわけではなさそうです。なのでこのまま進めます。

Palantirのcore関数を実行

このライブラリのかなめとなる関数を実行します。ここで疑似時間や細胞運命の計算をします。

pr_res = palantir.core.run_palantir(
    adata, start_cell, num_waypoints=500
)

結果は以下となります。

Pseudotimeは疑似時間であり、どれだけ初期細胞から時間がたったかを表します。これはつまり疑似時間が違えば違う細胞であることから、細胞のクラスタリングと似た結果になります。

Entropyは細胞の可塑性（他の運命への柔軟性）であり、Pseodptimeとおおむね逆の結果になることがわかります。つまり、疑似時間がたてばたつほどおおむね他の細胞への柔軟性は減るということです。プロットの結果は矛盾していません。

左下の２つはターミナル状態（終端状態）となる確率を表します。それぞれの細胞がW31350(Neu、好中球)になる確率、W38920(Ery、赤血球)になる確率を表します。

また、ここでのPseudotimeが0とは初期細胞を表します。しかしこれは上で示した初期細胞の位置と一致していません。したがってpseudotimeにおける初期細胞はpseudotime独自に算出されることがわかります。
どのように算出されているか見てみましょう

    # Determine the boundary cell closest to user defined early cell
    dm_boundaries = pd.Index(set(data_df.idxmax()).union(data_df.idxmin()))
    dists = pairwise_distances(
        data_df.loc[dm_boundaries, :], data_df.loc[early_cell, :].values.reshape(1, -1)
    )
    start_cell = pd.Series(np.ravel(dists), index=dm_boundaries).idxmin()
    if use_early_cell_as_start:
        start_cell = early_cell

まず各拡散成分の最大・最小値を持つ細胞（dm_boundaries）を抽出し、ユーザーで指定されたearly_cellとのユークリッド距離を計算した上で、最も近い境界細胞をstart_cellとして決定しています。もしオプションuse_early_cell_as_startがTrueの場合は、単にearly_cellがstart_cellとして使われます。
拡散成分とは前で計算した拡散マップの各軸の値をさします。各拡散成分の最大・最小値を持つ細胞（dm_boundaries）とは極端な細胞であり、つまり初期細胞やterminal細胞などをあらわす可能性が高い細胞です。この細胞の中から指定した細胞に最も近い細胞を初期細胞として選ぶということになっています。

次は疑似時間と終端状態となる確率のプロットを見てみます。

終端状態はW31350(Neu、好中球)、W38920(Ery、赤血球)の２つであり、実際に２つに分岐していることが確認できます。W31350(Neu、好中球)、W38920(Ery、赤血球)ともに0.1のあたりで大きな遺伝子発現の変化が起きていることがわかります。

それぞれの細胞の位置を確認してみましょう。左の点がW31350(Neu、好中球)、真ん中が設定した初期細胞、右の点がW38920(Ery、赤血球)です。

エントロピーを可視化したumap上に軌跡を視覚化します。

矢印は10個あるので１つが0.1に相当します。ここでのエントロピーとは他の運命への柔軟性を表しました。エントロピーの等高線と垂直に軌跡が描かれていることがわかります。また当然ですがエントロピーの低い方に向かっていくのがわかります。細胞が生まれて時間がたつほどほかの細胞になる確率はさがるという当然の結果です。

gene_trendとは

唐突ですがここでgene_trendというものを説明します。これは今からの議論で必要となるものです。
教科書的な説明は以下のようになります。
gene_trendsとは、細胞の擬似時間軸（pseudotime）に沿った各遺伝子の発現変化パターンを表現するデータ構造で、細胞が特定の分岐（branch）をたどる際に、各遺伝子がどのように発現しているかの時間的変化をモデル化し、定量的なトレンドとして抽出します。
しかしこれでは何を言っているかわからないと思います。
遺伝子発現量と比較しながら理解していきましょう。

gene_trendsと実際の遺伝子発現量には、以下の違いがあります

遺伝子発現量とはRNA-seqやscRNA-seqで直接計測される生データであり、個々の細胞/サンプルにおける遺伝子のmRNA量を反映します。
一方でgene_trends とは擬似時間軸に沿った発現パターンを数学的にモデル化した推定値です。mellon.FunctionEstimatorなどの平滑化アルゴリズムで算出され、細胞分化の動的プロセスを抽象化しています。mellonのgithubは以下です。

違いまとめ

	実測発現量	gene_trends
データ型	離散値（カウントデータ）	連続関数（平滑化曲線）
時系列解像度	個々の細胞の瞬間値	擬似時間軸上の連続的変化
遺伝子選択	発現量の絶対値で判定	トレンドの形状（傾き）で判定
生物学的解釈	瞬間的な細胞状態	細胞分化の動的プロセス

イメージとしてはgene_trendとは遺伝子発現量を時系列データにしたものと考えればよいでしょう。

それではgene_trendがどのように用いられているか見ていきましょう。
下は各細胞において疑似時間を横軸に、gene_trendを縦軸にプロットしたものです。

W31350(Neu)でCD34が時間とともに減少していません。これもCD34が初期細胞のマーカーとして適していないことを表しています。つまり本当の初期細胞のマーカーであれば時間がたつのに従ってgene_trendは下がってほしいわけです。エントロピーと軌跡のプロットを見ると、5,6番目の矢印でCD34がよく発現しているエリアを通ることがわかります。そしてこれはこのgene_trendと疑似時間のプロットとよく一致します。
そのほかのマーカーでは予想通りの結果となっています。例えばMpoを見てみましょう。W31350(Neu、好中球)は時間がたつのに従ってgene_trendが増えていきます。一方でW38920(Ery、赤血球)を見てみると時間がたつとgene_trendが減っていきます。つまりMpoはW31350(Neu、好中球)でよく発現しW38920(Ery、赤血球)で発現しないマーカーということになります。

マーカー遺伝子	細胞
CD34	造血幹細胞（HSC）
MPO	好中球
GATA1	赤血球系
IRF8	単球系前駆細胞

そしてこれは上の表に合致します。
一方IRF8ではW31350(Neu、好中球)で途中まで大きく発現し最後にはgene_trendが減ります。これもIRF8が単球系前駆細胞のマーカー遺伝子であることを考えると納得できます。

CD34とW31350(Neu、好中球)を詳しく見てみます。
発現量とgene_trendの関係が興味深いです。

次にNeu（好中球）に関連するgene_trendを可視化します。

遺伝子を1000個名前順に選択し、Pseudotimeと発現量のカーブを描きそれによってクラスタリングしてプロットします。クラスタリングはLeidenによるものです。様々なタイプの遺伝子発現動態がわかります。多くのクラスターでは発現量が増えていっています。これは分化が進み様々な遺伝子が発現していっていることを示しています。一方でcluster6では発現が少なくなっていっているのがわかります。従ってこの中に初期細胞のマーカー遺伝子が含まれているであろうことがわかります。
なおleidenアルゴリズムについては以下が参考になります。

まとめ

初期細胞選択には十分気を遣う必要があると感じました。初期細胞のマーカーを何個か試して適切なマーカーを選択する必要がありそうです。しかしそれさえできれば、高解像度で細胞運命を解析できて非常に使いやすいと感じました。今後の課題としては複数マーカーを組み合わせた初期細胞選択手法の検討や、分岐点近傍の細胞状態を詳細に解析するための補助ツールの必要性などがあると考えます。

(この記事は研究室インターンで取り組みました：https://kojima-r.github.io/kojima/)