RNA Velocity解析をやってみたくないでしょうか。RNA Velocity解析はシングルセルデータから細胞の状態遷移を把握することが可能です。この記事ではRNA Velocityの初歩的な操作を学びます。この記事によってRNA Velocityに入門してみましょう！
Google Colab, Python 3.10.12, scvelo 0.3.3
公共データを用いたSingle Cell RNA-seq解析に関する初心者向け技術書を販売中
プログラミング初心者でも始められるわかりやすい解説！
RとSeuratで始めるSingle Cell RNA-seq解析！
公共データを用いたシングルセル ダイナミクス解析に関する初心者向け技術書を販売中
シングルセルデータの高度な解析であるTrajectory解析、RNA Velocity解析、空間トランスクリプトーム解析が環境構築方法から詳しく解説されています！
RNA Velocity解析とは？
RNA Velocity解析は、一言で言うと、個々の細胞がどのような状態に変化していくのか（分化や状態変化の方向と速度）を予測する技術です。従来のRNA-Seqは細胞のスナップショット（静的な情報）しか捉えられませんでしたが、RNA Velocityは細胞内のmRNAのスプライシング状態に着目することで、時間的な変化（動的な情報）を捉えることを可能にしました。
https://www.nature.com/articles/s41586-018-0414-6
Bergen V, et al. Generalizing RNA velocity to transient cell states through hidden time. Nat Biotechnol. 2020;38(12):1408-1414. より引用
RNA Velocity解析では、スプライシングの過程に着目し、未成熟なRNA（unspliced RNA）と成熟したRNA（spliced RNA）の比率を分析します。unspliced RNAとspliced RNAの比率を比較することで、各遺伝子の発現量の変化の方向（増加傾向か減少傾向か）と速度を推定し、細胞の状態変化の方向と速度をベクトルとして表現することができます。これにより、細胞の分化過程や状態変化のダイナミクスを可視化し、理解を深めることができるのです。
RNA Velocity解析によって以下のようなことがわかります。
- 細胞の分化経路の予測: 例えば、幹細胞がどのような細胞に分化していくのか、その経路を予測することができます。
- 細胞の状態変化の方向と速度の把握: 細胞がどのような状態に変化していくのか、その変化のスピードを把握することができます。
- トランジション状態にある細胞の特定: 分化の途中段階など、過渡的な状態にある細胞を特定することができます。
scVeloとは？
scVeloは、シングルセルRNAシーケンスデータを用いたRNA velocity解析のためのPythonパッケージです。従来のRNA velocity解析手法を拡張し、より高度な解析を可能にしています。
https://scvelo.readthedocs.io/en/stable/perspectives/index.html
より引用
従来のRNA velocity解析では、定常状態（steady-state）を仮定していましたが、scVeloでは動的モデルを用いることで、より一般的な細胞状態の変化を捉えることができます。これにより、トランジション状態にある細胞や、複雑な分化経路を持つ細胞集団の解析が可能になっています。
scVeloはPythonパッケージとして提供されており、scVeloのウェブサイトやドキュメントで詳細な情報やチュートリアルが公開されています。こちらをフォローする形で進めていきます。
https://scvelo.readthedocs.io/en/stable/VelocityBasics.html
また上記チュートリアルでのGoogle Colabでエラーが出るところを直したColabリンクはこちらになります。ご自身の環境でもColab上で試されてもどちらでも構いません。
https://colab.research.google.com/drive/1fiTfbrYkoIk6bSf3eWDSBynyVKTf2n3U?usp=sharing
RNA Velocity解析の実装方法
それでは実際にRNA Velocity解析を実装していきましょう。
ライブラリのインストール
pip install scvelo --upgrade --quietimport scvelo as scvscVeloの設定
scv.settings.verbosity = 3 # show errors(0), warnings(1), info(2), hints(3)
scv.settings.presenter_view = True # set max width size for presenter view
scv.set_figure_params('scvelo') # for beautified visualizationscv.settings.verbosity = 3- scVeloのログレベルを設定します。
0: エラーのみ表示1: 警告とエラーを表示2: 情報、警告、エラーを表示3: ヒント、情報、警告、エラーを表示
- この設定では、ヒントを含む全てのログが表示されます。
- scVeloのログレベルを設定します。
scv.settings.presenter_view = True- プレゼンタービューを有効にします。
- プレゼンタービューでは、プロットの最大幅が制限され、プレゼンテーションなどで見やすくなります。
scv.set_figure_params('scvelo')- scVeloのデフォルトの図パラメータを設定します。
- これにより、scVeloのビジュアライゼーションがより美しくなります。
Load the Data
それでは、データをロードしていきます。
scv.datasets.pancreas() で読み込まれるデータは、Baron et al. (2016) によって公開されたマウスの膵臓のシングルセルRNAシーケンスデータです。このデータには、様々な種類の細胞（α細胞、β細胞、δ細胞など）が含まれており、膵臓の発生や機能を研究するために広く使われています。
adata = scv.datasets.pancreas()
adata読み込まれたデータは adata という変数に格納されます。 adata はAnnDataオブジェクトで、シングルセルデータとそのメタデータを格納するための標準的なデータ構造です。このデータには、遺伝子発現量 (X)、細胞のクラスタ情報 (obs)、遺伝子の情報 (var)、スプライシングされたmRNAとスプライシングされていないmRNAの情報 (layers) などが含まれています。
adataが気になる方はこちらのanndataのドキュメントをご参照ください。
https://anndata.readthedocs.io/en/stable
scv.pl.proportions(adata)scv.pl.proportions(adata): AnnDataオブジェクト(adata)に格納されているデータのスプライシングされたmRNAとスプライシングされていないmRNAの割合を計算し、円グラフと棒グラフで可視化する関数です。
出力は2つのグラフで構成されています。
- 円グラフ: データ全体におけるスプライシングされたmRNAとスプライシングされていないmRNAの割合を示しています。
- 棒グラフ: 各細胞クラスターにおけるスプライシングされたmRNAとスプライシングされていないmRNAの割合を示しています。
一般的に、スプライシングされていないmRNAの割合は、使用されるプロトコル (Drop-Seq、Smart-Seqなど) によって異なり、10〜25% 程度です。
Preprocess the Data
scVeloでは、RNA velocityを正確に推定するために、前処理として遺伝子発現データのフィルタリング、正規化、モーメントの計算などを行います。
scv.pp.filter_and_normalize(adata, min_shared_counts=20, n_top_genes=2000)
scv.pp.moments(adata, n_pcs=30, n_neighbors=30)
scv.pp.filter_and_normalize(adata, min_shared_counts=20, n_top_genes=2000) について
min_shared_counts=20: 整数値 (デフォルト値: None)。遺伝子フィルタリングの閾値を指定します。この値以上のカウント数（スプライシングされたmRNAとスプライシングされていないmRNAの合計）を持つ遺伝子のみが保持されます。低発現の遺伝子を除去することで、ノイズを減らし、計算の効率と精度を向上させることができます。n_top_genes=2000: 整数値 (デフォルト値: None)。保持する変動の大きい遺伝子の数を指定します。上位n_top_genes個の変動の大きい遺伝子が選択されます。変動の大きい遺伝子は、細胞の状態や変化を捉えるのに重要であると考えられます。多くの場合、2000〜5000個の遺伝子を選択します。
scv.pp.moments(adata, n_pcs=30, n_neighbors=30) について
- 主成分分析 (PCA): 高次元データである遺伝子発現データを、低次元の主成分空間に射影します。
n_pcs=30は、上位30個の主成分を使用することを意味します。 - 最近傍探索: 各細胞の近傍にある細胞を探索します。
n_neighbors=30は、各細胞について30個の近傍細胞を探索することを意味します。
これらの前処理により、ノイズが除去され、細胞間の関係が明確化されたデータが得られます。これにより、より正確なRNA velocity解析が可能になります。
ベストプラクティスが知りたい場合は、こちらの論文を参考にしてください。
https://www.embopress.org/doi/full/10.15252/msb.20188746
Estimate RNA velocity
scv.tl.velocity(adata)
scv.tl.velocity_graph(adata)
このコードは、RNA velocityを計算しています。scVeloでは、各遺伝子について、スプライシングされていないmRNAとスプライシングされたmRNAの量比が一定になる定常状態を仮定します。各細胞の遺伝子発現量を、定常状態における量比と比較することで、RNA velocityを推定します。
その後、scv.tl.velocity_graph(adata)で速度グラフを計算しています。速度グラフは、細胞間遷移の可能性を表現したもので、RNA velocity解析において重要な役割を果たします。
Project the velocities
それではRNA Velocityを可視化してみましょう。
scv.pl.velocity_embedding_stream(adata, basis='umap')
scVeloで計算されたRNA velocityを、UMAPでストリームラインとして可視化しています。以下の出力が得られます。
このRNA velocity解析の結果から、膵臓の細胞分化は、Ductal細胞から始まり、Ngn3を発現する前駆細胞を経て、様々な内分泌細胞へと分化していく複雑な過程であることが示唆されます。
次に、RNA velocityを、ベクトル場として可視化して見ましょう。
scv.pl.velocity_embedding(adata, arrow_length=3, arrow_size=2, dpi=120)
引数は以下の解説になります。
arrow_length=3: 矢印の長さを調整します。arrow_size=2: 矢印のサイズを調整します。dpi=120: 画像の解像度を調整します。
ストリームラインとして表示された速度ベクトル場は、発生過程を詳細に理解するのに役立ちます。これは、導管細胞と内分泌前駆細胞の循環集団を正確に描写します。さらに、細胞系譜の決定、細胞周期の終了、および内分泌細胞の分化における細胞状態を明らかにします。
最後に
いかがだったでしょうか。RNA Velocityを使えるようになるとシングルセルデータの解析がぐっと広がります。次はRNA Velocityのデータ準備の方法について解説します！ぜひともじっくり学んでみてください！
参考
https://scvelo.readthedocs.io/en/stable/VelocityBasics.html
公共データを用いたSingle Cell RNA-seq解析に関する初心者向け技術書を販売中
プログラミング初心者でも始められるわかりやすい解説！
RとSeuratで始めるSingle Cell RNA-seq解析！
公共データを用いたシングルセル ダイナミクス解析に関する初心者向け技術書を販売中
シングルセルデータの高度な解析であるTrajectory解析、RNA Velocity解析、空間トランスクリプトーム解析が環境構築方法から詳しく解説されています！
