Research Wiki

Integrating single-cell transcriptomic data across different conditions, technologies, and species

異なる実験条件、測定技術、あるいは生物種間で得られた単一細胞RNAシーケンシングデータを、正準相関分析と動的時間伸縮法を用いて共通の低次元空間に統合するアライメント手法を開発した。

20 min read

  • 著者: Andrew Butler, Paul Hoffman, Peter Smibert, Efthymia Papalexi, Rahul Satija
  • Corresponding author: Rahul Satija (New York Genome Center, New York, NY, USA; Center for Genomics and Systems Biology, New York University, New York, NY, USA)
  • 雑誌: Nature Biotechnology
  • 発行年: 2018
  • Epub日: 2018-04-02
  • Article種別: Original Article (Method)
  • PMID: 29608179

背景

scRNA-seq (single-cell RNA sequencing: 一細胞RNAシーケンシング) 技術は、複雑な組織内に存在する細胞のヘテロ性を一細胞レベルで高解像度に解明するための強力なツールとして急速に普及した。特に、Zheng et al. NatCommun 2017 らの報告にあるように、液滴ベースのハイスループットなシーケンシング技術の登場により、数万個規模の細胞プロファイリングが日常的に行われるようになっている。しかし、異なる実験条件 (刺激 vs 対照、健常 vs 疾患)、異なる測定技術 (Smart-seq2 vs 10x Chromium vs Drop-seq)、あるいは異なる生物種 (ヒト vs マウス) から得られたデータを横断的に統合して解析することは、依然として極めて困難な計算科学的課題であった。単純にデータを結合 (merge) するだけでは、技術的なバッチ効果 (batch effect) が生物学的なシグナルを完全に覆い隠してしまい、次元削減空間上で技術的差異が主たる変動軸となってしまう。

従来の bulk RNA-seq 向けに開発されたバッチ補正手法である ComBat や limma::removeBatchEffect Ritchie et al. NucleicAcidsRes 2015 は、各バッチ間で細胞型の構成比率が同一であると仮定しているため、細胞構成が異なるシングルセルデータに適用すると深刻な歪みが生じるという課題があった。また、条件特異的な生物学的応答を技術的バッチと混同して除去してしまう危険性も高かった。さらに、MNN (mutual nearest neighbor: 相互最近傍) 法などの初期のシングルセル向けバッチ補正法も提案されていたが、治療や刺激による「条件特異的な細胞状態の変化」を正確に保持しつつ、異なるデータセット間で共通の細胞集団を整列 (alignment) させる手法は未確立であった。

先行研究において、Satija et al. NatBiotechnol 2015 らは一細胞の空間再構成技術を開発し、Tirosh et al. Science 2016 らは腫瘍微小環境における細胞生態系の解明を進めてきたが、条件や技術、種を横断してデータを統合する一般的なアライメント戦略は不足していた。特に、異なるデータセット間で「共通する生物学的変動軸」を同定し、非線形な細胞密度の変化や遺伝子発現スケールの差異を補正する包括的なアプローチは未解明であり、大規模な統合解析の大きな障壁となっていた。このように、異種データを頑健に統合するための数理モデルやアルゴリズムが決定的に不足しているという gap が残されていた。

目的

本研究の目的は、画像処理分野における「多様体アライメント (manifold alignment)」の概念を scRNA-seq データ解析に応用し、異なる条件、測定技術、あるいは生物種間で得られた複数のデータセットを、共通 of 低次元空間に統合する新しい計算手法を開発することである。具体的には、複数データセット間で共有される遺伝子相関構造を CCA (canonical correlation analysis: 正準相関分析) によって同定し、DTW (dynamic time warping: 動的時間伸縮法) を用いて非線形に整列するアルゴリズムを構築する。

これにより、条件特異的な生物学的応答 (例: 薬物や刺激に対する反応) を保持したまま、共通の細胞亜集団 (subpopulation) を高精度に同定することを可能にする。開発した手法を R パッケージ Seurat v2 に実装・公開し、(i) 静止状態 vs IFN-β (interferon-beta: インターフェロンβ) 刺激下のヒト PBMC (peripheral blood mononuclear cell: 末梢血単核細胞)、(ii) 異なる scRNA-seq 技術で測定されたマウス造血前駆細胞、(iii) ヒトおよびマウスの膵島アトラス、という3つの代表的な統合シナリオにおいて、その有効性と頑健性を実証することを目的とする。

結果

