- 著者: Gioele La Manno, Ruslan Soldatov, Amit Zeisel, Emelie Braun, Hannah Hochgerner, Viktor Petukhov, Katja Lidschreiber, Maria E. Kastriti, Peter Lönnerberg, Alessandro Furlan, Jean Fan, Lars E. Borm, Zehua Liu, David van Bruggen, Jimin Guo, Xiaoling He, Roger Barker, Erik Sundström, Gonçalo Castelo-Branco, Patrick Cramer, Igor Adameyko, Sten Linnarsson, Peter V. Kharchenko
- Corresponding author: Sten Linnarsson (Karolinska Institutet) / Peter V. Kharchenko (Harvard Medical School)
- 雑誌: Nature
- 発行年: 2018
- Epub日: 2018-08-23
- Article種別: Original Article
- PMID: 30089906
背景
単一細胞RNAシーケンス (scRNA-seq; single-cell RNA sequencing) は個々の細胞の転写産物発現量を高精度・高感度・高スループットで定量する技術として急速に普及し、細胞生物学に革命をもたらした。10x Genomics Chromiumプラットフォームを用いたマッシブ並列デジタル転写プロファイリング (Zheng et al. NatCommun 2017) など、近年の技術的進歩により数千から数万細胞規模の網羅的解析が日常的に実施可能となっている。しかし、これらの技術が提供するデータは特定の測定時点における細胞の静的なスナップショットに過ぎず、細胞分化・組織発生・再生などの時間依存的プロセスにおける細胞状態の変化の方向や速度を直接的に解析する能力を本質的に欠いていた。
細胞の分化プロセスは数時間から数日のタイムスケールで進行し、これはmRNAの典型的な半減期と同等のオーダーである。バルクRNA-seqを用いた先行研究においては、一般的なoligo-dTプライマープロトコルでもリードの約15%がイントロン配列を含む未スプライスの前駆体mRNA (unspliced mRNA) に由来することが報告されており (Gaidatzis et al. 2015 Nat Biotechnol、14.6% in bulk)、この現象がscRNA-seqにも拡張適用可能であることが示唆されていた。また、連動する前駆体mRNAと成熟mRNAの動態が転写後調節の手がかりとなり得ることは、バルク解析においても理論的に提示されていた (Zeisel et al. 2011 Mol Syst Biol)。さらに、単一RNA標本からin vivoのmRNA動態を推定するSnapShot-Seq法も報告されていた (Gray et al. 2014 PLOS ONE (Public Library of Science ONE))。
しかしながら、これらの知見はバルク解析に限定されており、単一細胞レベルの不均一な細胞集団においてunsplicedとsplicedのmRNA比から転写動態を推定し、個々の細胞が将来的にどの発現状態へと遷移するかを予測する計算手法は確立されていなかった。がんゲノミクスや発生生物学において重要な細胞状態遷移の動的追跡——例えば腫瘍微小環境における細胞状態不均一性 (Schaff et al. CancerRes 2021) やエピゲノム不均一性を起源とする腫瘍進化 (Laisne et al. NatRevCancer 2025)——を追跡するためには、静的スナップショット解析を超えた動的情報の gap in knowledge が存在していた。代謝標識などの追加実験を必要とせず、既存の標準的なscRNA-seqデータのみから細胞の将来の発現軌跡を直接予測できる計算フレームワークが強く求められていた。
目的
本研究の目的は、標準的なscRNA-seqデータに内在するunspliced/spliced mRNA比の情報から遺伝子発現状態の時間微分である「RNA velocity (RNA速度)」を直接推定する数理モデルおよびソフトウェアフレームワーク「velocyto」を開発することである。SMART-seq2 (Switching Mechanism At the 5’ end of RNA Template sequencing 2)、inDrop、10x Genomics Chromiumを含む複数のscRNA-seqプロトコルに適用し、個々の細胞の将来の発現状態を数時間スケールで予測できることを実証する。さらに、マウス海馬発生の複雑な分岐系統樹の再構築と、ヒト胎児脳における転写キネティクスの多様性の解明を通じて、追加の実験的介入なしに細胞の分化方向と運命決定ダイナミクスを定量的に記述する手法の汎用性と正確性を確立することを目指す。
結果
技術的検証：unsplicedリードの普遍的検出と一次反応速度論との整合性: SMART-seq2、STRT/C1、inDrop、10x Genomics Chromiumの4プロトコルを網羅的に解析した結果、すべてのプラットフォームで全リードの15〜25%がunsplicedイントロン配列を含んでいた (Fig. 1a)。これはバルクRNA-seqにおける先行報告 (14.6%) および単一細胞RNA-seqにおける観察値 (約20%) と一致する。HEK293細胞を用いた4sU代謝標識実験 (n=998 genes、n_technical=2、n_biological=2) において、解析遺伝子の83.4%が一次反応速度論と一致する発現時間経過を示し (Fig. 1)、unsplicedリードが新生転写物を代表していることが独立して実証された。スプライシング速度 β と分解速度 γ の最尤推定から導出した特性時間定数 τ (63.2%の漸近値に到達するのに要する時間) の分布は、spliced/unspliced比の変化が10〜100分のタイムスケールで検出可能であることを示した。マウス肝臓の概日リズム時系列データでは、unspliced mRNA量が次の時点のspliced mRNAレベルを一貫して予測し (Fig. 1e)、Fgf1 (fibroblast growth factor 1) およびCbs (cystathionine beta-synthase) のフェーズポートレートにおいて転写誘導時の上方偏差・抑制時の下方偏差が確認された。一次微分方程式モデルを解いて外挿することで、24時間の概日サイクルの進行方向が正確に再現された (Fig. 1h)。
クロマフィン細胞分化における分化方向予測の実証: マウス胎生期E12.5のクロマフィン細胞データ (n=385 cells) に本手法を適用した結果、RNA velocityはSchwann cell precursor (SCP; シュワン細胞前駆体) からクロマフィン細胞への一方向性の分化経路を正確に予測した (Fig. 2d)。これはPLP1 (proteolipid protein 1)-CreERT2を用いた独立した系譜追跡実験結果と完全に一致した。転写抑制される遺伝子Serpine2 (serpin family E member 2) は分化に伴い約2-fold低下し、転写誘導される遺伝子Chga (chromogranin A) は約4-fold上昇した。それぞれのフェーズポートレートは定常平衡線からの期待される偏差パターンを明確に示した (Fig. 2b, c)。EdU (5-ethynyl-2’-deoxyuridine) パルス標識実験に基づく外挿時間の定量解析では、最適外挿時間 t* の分布モードが2.1時間であり、単一細胞レベルで2.5〜3.8時間先の将来状態を高い信頼性で予測できることが実証された (Fig. 2f, g)。また細胞周期の動態についても、PCAおよび細胞周期関連遺伝子の集中解析においてRNA velocityが細胞周期の位相を正確に捉えることが確認された。
マウス海馬発生における複雑な分岐系統樹の自動再構築: 発達中のマウス海馬から得られた大規模データ (n=18,213 cells、P0: 8,113 cells、P5: 10,100 cells) において、t-SNE空間上の複雑な分化マニフォールドにRNA velocityを適用した (Fig. 3)。既知のマーカー遺伝子を用いて、アストロサイト (Aqp4 (aquaporin 4) 陽性)、OPC (oligodendrocyte precursor cell; オリゴデンドロサイト前駆細胞、Pdgfra (platelet-derived growth factor receptor alpha) 陽性)、歯状回顆粒ニューロン (Prox1陽性)、海馬各領域 (CA1/CA2/CA3/subiculum/hilus) のニューロンを同定した。RNA velocityは放射状グリア (radial glia; Hes1 (hairy and enhancer of split 1)・Hopx (homeodomain-only protein homeobox) 陽性) を起点として各主要系統へと分岐する強い方向性フローを明確に示した (Fig. 3c)。マルコフランダムウォークモデルの正方向・逆方向シミュレーションにより、事前の生物学的知識なしに分化終点と起点 (radial glia) が自動的に同定された。特に、pre-OPCからOPCへと至る狭い分化通路においてコミットメントのダイナミクスが単一細胞遷移確率として可視化され、通路入口の細胞は複数の運命に向かう可能性を持つのに対し、通路内の細胞はほぼ完全にOPC運命に固定されていることが示された (Fig. 3d)。また歯状回顆粒ニューロン/CA pyramidal neuronの分岐点において、Prox1 (prospero homeobox 1) の発現活性化が運命決定の方向を既に傾けていることが捉えられた (Fig. 3e)。
ヒト胎児脳における転写キネティクスの多様性と持続転写現象の同定: 受胎後10週のヒト前脳グルタミン酸作動性ニューロン系譜 (n=1,720 cells) において、増殖性前駆細胞 (radial glia、SOX2 (SRY-box transcription factor 2) 陽性) から中間神経芽細胞 (EOMES (eomesodermin) 陽性) を経てSLC17A7 (solute carrier family 17 member 7; 小胞型グルタミン酸トランスポーター) 陽性の成熟興奮性ニューロンへと至る分化経路を方向性をもって記述した (Fig. 4a)。RNA velocityにより予測された新規マーカー遺伝子 (CLU (clusterin)・FBXO32 (F-box only protein 32): ventricular zone; UNC5D (unc-5 netrin receptor D): intermediate zone; SEZ6 (seizure-related 6 homolog)・RBFOX1 (RNA-binding Fox-1 homolog 1): cortical plate) の発現パターンを、ヒト胎児皮質切片のRNAscope多重in situ hybridizationにより検証し、空間的な皮質層構造とscRNA-seqの擬似時間軸上の分布が密接に対応することを確認した (Fig. 4b, c)。主曲線 (principal curve) を用いた擬似時間解析では、unspliced mRNAが上昇期・下降期ともに一貫してspliced mRNAに先行するキネティクスが観察された (Fig. 4d)。遺伝子ごとの転写キネティクスは多様であり、RNASEH2B (ribonuclease H2 subunit B) のように spliced-unspliced 乖離が小さい高速動態の遺伝子がある一方で、DCX (doublecortin)・ELAVL4 (ELAV-like RNA-binding protein 4)・STMN2 (stathmin 2) は初期転写バーストの後に転写が持続する「オーバーシュート (overshoot)」動態を示した。分解速度の遅い遺伝子を迅速に誘導するための持続転写調節機構が、ヒト胚の発生過程において機能していることが本研究で直接示された。
多プラットフォームにわたる頑健性の検証: マウス骨髄好中球成熟 (inDrop、n=3,018 cells)、光刺激後の視覚野ニューロン応答 (inDrop、n=952 cells)、腸管上皮成熟 (10x、n=2,683 cells)、オリゴデンドロサイト分化 (10x、n=6,307 cells) など複数の組織・生物学的システムでRNA velocity推定を検証した結果、すべてのシステムで一貫して妥当な分化方向が捉えられた。分化過程で発現変動を示す遺伝子のうち75%においてRNA velocity推定値と実測の擬似時間微分の間に正の Spearman 相関 (中央値 r=0.65、n=4186 genes across datasets) が確認され (Extended Data Fig. 8)、数理モデルの予測精度が独立して実証された。モデルパラメータや細胞・遺伝子のサブサンプリングに対する感度分析ではvelocity推定は全般に頑健であり、可視化の近傍サイズ (k) が最も影響の大きいパラメータと特定された。マウス多組織データセット (Tabula Muris) の解析では11%の遺伝子で組織サブセット間に異なるγ値が観察され (Extended Data Fig. 3)、組織特異的なオルタナティブスプライシングや分解速度の多様性を反映していた。
考察/結論
先行研究との違い: 本研究は、従来のscRNA-seq解析が捉えていた静的な細胞状態のスナップショットという限界を根本的に克服し、unspliced/spliced mRNAの比から遺伝子発現の時間微分を直接算出する「RNA velocity」という動的指標を導入した点で、これまでの研究と大きく異なる。既報のバルクRNA-seqに基づく転写速度推定は集団平均に限られており、個々の細胞の分化方向の多様性や運命決定の分岐点ダイナミクスを捉えることができなかった点で本手法と対照的である。従来の擬似時間 (pseudotime) 解析手法は細胞状態の連続性を静的に整列するのみで、細胞が移動する「方向」と「速度」というベクトル情報を直接提供できず、系統樹の根方向と先端方向の区別には外部知識が必要であった。本研究はこの課題を克服し、代謝標識実験を必要とせず既存のscRNA-seqデータのみからキネティクスパラメータを推定できる計算フレームワークを提供した点で、相違は根本的である。
新規性: 本研究の核心的な新規性は、unspliced/spliced mRNA比という既存データ内に潜在していた動的情報を利用してRNA velocityという概念を定式化し、各細胞の将来発現状態を数時間先まで予測する計算基盤を確立した点にある。これまで報告されていなかったヒト胎児脳グルタミン酸作動性ニューロン発生における「オーバーシュート型」持続転写調節を同定し、分解速度の遅い遺伝子の迅速な誘導機構をヒト胚で初めて直接示したことは、発達神経科学における新規の知見である。また、マルコフランダムウォークモデルとの統合により、系統樹の根状態・終点状態・コミットメント段階が生物学的事前知識なしに完全自動で同定できることを示した点も、従来にはない新規なアプローチである。さらに、novel な計算ツール velocyto をR/Pythonの両言語でオープンソースとして公開し、多様なscRNA-seqプロトコルへの汎用性を実証したことも重要な貢献である。
臨床応用: RNA velocityは発生生物学の枠を超え、がん微小環境における細胞状態遷移や薬剤耐性獲得プロセスの動的追跡に向けた臨床応用が期待される。再生医療においてはiPS細胞からの目的細胞への分化誘導プロセスの品質管理への橋渡しとなる技術としての意義を持つ。免疫チェックポイント療法や細胞療法 (CAR-T細胞療法など) においても、免疫細胞の活性化・疲弊状態間の遷移方向をRNA velocityで推定することで治療応答の動的予測への臨床的意義が見込まれる。先天性発達障害や神経疾患における異常な転写キネティクスの同定にも応用可能であり、ゲノム医療における診断・治療標的の探索に貢献することが期待される。
残された課題: 今後の検討課題として、現モデルが仮定する単純一次反応速度論の限界が挙げられる。組織・細胞状態依存的なスプライシング速度や分解速度の変化、オルタナティブスプライシング (11%の遺伝子で二重γが観察) など、より複雑な転写調節がvelo city推定のバイアスとなる limitation がある。また、数時間程度の予測時間窓を超える長期的な細胞運命予測には逐次外挿による誤差蓄積が課題となる。シーケンス深度が低いデータではunsplicedリードの検出感度が低下し推定の信頼性が損なわれる問題も残されており、低深度データへの統計的補正手法の確立が今後の重要な研究方向性となる。著者らはfuture researchとして、RNA velocityに基づいたマニフォールド学習アルゴリズムの開発により、マニフォールドとその上のキネティクスを同時に推定できるフレームワークの構築を展望している。
方法
数理モデルの構築: スプライス済み成熟mRNA量 (s) の時間微分であるRNA velocityを、unspliced mRNA (u) からの産生 (スプライシング速度 β) と成熟mRNAの分解 (分解速度 γ) のバランスとして記述する一次反応速度論モデルを構築した (ds/dt = βu - γs)。定常状態では転写速度 α が一定のとき u = γs という固定傾きの関係が成立し、各遺伝子の平衡係数 γ はextreme quantile regression (最極端分位回帰) によりロバストに推定した。観察されたu/s比が定常状態の平衡線から外れる方向と大きさに基づき、個々の細胞における各遺伝子のRNA velocityを算出した。定常状態から大きく外れた遺伝子に対しては、遺伝子構造に基づく代替フィット手法 (structure-based fit) も開発した。
ソフトウェア開発とデータ処理: BAMファイルからリードをspliced・unspliced・ambiguousの3カテゴリに分類してカウントするコマンドラインツール velocyto.py (Python版) および velocyto.R (R版) を開発した。ヒトゲノムアノテーション (GRCh37.82) およびマウスゲノムアノテーション (GRCm38 (Genome Reference Consortium Mouse Build 38), Ensembl release 84) を使用し、UMI (unique molecular identifier) ベースのプロトコルでは細胞バーコードとUMI配列の組み合わせで重複除去を実施し、SMART-seq2プロトコルでは各リードを独立した分子として処理した。
検証データセット: 以下の実験系を用いた。(1) マウス肝臓の概日リズム時系列バルクRNA-seq (SRA (Sequence Read Archive): SRA025656、24時間タイムコース)、(2) マウスクロマフィン細胞分化データ (SMART-seq2、GEO (Gene Expression Omnibus): GSE99933、n=385 cells、胎生期E12.5)、(3) マウス海馬発生データ (10x Genomics Chromium、GEO: GSE104323、n=18,213 cells、P0+P5)、(4) ヒト胎児前脳データ (10x Chromium V2、SRA: SRP129388、n=1,720 cells、受胎後10週)。ヒト胎児前脳組織は倫理承認 (Regional Ethics Vetting Board Stockholm 2007/1477-31/3) に基づいて取得した。検証実験として、ヒトHEK293 (Human Embryonic Kidney 293 cell line, ATCC CRL-1573) に4-チオウリジン (4sU) を5・15・30分間添加した代謝標識実験 (GEO: GSE115813) を実施した。マウスはCD1 (CD-1 outbred albino mouse model, Charles River) 系統およびHtr3a-EGFPトランスジェニックマウスを使用した。
統計解析と可視化: RNA velocityベクターをPCA (principal component analysis; 主成分分析) またはt-SNE (t-distributed stochastic neighbor embedding) の低次元空間に射影して可視化した。Spearman順位相関係数によりγ推定の頑健性を評価し、Pearson相関係数を用いて細胞間距離を計算した。Louvain community detection algorithmによりクラスタリングを実施した。視覚野データセットでは負の二項一般化線形モデルと尤度比検定を用いて刺激関連遺伝子を選択し、有意水準はspliced分子でp=0.001、unspliced分子でp=0.03とした。分化終点および根状態の自動同定にはマルコフランダムウォーク (Markov random walk) モデルを適用した。