- 著者: Sara Aibar, Carmen Bravo González-Blas, Thomas Moerman, Vân Anh Huynh-Thu, Hana Imrichova, Gert Hulselmans, Florian Rambow, Jean-Christophe Marine, Pierre Geurts, Jan Aerts, Joost van den Oord, Zeynep Kalender Atak, Jasper Wouters, Stein Aerts
- Corresponding author: Stein Aerts (stein.aerts@kuleuven.vib.be, KU Leuven / VIB Center for Brain & Disease Research)
- 雑誌: Nature Methods
- 発行年: 2017
- Epub日: 2017-10-09
- Article種別: Brief Communication (Methods)
- PMID: 28991892
背景
単細胞RNA-seq (scRNA-seq) 技術の急速な進歩により、個々の細胞の転写状態と分化過程の高解像度同定が可能になった (Linnarsson et al. 2016; Stegle et al. 2015)。しかし細胞状態の根底にある GRN (gene regulatory network; 遺伝子調節ネットワーク) — TF (transcription factor; 転写因子) が互いや下流遺伝子を制御する調節回路 — を単細胞解像度で推定する手法は根本的に不足していた。既存の計算法は機械学習による共発現パターンの推定にとどまり (Moignard et al. 2015; Pina et al. 2015)、シングルセルレベルの転写バーストに起因する確率的遺伝子発現変動 (Raj 2008) を反映した大量の擬陽性間接標的が混入する問題があった。TFとその標的遺伝子間のcis調節配列 (シス調節; プロモーターや近位エンハンサーへのTF結合) を利用して直接的な調節関係のみを抽出する手法は未解明な課題として残されており (Trapnell et al. NatBiotechnol 2014)、腫瘍データでは腫瘍固有のバッチ効果が細胞状態解析をさらに困難にしていた (Tirosh et al. Science 2016)。GRNに基づいて単細胞データから細胞状態を同定・クラスタリングする統合的パイプラインの不足が、分野全体の瓶頸となっていた。
目的
SCENIC (Single-Cell Expression Network Inference and Clustering; 単細胞発現ネットワーク推定とクラスタリング) を開発し、(1) scRNA-seqデータからTFとその直接標的遺伝子群「レギュロン (regulon)」を高精度で同定し、(2) 各細胞でのレギュロン活性に基づいて腫瘍横断的な細胞状態をクラスタリングする統合的計算フレームワークを構築・検証すること。
結果
SCENICアルゴリズム — GRN推定・cis調節フィルタ・AUCスコアリングの3ステップパイプライン：
SCENICは3つの連続ステップで構成される (Fig 1a)。ステップ1では GENIE3 (gene network inference with ensemble of trees; ランダムフォレスト回帰に基づくGRN推定ツール) が各TFと候補標的遺伝子間の共発現関係を大規模にスコアリングし共発現モジュールを生成する。ステップ2では RcisTarget (Regulatory cis-target enrichment tool; シス調節モチーフエンリッチメント解析パッケージ) が各共発現モジュールについてシス調節配列のモチーフエンリッチメント (FDR p < 0.05) を実施し、対応するTFの結合モチーフが有意に濃縮されたモジュールのみを保持・精製してレギュロン (regulon; 転写因子とその直接cis調節標的遺伝子群) を定義する。ステップ3では AUCell (area under curve-based cell scoring; 各細胞のレギュロン活性AUCスコアリングアルゴリズム) が各細胞においてレギュロンを構成する遺伝子のランキングから AUC (area under the curve; 曲線下面積) スコアを算出し、バイナリ活性マトリクスに変換する。このマトリクスから細胞状態のクラスタリングとGRN活性の可視化が行われる。SCENICはhttp://scenic.aertslab.orgで公開されている。
マウス脳での検証 — 感度・特異度・ARI指標で既存法を上回るクラスタリング精度：
成体マウス脳scRNA-seqデータ (Zeisel et al. Science 2015) に適用したところ、初期1,046共発現モジュールから151レギュロン (7% of initial TFs) が同定され (Fig 1b)、 ARI (adjusted Rand index; 調整ランド指数) >0.80、sensitivity 0.88、specificity 0.99 を達成し、専用単細胞クラスタリング手法を上回った (Fig 1e)。バイナリレギュロン活性マトリクスの t-SNE (t-distributed stochastic neighbor embedding; t分布型確率的近傍埋め込み) では、microglia・astrocyte・interneuron・pyramidal neuronなど主要細胞型が明確に分離された (Fig 1d)。各細胞型のマスター制御因子 (例: microglia — Cebpb/Stat6/Nfkb1; interneuron — Dlx1/Dlx2) が文献と一致して同定された (Fig 1c)。100細胞サブサンプリングや1/3シーケンシングカバレッジでも2-6細胞しかいない稀少microglia・astrocyte・interneuronが堅牢に検出された。さらにヒト脳データ (Lake et al. Science 2016) との比較では、DLX1 (distal-less homeobox 1) / DLX2 regulonのinterneuronにおける活性が種間で保存されることを確認し (Fig 2a, b)、マウス+ヒト合同クラスタリングでも種ドリブンではなく細胞型ドリブンの分離が得られた (Fig 2c)。
腫瘍データへの応用 — バッチ効果を克服した腫瘍横断的細胞状態の同定：
腫瘍scRNA-seqデータではSCENICの特性が最も際立った。乏突起膠腫データ (Tirosh et al. Nature 2016) では n=4043 cells from 6 tumors を解析し (Fig 3a-e)、通常の発現クラスタリングが腫瘍起源別に細胞を分離するのに対し、SCENICのバイナリレギュロン活性マトリクスは3腫瘍横断的癌細胞状態 — oligodendrocyte様 (SOX10 (SRY-box transcription factor 10) / SOX4/SOX8/OLIG1 (oligodendrocyte lineage-specific factor 1) /OLIG2/ASCL1駆動)、astrocyte様 (SOX9/NFIB/AP-1駆動)、cycling cells (E2F/FOXM1駆動) — を明確に識別した。黒色腫データ (n=1252 cells from 14 lesions) では MITF^high (増殖性; MITF/STAT/IRF制御) とMITF^low (浸潤性; WNT5A/LOXL2/ZEB1高発現) の2大細胞状態に加え、MITF^low 状態に特異的な新規TFとして NFATC2 (Nuclear Factor of Activated T-Cells 2; 114 predicted target genes) と NFIB (Nuclear Factor I-B; 15 predicted target genes) を発見した (Fig 3f-h)。SCENICは Combat や Limma など専用バッチ補正手法と異なりバッチ効果の事前指定を必要とせず、生物学的GRN特徴量に基づいて腫瘍効果を自動除去する。
新規調節因子NFATC2/NFIBの実験的検証 — IHCとsiRNAノックダウン：
計算的に発見したNFATCおよびNFIBの臨床的意義を検証した (Fig 3i, j)。n=25 specimens の human melanoma IHC 解析では、NFIB/NFATC2の発現がsentinel lymph nodeで最高値を示し、浸潤マーカーZEB1との共局在が観察された。この所見は転移初期イベントとの関連を示唆する。melanoma A375 cell line (NFATC2/NFIB高発現) を用いたNFATC2 siRNAノックダウン実験では、NFATC2レギュロン遺伝子群の有意な上方制御と、黒色腫浸潤状態シグネチャ遺伝子・細胞接着/細胞外マトリクス調節遺伝子の活性化が確認された (Supplementary Table 1)。さらにMITFおよびSTAT1の予測標的は ChIP-seq データと高い一致を示した (Fig 3j)。スケールアップとして >40000 single cells のマウス網膜データでのサブサンプリング戦略と、GENIE3のScala/Apache Spark実装 GRNBoost (Gene Regulatory Network Boost; 勾配ブースティングによる高速GRN推定) による計算時間の劇的短縮も実証された (Bray et al. NatBiotechnol 2016 による高速RNA-seq定量とともに大規模scRNA-seq解析基盤を支える)。
考察/結論
① 先行研究との違い: 共発現ネットワーク推定のみを行う従来GRN手法 (GENIE3単独、Moignard et al. 2015) とは異なり、SCENICはRcisTargetによるcis調節配列エンリッチメントで間接標的を系統的に除去しレギュロンの特異度を劇的に向上させる。発現正規化ベースのクラスタリング (Seurat: Cell 2019-Stuart) と対照的に、SCENICはレギュロン活性という生物学的に意味ある特徴量を用いるため、腫瘍起源バッチ効果を事前指定なしに自動除去し腫瘍横断的な真の細胞状態を抽出できる点が根本的に相違する。
② 新規性: cis調節配列解析とscRNA-seqデータを統合し、バイナリレギュロン活性マトリクスによる細胞状態クラスタリングというアプローチは本研究で初めて提案された。腫瘍データでのバッチ効果克服と、NFATC2/NFIBという新規なメラノーマ浸潤性転写調節因子を計算的に発見してIHC・ノックダウンで実験的に検証した点は新規な貢献である。
③ 臨床応用: 腫瘍横断的細胞状態の客観的同定はがん免疫療法の治療応答予測、サブタイプ分類、薬剤標的同定に直接利用可能である。NFATC2/NFIBはメラノーマ転移初期段階の新規バイオマーカーとして臨床的意義があり、JNK/MAPKシグナルを介した転移制御機構の標的としての臨床応用が期待される。腫瘍内不均一性の客観的スコアリング指標として臨床現場への応用も見込まれる。
④ 残された課題: SCENICは既知TF結合モチーフを前提とするため、モチーフ未定義のオーファンTFには適用困難であり、今後の研究でモチーフデータベースの拡充と動的ネットワーク推定への拡張が求められる。Human Cell Atlas など百万細胞スケールへのスケールアップ、単細胞マルチオミクス (ATACseq+RNAseq) との統合、時系列scRNA-seqを用いた動的GRN推定など今後の検討課題は多く、limitation を含むさらなる検証が必要である。
方法
scRNA-seqデータ前処理後にSCENIC 3ステップを実行: (1) GENIE3: ランダムフォレスト回帰でTF-遺伝子共発現スコアを計算; (2) RcisTarget: cis調節配列のゲノムランキングとTF結合モチーフライブラリを用いてFDR p < 0.05でモチーフエンリッチメントを評価し間接標的を除去; (3) AUCell: 各細胞でレギュロン遺伝子AUCを算出しバイナリ活性マトリクスを生成。解析データ: 成体マウス脳 (Zeisel et al. Science 347:1138, 2015)、ヒト脳 (Lake et al. Science 352:1586, 2016)、乏突起膠腫 (Tirosh et al. Nature 539:309, 2016)、黒色腫 (Tirosh et al. Science 352:189, 2016)、マウス網膜 (>40000細胞, Macosko et al.)。IHC: n=25 human melanoma specimens、抗NFATC2/NFIB/ZEB1/EPHA2抗体。siRNA knockdown: NFATC2 siRNA in A375 melanoma cell line、下流遺伝子発現プロファイリングで検証。ChIP-seq検証: MITF/STAT1 ChIP-seqシグナルの予測標的領域と対照ランダム領域への集積を比較 (Fig 3j)。統計: Fisher正確検定ベースのモチーフエンリッチメント、AUCellスコアのKS検定でOn/Off閾値決定。GRNBoost: GENIE3のScala実装でgradient boosting使用、Apache Sparkで並列化。