- 著者: Ilya Korsunsky, Nghia Millard, Jean Fan, Kamil Slowikowski, Fan Zhang, Kevin Wei, Yuriy Baglaenko, Michael Brenner, Po-ru Loh, Soumya Raychaudhuri
- Corresponding author: Soumya Raychaudhuri (Brigham and Women’s Hospital / Broad Institute / University of Manchester)
- 雑誌: Nature Methods
- 発行年: 2019
- Epub日: 2019-11-18
- Article種別: Original Article
- PMID: 31740819
背景
single-cell RNA-seq (scRNA-seq) 技術の急速な進歩により、数千から数百万細胞の網羅的転写プロファイリングが可能となり、Human Cell Atlas (HCA) や Accelerating Medicines Partnership (AMP) のような大規模な参照データセットが急増している。これらのデータセットは、多様な生物学的・臨床的条件下での細胞タイプの包括的な特徴付けを可能にする。しかし、異なるドナー、研究、およびプラットフォーム間で取得された複数のデータセットを統合して解析しようとすると、技術的なバッチ効果が生物学的シグナルと交絡し、共通の細胞型を正確に同定することが困難になるという課題が浮上している。この問題に対処するため、これまでにもMNN Correct、BBKNN、MultiCCA/Seurat3、Scanoramaといった複数のデータセット統合アルゴリズムが提案されてきた。例えば、Butler et al. NatBiotechnol 2018 はSeurat3を開発し、Zheng et al. NatCommun 2017 は大規模な単一細胞プロファイリング技術を報告している。また、Azizi et al. Cell 2018 は乳腺腫瘍微小環境における免疫細胞の多様な表現型を単一細胞レベルで解析した。しかし、これらの先行手法は、(i) 大規模データセットへのスケーラビリティ、(ii) 広範な細胞集団と微細なサブ集団の同時同定、(iii) 複数バッチ、組織、ドナーといった複雑な実験デザインへの柔軟な対応、(iv) scRNA-seqと空間トランスクリプトミクスのような異なるモダリティを跨ぐ統合、という4つの主要な課題を同時に解決するには至っていなかった。特に、大規模な参照アトラスの構築には、計算効率と精度を両立する新たな統合手法が不可欠であり、既存のツールではメモリや計算時間のボトルネックが顕著に現れることが報告されている。このような背景から、技術的差異を効率的に補正しつつ、細胞タイプ間の生物学的差異を維持できる、より高性能なデータ統合ツールの開発が強く求められていたが、そのための効果的なアルゴリズムは未解明な点が多かった。既存の手法では、これらの複雑な統合シナリオに十分に対応できる柔軟性が不足しており、特に大規模データセットにおける計算資源の制約が大きな課題として残されていた。
目的
本研究の目的は、上記背景で述べた4つの課題（スケーラビリティ、広範・微細な細胞集団の同定、複雑な実験デザインへの柔軟性、モダリティ間統合）を同時に解決する、新規のscRNA-seqデータ統合アルゴリズムHarmonyを開発することである。Harmonyは、PCAなどの低次元埋め込みを起点とし、ソフトクラスタリングに基づく細胞特異的な線形補正関数を反復的に学習する手法である。本研究では、Harmonyのアルゴリズムを詳細に記述し、その性能を既存の主要な統合手法（MNN Correct、BBKNN、MultiCCA、Scanorama）と比較することで、統合精度と計算効率の両面における優位性を実証することを目指す。また、Harmonyが個人用コンピュータ上で最大100万細胞規模のデータ統合を可能にすることを示し、Human Cell Atlasのような大規模プロジェクトへの適用可能性を提示する。
結果
アルゴリズム概要とLISI評価指標の導入: Harmonyは、PCAで得られた低次元空間において、細胞型・状態の連続的な遷移を保持しつつ、データセット間の技術的差異を効率的に補正する反復ソフトクラスタリングアルゴリズムである。既存のMNN Correct (Mutual Nearest Neighbors Correct) やMultiCCA (Multiple Canonical Correlation Analysis) とは異なり、Harmonyはソフトクラスタリングを用いて細胞型固有の線形補正を学習し、計算はPC空間上の線形操作のみで完結するため軽量である。本研究では、統合性能と精度を定量的に評価するため、新たにlocal inverse Simpson’s Index (LISI) を導入した。iLISIは局所近傍におけるデータセットの混合度を、cLISIは細胞型の純度をそれぞれ示し、両指標を組み合わせることで「過剰混合」や「不十分混合」を区別できるようになった。
Cell lineベンチマーク (Jurkat/293T): 3つのデータセット（純粋Jurkat細胞 n=3,255 cells、純粋293T細胞 n=2,859 cells、50/50混合Jurkat細胞 n=1,799 cellsと293T細胞 n=1,565 cells）を用いた評価では、Harmony適用後、median iLISI 1.59 (95% CI 1.27–1.97) を達成し、理論的な最大値 (Jurkat細胞で1.8、293T細胞で1.5) に近接する高いデータセット混合度を示した。同時に、median cLISI 1.00 (95% CI 1.00–1.02) を達成し、Jurkat細胞と293T細胞の正確な分離が維持されていることを確認した。これは、既存のMNN Correct、BBKNN (Batch-Balanced K-Nearest Neighbors)、MultiCCA、Scanoramaといったアルゴリズムと比較して、統計的に優位な統合性能であった (Supplementary Results, Fig 1b, Table 1)。
HCA大規模データでのスケーラビリティ: Human Cell Atlas (HCA) の528,688細胞データセットを30,000細胞から500,000細胞にダウンサンプリングして計算効率を比較した。Harmonyは、30,000細胞で4分・0.9 GB、500,000細胞で68分・7.2 GBと、細胞数に対してほぼ線形に近いスケーリングを示した (Fig 3a,b)。これはMultiCCAやMNN Correctと比較して30〜200倍高速であり、メモリ使用量も125,000細胞時点で他手法より30〜50倍少なかった。MNN Correct、MultiCCA、Scanoramaは125,000細胞を超えると実行不可能となるか、過剰なメモリを要求した。Harmonyは10⁶細胞オーダーのデータセットを個人用コンピュータ上で処理可能であることを実証した。
PBMCマルチテクノロジー統合: 10x Genomics v2/v3、Smart-seq2、CEL-Seq2など異なる技術で取得されたPBMC (Peripheral Blood Mononuclear Cell) データセットの統合において、Harmonyは技術差を効果的に補正しつつ、T細胞、B細胞、NK細胞、単球といった主要な細胞型構造を保持した (Fig 4a,b)。統合後、median iLISIは1.96 (95% CI 1.36–2.56) と高い混合度を示し、他の手法 (median iLISI 1.02, 95% CI 1.00–2.11) を上回った。また、HarmonyはCD4陽性ナイーブT細胞、エフェクターメモリーCD4陽性T細胞、Treg、メモリーCD8陽性T細胞、エフェクターCD8陽性T細胞、ナイーブB細胞、メモリーB細胞といった微細なサブタイプも、技術的な交絡なく同定可能であった (Fig 4c,d)。
膵島メタ解析 (5研究): 5つの独立した研究から得られた36ドナーのヒト膵島scRNA-seqデータ (n=14,746 cells) を統合した。Harmonyは、異なる技術プラットフォームと複数ドナーの両方のバッチ効果を同時に補正し、α細胞、β細胞、δ細胞、PP細胞、外分泌細胞などの主要な細胞型を一貫して同定した (Fig 5a-d)。Harmonyは、技術 (median iLISI 2.17, 95% CI 1.02–3.91) とドナー (median iLISI 5.05, 95% CI 1.24–10.05) の両方で高い統合度を達成し、細胞型分類精度は98%以上であった。他のアルゴリズム（BBKNN, Scanorama, MNN Correct, MultiCCA, Ritchie et al. NucleicAcidsRes 2015）では、ドナーと技術の両方を同時に効果的に混合することはできなかった。さらに、Harmonyは、ER (Endoplasmic Reticulum) ストレス下のβ細胞や、これまで報告されていないERストレス下のα細胞といった稀な細胞サブタイプ（全細胞の2%未満）を、技術的交絡なく抽出することに成功した (Fig 5e-j)。これらのERストレス細胞集団の割合はドナー間で有意に相関しており (Spearman r=0.46, p=0.004)、糖尿病におけるα細胞損傷とβ細胞機能不全の並行性を示唆した (Fig 5k)。
マウス胚発生 (縦断): マウスの造血発生における8つの異なる時点 (E6.75からE8.5) からのn=15,875 cellsを統合した。Harmonyは、各サンプルが異なる発生段階に由来し、全ての細胞型が全てのサンプルに存在しないという課題にもかかわらず、バッチ効果を除去しつつ、発生軌跡の連続的な細胞状態遷移を保持した (Supplementary Fig 16a-e)。Harmonyのソフトクラスタリングは、共通の前駆細胞から内皮系および造血系への分岐といった、連続的な構造を維持する上で重要であった。
Cross-modality統合 (scRNA-seq + spatial transcriptomics): マウス視床下部視索前野のscRNA-seqデータ (10x Genomics, n=30,370 cells) と空間トランスクリプトミクス (MERFISH, n=64,373 cells) データを、154の共通遺伝子を用いて統合した。Harmonyは、限られた共通遺伝子数と検出効率の大きな違いにもかかわらず、両モダリティの細胞を共通の埋め込み空間に混合し (median iLISI 1.15, 95% CI 1.00–1.99)、12の主要な神経細胞およびグリア細胞集団を正確に統合した (予測精度93.7%) (Fig 6a,b)。これにより、MERFISHで測定されていない遺伝子の空間的局在を推定することが可能となり、これまで報告されていない神経転写因子の空間パターンを特定した。例えば、クロマチンオーガナイザーであるSatb1は、抑制性ニューロンにおいて強く空間的自己相関を示し (Moran’s I = 0.44)、その予測発現パターンはAllen Brain Atlasの画像と定性的に一致した。
考察/結論
先行研究との違い: Harmonyは、単一細胞RNAシーケンスデータ統合における主要な課題を解決する、高速かつ高精度なアルゴリズムである。先行研究であるMNN CorrectやMultiCCAといった既存手法と異なり、HarmonyはPCA等の低次元埋め込みを起点に、ソフトクラスタリングを用いて細胞型固有の線形補正を反復的に学習する点で、既存手法と異なる。特に、情報量ペナルティを導入することで、複数の技術的・生物学的因子を同時に補正できる柔軟性を持つ。また、PC空間上の線形操作のみで完結する軽量設計により、既存手法が抱えていた大規模データにおけるメモリボトルネックを克服した。
新規性: 本研究で初めて、Harmonyが6つの多様なベンチマーク解析において、既存アルゴリズムと比較して統合精度（iLISI中央値 1.59 vs 1.01）と計算効率（500,000細胞でHarmonyは68分、7.2GBメモリに対し、他は125,000細胞でメモリ超過）の両方で優れていることを初めて実証した。これにより、個人用PCで最大100万細胞の統合が可能となる。また、LISIという新規評価指標の導入により、統合の質をより客観的に評価できるようになった点も重要である。さらに、膵島細胞のメタ解析では、これまで報告されていないERストレス下のα細胞集団を同定し、糖尿病におけるα細胞損傷とβ細胞機能不全の並行性を示唆するなど、生物学的な新規知見も得られた。
臨床応用: Harmonyの計算効率と柔軟性は、Human Cell AtlasやAccelerating Medicines Partnershipのような大規模な参照アトラス構築に不可欠である。複数施設・複数バッチの患者検体を統合することで、疾患特異的な細胞集団や稀なサブタイプを同定し、疾患メカニズムの解明や新たな治療標的の発見に繋がる可能性がある。例えば、関節リウマチやループス腎炎のような炎症性疾患において、多様なドナーや組織から得られた免疫細胞のデータを統合することで、疾患関連細胞の同定が加速されることが期待される。また、異なるモダリティ（scRNA-seqと空間トランスクリプトミクス）の統合能力は、細胞の転写プロファイルと空間的文脈を結びつけ、疾患病態の理解を深めるための強力なツールとなる。
残された課題: 今後の検討課題としては、データセット間で細胞型構成が大きく異なる場合のペナルティ調整の最適化や、極端にサイズ差のあるデータセットにおけるLISI解釈の困難性が挙げられる。また、Harmonyは既知の変動源の影響をモデル化し除去するが、未知の潜在的変動源に対処するための拡張も今後の検討課題である。将来的には、Harmonyを遺伝子カウントの正確なモデリングに拡張し、La et al. Nature 2018のような全遺伝子発現プロファイルを必要とする手法の前処理に適用することや、数億細胞規模の参照アトラスへの高速マッピング機能の開発が期待される。Harmonyは、その優れた性能と柔軟性から、単一細胞解析の標準ツールとして広く採用され、Seuratパイプラインにも統合されている。
方法
Harmonyアルゴリズムは、PCAで得られた低次元座標から開始し、以下の4ステップを収束するまで反復する。まず、(a) ソフトk-meansクラスタリングを用いて、各細胞を複数のクラスターに確率的に割り当てる。この際、データセットの多様性を最大化する情報量ペナルティを課すことで、クラスターが単一のデータセットに偏ることを防ぐ。次に、(b) 各クラスターにおけるデータセット別のセントロイドとグローバルセントロイドを計算する。続いて、(c) これらのセントロイド間の差から、データセット別の線形補正因子を算出する。最後に、(d) 各細胞はソフトメンバーシップで重み付けされた線形和としてこれらの補正因子を受け取り、その細胞特異的な線形因子で座標が補正される。この4ステップは、細胞のクラスター割り当てが安定するまで繰り返される。
アルゴリズムの評価指標として、本研究では新たにlocal inverse Simpson’s Index (LISI) を定義した。LISIは、各細胞の局所近傍におけるデータセットの混合度を定量化するiLISI (integration LISI) と、細胞型の純度を定量化するcLISI (cell-type LISI) の2種類で構成される。iLISIは近傍内の有効データセット数を、cLISIは細胞型純度をそれぞれ示し、両指標で同時に良好な値を示すことで、「過剰混合」と「不十分混合」を区別できる。
Harmonyの性能は、以下の6つの異なるシナリオでベンチマークされた。(1) Jurkat/293T細胞株データセット: 3つのデータセット（純粋Jurkat細胞、純粋293T細胞、50/50混合細胞）を用いて、統合と細胞型保持の精度を評価した。 (2) 末梢血単核球 (PBMC) マルチテクノロジー統合: 10x Genomics v2/v3、Smart-seq2、CEL-Seq2など異なる技術で取得されたPBMCデータセットを統合し、技術差の補正能力と細胞型構造の保持能力を評価した。 (3) 膵島細胞メタ解析: 5つの独立した研究から得られたヒト膵島scRNA-seqデータを統合し、複数ドナー・複数研究にわたる細胞型の一貫した同定能力を評価した。 (4) Human Cell Atlas (HCA) 大規模データセット: 528,688細胞からなるHCAデータをダウンサンプリングし、30,000細胞から500,000細胞までの規模で、Harmonyの計算効率（ランタイムとメモリ使用量）とスケーラビリティを既存手法と比較した。 (5) マウス胚発生データセット: 複数時点・複数バッチをまたぐ発生軌跡データセットを統合し、時間軸に沿った連続的な細胞状態遷移を保持する能力を評価した。 (6) scRNA-seqと空間トランスクリプトミクス (MERFISH) のクロスモダリティ統合: 異なるモダリティ間で共通の埋め込み空間に細胞を射影し、空間情報と細胞型同定の橋渡しが可能であることを示した。
比較対象として、MNN Correct (scran Rパッケージ v1.9.4)、BBKNN (Teichlab/bbknn)、MultiCCA (Seurat Rパッケージ v2.3.4)、Scanorama (brianhie/scanorama) の各アルゴリズムが用いられた。HarmonyはRパッケージとして公開されており (github.com/immunogenomics/harmony)、Seurat (single-cell RNA sequencing analysis toolkit) パイプラインへの統合も可能である。LISIの計算はgithub.com/immunogenomics/LISIで提供されている。統計解析には、LISIスコアの比較にブートストラップ推定フレームワークを用い、差次的発現解析にはLIMMA Rパッケージの二側性調整t検定を、相関解析にはSpearman相関を用いた。