前回の記事によりTrajectory解析に必要なMonocle3の環境ができました。いよいよこちらの環境を使ってTrajectory解析に進んでいこうと思います。この記事ではTrajectory解析の基本的な解析を学びます。ぜひやってみましょう。
macOS Monterey (12.4), Quad-Core Intel Core i7, Memory 32GB, R 4.3.0, Monocle3 1.3.7
公共データを用いたSingle Cell RNA-seq解析に関する初心者向け技術書を販売中
プログラミング初心者でも始められるわかりやすい解説！
RとSeuratで始めるSingle Cell RNA-seq解析！
公共データを用いたシングルセル ダイナミクス解析に関する初心者向け技術書を販売中
シングルセルデータの高度な解析であるTrajectory解析、RNA Velocity解析、空間トランスクリプトーム解析が環境構築方法から詳しく解説されています！
Trajectory解析とは？
Trajectory解析（軌跡解析、軌道解析とも呼ばれる）は、シングルセルRNAシーケンス（scRNA-seq）データの解析において、細胞が時間経過とともにどのように状態を変化させていくのか、その過程（軌跡）を推定する手法です。従来のバルクRNA-seqでは細胞集団の平均的な遺伝子発現しか見えませんでしたが、scRNA-seqによって個々の細胞のトランスクリプトームを捉えることができるようになり、Trajectory解析によって細胞の動的な変化をより詳細に理解することが可能になりました。
Hussein A et al., 2021 より引用
https://www.nature.com/articles/s41467-021-26282-z
Trajectory解析は、以下のような疑問に答えるのに役立ちます。
- 幹細胞はどのような細胞に分化していくのか？その経路は？
- 細胞はどのような中間状態を経て最終的な状態に到達するのか？
- 細胞の状態変化はどのような順序で起こるのか？
Trajectory解析では、各細胞の遺伝子発現プロファイルを基に、細胞間の類似性に基づいて細胞を並べ、擬似的な時間軸（pseudotime）を構築します。この擬似時間軸に沿って細胞を配置することで、細胞がどのような経路をたどって状態変化していくのかを可視化することができます。
Monocle 3とは？
Monocle 3 は、単一細胞トランスクリプトームデータ (scRNA-seq) の解析に特化したオープンソースのRパッケージです。主に、単一細胞遺伝子発現データのクラスタリングや分類、細胞の軌跡推定（擬似時間解析）、および差次的発現解析をサポートします。これにより、細胞の多様性や動的な発達過程を深く理解するための重要な洞察を得ることができます。
公式サイトより引用
https://cole-trapnell-lab.github.io/monocle3
論文
Monocle3はウェブサイトやドキュメントで詳細な情報やチュートリアルが公開されています。こちらをフォローする形で進めていきます。
https://cole-trapnell-lab.github.io/monocle3/docs/trajectories
Monocle3の基本的な使い方
それでは実際にMonocle3を使っていきます。まずはMonocle3でUMAPとt-SNEを描画するところまでやってみます。
データロード
Cao & Packerらの線虫データをロードしていきます。
https://www.science.org/doi/10.1126/science.aam8940
expression_matrix <- readRDS(url("<https://depts.washington.edu>:/trapnell-lab/software/monocle3/celegans/data/cao_l2_expression.rds"))
cell_metadata <- readRDS(url("<https://depts.washington.edu>:/trapnell-lab/software/monocle3/celegans/data/cao_l2_colData.rds"))
gene_annotation <- readRDS(url("<https://depts.washington.edu>:/trapnell-lab/software/monocle3/celegans/data/cao_l2_rowData.rds"))expression_matrix <- readRDS(...): 遺伝子発現行列を読み込みます。行が遺伝子、列が細胞です。cell_metadata <- readRDS(...): 細胞のメタデータ（細胞ID、実験条件など）を読み込みます。gene_annotation <- readRDS(...): 遺伝子の注釈情報（遺伝子ID、遺伝子名など）を読み込みます。
# Make the CDS object
cds <- new_cell_data_set(expression_matrix,
cell_metadata = cell_metadata,
gene_metadata = gene_annotation)
new_cell_data_set()関数がcell_data_setオブジェクトを作成する役割を担っています。この関数に、必要な情報を引数として渡しています。
前処理
Monocle 3におけるデータの前処理について説明します。PCAを行って次元削減につかう主成分の数を見極めていきます。このあたりはSeuratと同じですね。
# Pre-process the data
cds <- preprocess_cds(cds, num_dim = 100)
plot_pc_variance_explained(cds)
preprocess_cds(cds, num_dim = 100): この関数が前処理を行います。cds: 前に作成したcell_data_setオブジェクトを指定します。num_dim = 100: PCAで計算する主成分の数を100に指定しています。つまり、高次元の遺伝子発現データが100個の主成分に要約されます。
次元削減と可視化
それでは次元削減していきます。Monocle 3では、次元削減にデフォルトでUMAPが使用されています。reduction_method引数で"PCA"、"tSNE"などを指定することもできます。
cds <- reduce_dimension(cds)
データをプロットしてみます。
plot_cells(cds)
今回のデータセットでは、すでにラベルをつけてくれているので、そちらで色を塗り分けてみます。
plot_cells(cds, color_cells_by="cao_cell_type")
たくさんの細胞種がUMAPプロット上で確認できるようになりました。
特定の遺伝子または遺伝子群の発現量に基づいて細胞を色分けするには（SeuratでいうFeature Plot）、Monocle3ではplot_cellsにgenes引数をつけ、遺伝子名を指定します。
plot_cells(cds, genes=c("cpna-2", "egl-21", "ram-2", "inos-1"))
無事UMAPが可視化できるようになりました。t-SNEのコードも紹介しておきます。reduce_dimension(cds, reduction_method="tSNE") でできますので試してみてください。
cds <- reduce_dimension(cds, reduction_method="tSNE")
plot_cells(cds, reduction_method="tSNE", color_cells_by="cao_cell_type")
Trajectory解析の実装方法
ここまででMonocle3の使い方を説明しました。いよいよTrajectory解析に入っていきます。
データの読み込み
データを読み込みます。こちらで使っているのは線虫（C. elegans）の胚発生における神経細胞のデータセットになります。
expression_matrix <- readRDS(url("<https://depts.washington.edu>:/trapnell-lab/software/monocle3/celegans/data/packer_embryo_expression.rds"))
cell_metadata <- readRDS(url("<https://depts.washington.edu>:/trapnell-lab/software/monocle3/celegans/data/packer_embryo_colData.rds"))
gene_annotation <- readRDS(url("<https://depts.washington.edu>:/trapnell-lab/software/monocle3/celegans/data/packer_embryo_rowData.rds"))
cds <- new_cell_data_set(expression_matrix,
cell_metadata = cell_metadata,
gene_metadata = gene_annotation)
前処理
先ほどと同じように、preprocess_cdsで前処理を行っていきます。今回はバッチ補正も行いますalign_cds()関数を用いてバッチ補正を行っています。
cds <- preprocess_cds(cds, num_dim = 50)
cds <- align_cds(cds, alignment_group = "batch", residual_model_formula_str = "~ bg.300.loading + bg.400.loading + bg.500.1.loading + bg.500.2.loading + bg.r17.loading + bg.b01.loading + bg.b02.loading")
cds <- align_cds(cds, alignment_group = "batch", residual_model_formula_str = "~ bg.300 + bg.400 + loading")
cds: 前処理後のcell_data_setオブジェクトを指定します。alignment_group = "batch": バッチ情報を格納しているcolData(cds)の列名を指定します。この例では、”batch”という列にバッチ情報が格納されていることを示しています。align_cds()関数は、この情報に基づいてバッチ間の違いを補正します。residual_model_formula_str = "~ bg.300 + bg.400 + loading": この引数は、連続的な影響を補正するための数式を指定します。この例では、”bg.300″、”bg.400″、”loading”という3つの変数の影響を補正しています。
次元削減と可視化
それでは次元削減していきます。Monocle3ではUMAPによる次元削減が強く推奨されていますので、デフォルトの引数のままreduce_dimension(cds)を動かしていきます。Trajectory解析では、細胞の分化経路という大域的な構造を捉えることが重要であるため、UMAPの方が適していると考えられます。
cds <- reduce_dimension(cds)
plot_cells(cds, label_groups_by_cluster=FALSE, color_cells_by = "cell.type")
推定した細胞の分化軌跡において、繊毛ニューロンに関連する遺伝子群 (che-1, hlh-17, nhr-6, dmd-6, ceh-36, ham-1) の発現がどのように変化するかを plot_cells() 関数を使って可視化してみます。
ciliated_genes <- c("che-1",
"hlh-17",
"nhr-6",
"dmd-6",
"ceh-36",
"ham-1")
plot_cells(cds,
genes=ciliated_genes,
label_cell_groups=FALSE,
show_trajectory_graph=FALSE)
細胞をクラスター化
細胞を異なる軌跡に分割するためにクラスタリングを使用します。これは、異なる種類の細胞が異なる発生軌跡を持つ可能性を考慮するためです。
cds <- cluster_cells(cds)
plot_cells(cds, color_cells_by = "partition")
以下のようなグラフが出ればクラスタリングがきちんとできています。
Trajectoryグラフを理解する
次に、learn_graph()関数を使用して、各パーティション内に主グラフを当てはめます。learn_graph(cds) では前のステップでクラスタリングされた細胞データ（cds）をもとに、細胞の遷移関係をグラフとして学習します。
cds <- learn_graph(cds)
plot_cells(cds,
color_cells_by = "cell.type",
label_groups_by_cluster=FALSE,
label_leaves=FALSE,
label_branch_points=FALSE)
plot_cells(...):color_cells_by = "cell.type": 細胞の種類ごとに色分けして表示します。label_groups_by_cluster=FALSE: クラスタごとにラベルを付けません。label_leaves=FALSE: グラフの末端（葉）にラベルを付けません。label_branch_points=FALSE: グラフの分岐点にラベルを付けません。
擬似時間に並べ替える
それでは、発生プログラムの進行状況に応じて細胞を順序付けてみましょう。
💡擬似時間について
細胞はそれぞれ異なるタイミングで変化するため、従来の時間軸だけではその変化を正確に捉えきれません。疑似時間という概念を用いることで、この非同期性を考慮し、各細胞の進行度合いを正確に把握することができます。
疑似時間に基づいて細胞を順序付けることで、細胞がどのように変化していくのか、その軌跡を再構築し、分化経路や発生過程を明らかにすることができます。さらに、疑似時間と遺伝子発現の変化を関連付けることで、発生や分化に関わる重要な遺伝子を特定することができ、疾患のメカニズム解明や治療法開発に役立ちます。
plot_cells(cds,
color_cells_by = "embryo.time.bin",
label_cell_groups=FALSE,
label_leaves=TRUE,
label_branch_points=TRUE,
graph_label_size=1.5)
このグラフは、細胞がどのように状態変化していくのかを示すもので、ノード（点）とエッジ（線）で構成されています。
- ノード: 細胞の状態を表します。
- エッジ: 細胞の状態遷移を表します。
グラフには、特別なノードがいくつか存在します。
- リーフノード（薄い灰色の円）: 軌跡の末端、つまり細胞の最終的な状態を表します。
- 分岐ノード（分岐ノード）: 細胞が複数の状態へと分岐する可能性のある地点を表します。
グラフは完全には接続されていないことに注意してください。異なるパーティション内の細胞は、グラフの異なるコンポーネントにあります。番号が付けられた円は、グラフ内の特別な点を示しています。
初期の細胞がどこに位置するのかがわかったので、order_cells() を呼び出します。これは、各細胞が疑似時間でどこに位置するのかを計算します。そのためには、order_cells() は軌跡グラフのルートノードを指定する必要があります。1つ以上のルートノードを選択するためのグラフィカルユーザーインターフェースが起動します。
ちょっと分かりづらいですが、ルートノード1は赤丸あたりです。ここを選んでいきます。
cds <- order_cells(cds)ルートノードを選んだ後に、「Done」をクリックします。
その後”pseudotime”でプロットすると、ルートノード起点にトラジェクトリーがプロットされます。
plot_cells(cds,
color_cells_by = "pseudotime",
label_cell_groups=FALSE,
label_leaves=FALSE,
label_branch_points=FALSE,
graph_label_size=1.5)
最後に
いかがだったでしょうか。Trajectory解析はそこまで解析方法は難しくないですが、解釈にいくつか難しさが存在するかもしれません。じっくりと各手法の意味を理解することが重要です。今後もしっかり学んでいきましょう。
公共データを用いたSingle Cell RNA-seq解析に関する初心者向け技術書を販売中
プログラミング初心者でも始められるわかりやすい解説！
RとSeuratで始めるSingle Cell RNA-seq解析！
公共データを用いたシングルセル ダイナミクス解析に関する初心者向け技術書を販売中
シングルセルデータの高度な解析であるTrajectory解析、RNA Velocity解析、空間トランスクリプトーム解析が環境構築方法から詳しく解説されています！