IFN-β 刺激による PBMC の条件横断的アライメント: 対照群と IFN-β 刺激群の PBMC データセット (計 n=14039 cells) を Seurat アライメントを用いて統合した結果、刺激の有無によるバッチ効果が完全に排除され、細胞型に基づいて共通の空間に整列された (Fig. 2b)。共同クラスタリングにより、CD4 記憶 T細胞、CD4 ナイーブ T細胞、CD8 T細胞、B細胞、NK (natural killer: ナチュラルキラー) 細胞、CD14 陽性単球、CD16 陽性単球、cDC、pDC、巨核球、赤芽球など、計 13 個の aligned subpopulation を同定した (Fig. 2c)。対照的に、単純な結合 (naïve merge) では刺激の有無が主たるクラスタリング軸となり、細胞型の分離が著しく損なわれていた (Fig. 2a)。アライメント後の各クラスタ内において、IFN-β 刺激による遺伝子発現変化を解析したところ、すべての細胞型で共通して発現上昇する ISG15 (interferon-stimulated gene 15) や IFIT1 (interferon-induced protein with tetratricopeptide repeats 1) などの一様マーカーが同定された一方、単球特異的に 1.98-fold の発現低下を示す CD14 など、細胞型特異的な応答プログラムが明らかになった (Fig. 2d)。さらに、pDC は他の骨髄系・リンパ系細胞とは異なる極めてユニークな IFN-β 応答シグネチャーを示すことが判明し、この所見はソートした細胞を用いた n=3 replicates のバルク RNA-seq 検証実験によって強く支持された (Fig. 2g)。

異なる scRNA-seq 技術間での造血前駆細胞データの統合: 測定技術、シーケンス深度、および UMI (unique molecular identifier) の dropout 率が大きく異なる造血前駆細胞データセット (MARS-seq: n=2686 cells、SMART-Seq2: n=765 cells) を Seurat アライメントにより統合した (Fig. 3b)。アライメント後、SMART-Seq2 の各細胞を MARS-seq の 18 個の既知クラスタにマッピングした結果、両技術間で系統特異的マーカー (GATA2、APOE、FOS など) の発現パターンが極めて高い精度で一致した (Fig. 3e, f)。さらに、赤血球分化経路にコミットした細胞群に対して拡散マップ (diffusion map) を適用し、共同発達軌跡 (joint developmental trajectory) を再構成した (Fig. 3g)。この軌跡に沿った擬時間 (pseudotime) 解析により、両技術間で保存された遺伝子発現動態が明らかになった一方、SMART-Seq2 において細胞ストレスに起因する JUN や FOS の特異的な 2.5-fold の発現上昇などの技術特異的効果も高解像度で検出された (Fig. 3h, i)。

種を超えたヒトとマウスの膵島アトラスの統合と希少細胞の同定: 数百万年の進化の歴史を経てグローバルな転写プロファイルが大きく乖離しているヒト (n=8424 cells) とマウス (n=1767 cells) の膵島細胞データを、1:1 のオーソログ遺伝子セットを用いて統合した (Fig. 4b)。Seurat アライメントは、種間の発現差異を乗り越えて、α細胞、β細胞、δ細胞、γ細胞、導管細胞、腺房細胞、星細胞、内皮細胞、免疫細胞などの主要な膵島細胞型を正確に対応付けた (Fig. 4c)。特筆すべき所見として、共同クラスタリングにより、ヒト (n=126 cells、β細胞全体の約 5.2%) およびマウス (n=10 cells、β細胞全体の約 1.3%) の両方において、ER (endoplasmic reticulum: 小胞体) ストレス応答遺伝子 (HERPUD1、GADD45A、ATF3、DDIT3) を特異的に高発現する希少な β細胞亜集団 (Beta_er_stress) が新規に同定された (Fig. 4f)。この亜集団の発現遺伝子は、GO (Gene Ontology) 解析において「ER unfolded protein response」などの生物学的プロセスに極めて強く富化していた (Fig. 4g)。

既存のバッチ補正手法との性能比較ベンチマーク: Seurat アライメントの性能を、広く用いられている ComBat および limma::removeBatchEffect と比較評価するため、近傍細胞のデータセット混合度を測定する「アライメントスコア」を算出した (Fig. 5)。PBMC、造血前駆細胞、および膵島データのすべてのシナリオにおいて、Seurat アライメントは ComBat や limma を大幅に上回るアライメントスコア (Seurat: >0.80 vs ComBat: <0.40, limma: <0.30) を達成した (Fig. 5c, f, i)。ComBat や limma では、細胞構成の異なるデータセット間で細胞の分布が著しく歪む現象が観察されたが、Seurat は細胞密度のシフトに対しても極めて頑健に動作した。さらに、陰性対照として異なる組織由来のデータセットをアライメントした際には、アライメントスコアが著しく低下し、生物学的な妥当性のない無理な整列が自動的に回避されることも確認された。

考察/結論

先行研究との違い: 本研究で開発された Seurat アライメント手法は、従来の bulk RNA-seq 向けの線形バッチ補正手法 (ComBat や limma) と異なり、各データセットの細胞型構成が不均一である場合や、細胞密度に大幅なシフトが存在する場合でも、生物学的なシグナルを歪めることなく頑健にデータを統合できる。また、MNN などの初期のシングルセル統合手法と比べても、条件特異的な活性化状態や刺激応答 (例: IFN-β 応答) を技術的バッチ効果と混同して消去することなく、正確に保持したまま整列できる点で決定的に優れている。

新規性: 本研究は、画像処理分野の多様体アライメント技術をシングルセルゲノミクスに初めて導入し、正準相関分析 (CCA) と動的時間伸縮法 (DTW) を組み合わせた非線形アライメントアルゴリズムを新規に開発した。これにより、測定技術 (Smart-seq2 vs MARS-seq) や生物種 (ヒト vs マウス) の壁を越えて、保存された遺伝子共発現ネットワークを抽出し、共通の低次元空間に投影することが可能となった。特に、種を超えて保存された小胞体 (ER) ストレス応答を示す希少な β細胞亜集団を本研究で初めて自動的に同定したことは、本手法の極めて高い検出能と新規性を示している。

臨床応用: 本手法は、トランスレーショナルリサーチ (translational research) や臨床医学研究において極めて高い臨床的有用性を持つ。例えば、がん免疫療法における治療前後の腫瘍微小環境の比較、健常者と疾患患者のコホートデータの統合、あるいはヒト疾患の動物モデル (マウスなど) から得られたシングルセルデータの臨床応用への橋渡し (bench-to-bedside) において、共通の治療標的細胞や希少な病態関連細胞 (例: ERストレス下にあるβ細胞や特定の活性化T細胞亜集団) を高精度に同定するための標準的プラットフォームとなる。

残された課題: 今後の検討課題として、数十から数百に及ぶ大規模なマルチバッチデータセットや、数百万細胞スケールの大規模アトラスデータに対する計算スケーラビリティの向上が挙げられる。また、本手法の limitation として、データセット間で完全に非重複な細胞集団 (片方のデータセットにしか存在しない独自の細胞型) が存在する場合、それらを完全に分離しつつ共通部分だけを整列させるための閾値設定 ( の制御など) には依然として改善の余地が残されている。今後の研究方向性として、空間トランスクリプトミクスデータやシングルセルエピゲノムデータ (scATAC-seq) など、異種オミクスデータの統合 (multi-omics integration) への拡張が期待される。

方法

本研究で開発された Seurat アライメントワークフローは、以下のステップで構成される。

  1. データの前処理と高変動遺伝子の選択: 各データセットに対して独立に対数標準化 (log-normalization) を行い、各遺伝子の分散・平均比 (dispersion) を算出する。各データセットから上位 1,000 個の HVG (highly variable gene: 高変動遺伝子) を選択し、それらの論理和 (union) を CCA の入力遺伝子セットとする。
  2. CCA による共有相関空間の構築: 選択された遺伝子セットを用いて、2つのデータセット と の間で最大相関を持つ CCV (canonical correlation vector: 正準相関ベクトル) を算出する。本手法では、細胞数が遺伝子数よりも大幅に多い高次元データに対応するため、共分散行列を対角行列と仮定する対角 CCA (diagonal CCA) を採用し、irlba パッケージを用いた部分 SVD (singular value decomposition: 特異値分解) により高速に計算を行う。
  3. DTW による CCV の非線形アライメント: 各 CCV について、相関の高い上位 30 個 of 遺伝子の加重平均として「メタ遺伝子 (metagene)」を定義する。メタ遺伝子の 95% 参照範囲を一致させる線形変換を行った後、DTW アルゴリズムを適用して非線形な歪み (細胞密度のシフトなど) を補正し、両データセットの細胞を共通の低次元空間に整列させる。
  4. 非重複集団の同定: CCA 空間が各細胞の分散をどの程度説明できるかを、独立に実行した PCA (principal component analysis: 主成分分析) の説明分散と比較し、その比率 が 0.5 未満の細胞を「非重複集団 (non-overlapping population)」として自動的に検出・除外する。
  5. 統合空間での下流解析: 整列された CCV 空間において、SNN (shared nearest neighbor) グラフを構築し、SLM (smart local moving) アルゴリズムによるモジュラリティベースのクラスタリングを実行する。可視化には t-SNE (t-distributed stochastic neighbor embedding) を用いる。
  6. 検証用データセット:
    • PBMC データセット: 8人の健常ドナー由来の PBMC を対照群と IFN-β 刺激群に分割し、10x Chromium 技術で測定したデータ (計 n=14039 cells)。
    • 造血前駆細胞データセット: C57BL/6 マウスの骨髄から採取した造血前駆細胞を、SMART-Seq2 (Switching Mechanism at 5’ end of RNA Template Sequencing 2) 技術 (n=765 cells) および MARS-seq (Massively Parallel RNA Single-Cell Sequencing) 技術 (n=2686 cells) の異なる技術で測定したデータ。
    • 膵島データセット: ヒト (n=8424 cells) およびマウス (n=1767 cells) の膵島細胞を inDrops 技術で測定したデータ。
  7. 統計解析および検証実験: クラスタ間の差分発現遺伝子の同定には Wilcoxon rank-sum test を使用し、有意水準は とした。また、メタ遺伝子の発現レベルの要約統計量として mean ± SD を算出した。さらに、pDC (plasmacytoid dendritic cell: プラズマサイトイド樹状細胞) および cDC (conventional dendritic cell: 通常樹状細胞) をフローサイトメトリーでソートし、IFN-β 刺激の有無によるバルク RNA-seq 検証実験 (n=3 replicates) をトリプリケートで実施した。

グラフビュー

バックリンク