- 著者: Tim Stuart, Andrew Butler, Paul Hoffman, Christoph Hafemeister, Efthymia Papalexi, William M. Mauck III, Yuhan Hao, Marlon Stoeckius, Peter Smibert, Rahul Satija
- Corresponding author: Rahul Satija (New York Genome Center / NYU Center for Genomics and Systems Biology, New York, NY, USA)
- 雑誌: Cell
- 発行年: 2019
- Epub日: 2019-06-13
- Article種別: Original Article
- PMID: 31178118
背景
シングルセルゲノミクス技術の急速な進歩により、個々の細胞レベルでトランスクリプトームを網羅的に解析する scRNA-seq (single-cell RNA sequencing) は、細胞の多様性や状態を理解するための標準的な手法として定着した。先行研究である Zheng et al. NatCommun 2017 などの技術開発により、数万から数十万個の細胞を並行してプロファイリングすることが可能となった。さらに、近年ではクロマチンアクセシビリティを測定する scATAC-seq (single-cell assay for transposase-accessible chromatin using sequencing) や、細胞表面タンパク質を同時に定量する CITE-seq (cellular indexing of transcriptomes and epitopes by sequencing)、さらには組織内の位置情報を保持したまま遺伝子発現を測定する空間トランスクリプトーム技術など、多様なシングルセルモダリティが開発されている。
しかし、これらの異なる技術やモダリティから得られるデータは、それぞれ独自のバイアス、感度、ノイズ特性を有しており、単純に比較や統合を行うことは極めて困難であった。先行研究である Butler et al. NatBiotechnol 2018 で導入された CCA (canonical correlation analysis) や、独立に開発された MNN (mutual nearest neighbor) 法は、同一モダリティ内のバッチ効果補正において一定の成果を収めた。しかし、(1) 比較するデータセット間で細胞型の構成が異なる非対称な状況、(2) 技術的な差異が生物学的な差異を上回るケース、(3) トランスクリプトーム以外の異なるモダリティ（エピゲノム、プロテオーム、空間情報）を横断する統合、においては依然として有効なアプローチが未確立であった。このように、多様なシングルセルデータを調和させ、参照データからクエリデータへと情報をシームレスに転送するための包括的な計算フレームワークが決定的に不足していた。この技術的 gap を解消し、異種モダリティ間の関係性を明らかにする手法はいまだ未開拓であり、統合解析における大きな課題となっていた。
目的
本研究の目的は、異なるプラットフォームやバッチ間で得られた scRNA-seq データの調和（ハモナイゼーション）を達成するだけでなく、エピゲノム（scATAC-seq）、プロテオーム（CITE-seq）、空間トランスクリプトームなどの異種シングルセルモダリティを統一的に統合するための、新しい計算フレームワークを開発することである。具体的には、(1) 異なるデータセット間で共通の生物学的状態にある細胞ペアを「アンカー（anchors）」として同定するアルゴリズムの構築、(2) アンカーを介した複数データセットの共通空間への投影、(3) 参照データセット（reference）からクエリデータセット（query）への細胞型ラベルや連続値データの高精度な転送（transfer learning）、(4) 空間トランスクリプトームデータにおける全トランスクリプトーム情報の補完（imputation）を実現し、これらを R パッケージ Seurat v3 として実装・検証することを目指す。
結果
膵臓および骨髄データセットにおけるアンカーベース統合の精度検証:
膵臓の4つの異なるプラットフォーム（CelSeq、CelSeq2、Smart-Seq2、inDrop）から得られた n=6321 cells のデータセットを用いて、Seurat v3 のバッチ補正性能を検証した。一部の細胞型を意図的に特定のデータセットから排除した非対称な状況下においても、Seurat v3 は他手法（Seurat v2、MNN、Scanorama など）と比較して、共有されていない細胞集団の誤統合を有意に抑制した。シルエット係数を用いた評価において、Seurat v3 は細胞型の分離度とバッチ間の混合度の両面で最も優れたスコアを示した (Figure 2)。アンカーのスコアリングシステムにより、正しい細胞型ペアのアンカーは誤ったペアに比べて高いスコアを示し、これが誤統合の防止に寄与していることが確認された。さらに、網膜双極細胞の n=6 batches のデータセットを用いた検証でも、バッチ間の細胞混合度と局所構造の維持において高いロバスト性を示した。
scATAC-seqとscRNA-seqの統合による新規インターニューロンサブセットの同定:
マウス大脳皮質の n=3482 nuclei から得られた scATAC-seq データと、n=14249 cells の scRNA-seq 参照データを統合した。scATAC-seq の各遺伝子領域およびプロモーター領域のシグナルから算出した遺伝子活性マトリックス（gene activity matrix）を介してアンカーを同定し、転写産物ラベルを転送した。その結果、クロマチンアクセシビリティ単独の解析では分類困難であった稀少なインターニューロンサブセット（Vip-Mybpc1 (Vasoactive intestinal peptide - Myosin binding protein C, slow-type) や Pvalb-Vipr2 (Parvalbumin - Vasoactive intestinal peptide receptor 2) など）を高精度に同定することに成功した (Figure 3)。さらに、同定された PV (parvalbumin) 陽性細胞特異的なアクセシビリティピークにおいて、転写因子 Mef2c の結合モチーフが極めて高い有意差（p<10^-22）で濃縮されていることを発見した。これは、Heinz et al. MolCell 2010 等で示された系統決定因子の挙動と一致し、Mef2c が PV 細胞の運命決定において 2.5-fold 以上の活性化を誘導することを示唆する。
CITE-seqデータを用いたヒト骨髄アトラスへのタンパク質発現情報の転送:
ヒト骨髄細胞 n=33454 cells を用いた CITE-seq 実験により、25種類の細胞表面タンパク質とトランスクリプトームを同時にプロファイリングした。このデータを、HCA (Human Cell Atlas) が公開している n=274932 cells の巨大な骨髄 scRNA-seq データセットに統合し、タンパク質発現量を予測した。クロスバリデーションにおいて、25タンパク質中23個で予測値と実測値の間に高い相関（中央値 R=0.826、すなわち 82% 以上の相関度）を認めた。予測された CD69 発現に基づき、CD8+ メモリーT細胞においてインターフェロン（IFN）-g (interferon-gamma) 応答遺伝子群が log2FC 1.8 以上の高発現を示す活性化サブクラスを同定した (Figure 4)。この解析は、Love et al. GenomeBiol 2014 などの標準的な発現変動解析ツールと組み合わせることで、タンパク質レベルの表現型と遺伝子発現の相関を大規模コホートで明らかにする有用性を示している。
空間トランスクリプトームデータにおける全トランスクリプトーム情報の補完:
マウス体性感覚皮質の osmFISH (cyclic single-molecule fluorescence in situ hybridization) データ（33遺伝子、n=5300 spots）および STARmap データ（1020遺伝子、n=1539 cells）を、Allen Brain Atlas の scRNA-seq データ（n=14249 cells）と統合した。アンカーを介して、空間データに測定されていない数万遺伝子の発現パターンを補完（imputation）した。補完された発現パターンは、独立した高感度 FISH (fluorescence in situ hybridization) 実験データと高い相関を示し、層特異的マーカー（Cux2、Rorb、Foxp2など）の空間分布を正確に再現した (Figure 5)。遺伝子ダウンサンプリング検証により、わずか200から300個の遺伝子情報があれば、90% 以上の精度で全トランスクリプトームの空間補完が可能であることが示された。
考察/結論
先行研究との違い: 本研究で開発された Seurat v3 のアンカーフレームワークは、従来のバッチ効果補正ツール（Seurat v2 の CCA や MNN など）と異なり、データセット間で細胞型の構成が異なる非対称な状況下でも、誤統合を強力に抑制するスコアリングシステムを導入した。また、単一のモダリティ内での統合にとどまらず、エピゲノム、プロテオーム、空間トランスクリプトームといった異種モダリティ間を横断して情報を転送できる点で、これまでの手法と対照的である。
新規性:
本研究で初めて、相互最近傍法（MNN）と正準相関分析（CCA）を組み合わせた「アンカー」という概念を定義し、異なるシングルセル技術間で共通の生物学的状態を教師なしで同定する統一的アルゴリズムを確立した。これにより、CITE-seq から得られたタンパク質発現情報を n=274932 cells の scRNA-seq データへ転送することや、空間データにおける全トランスクリプトーム情報の補完を新規に実現した。
臨床応用: 本フレームワークは、がん微小環境（TME (tumor microenvironment)）の不均一性を解明するための臨床応用に極めて有用である。例えば、肺がんなどの胸部腫瘍学において、複数の患者コホートから得られた scRNA-seq データをバッチ効果を排除して統合し、腫瘍浸潤リンパ球（TIL (tumor-infiltrating lymphocyte)）の微細なサブセットを同定することが可能となる。さらに、空間トランスクリプトームデータと scRNA-seq アトラスを統合することで、がん細胞と免疫細胞の局所的な相互作用を単一細胞解像度で再構築し、免疫チェックポイント阻害剤の治療感受性予測などの臨床現場における意思決定を支援する。
残された課題: 今後の検討課題として、アンカーの同定が初期の CCA 投影に依存しているため、データセット間で共有される細胞型が極めて少ない、あるいは完全に異なる細胞型のみで構成されている場合には、過剰補正（over-correction）が生じるリスクが残されている。また、数百万細胞規模の超巨大データセットに対する計算コストの増大や、モダリティ間の因果関係（例えば、クロマチン構造の変化がどのように遺伝子発現を誘導するか）を直接推定できない点が limitation として挙げられる。後続の Seurat v4 では、これらの課題に対して weighted nearest neighbor (WNN (weighted nearest neighbor)) 法などのマルチモーダル統合アプローチが導入されることになる。
方法
データの標準化と高変動遺伝子の選択:
各データセットの log 正規化を行った後、遺伝子発現の平均と分散の関係を制御するために、Hafemeister et al. GenomeBiol 2019 に基づく分散安定化変換を適用した。具体的には、生カウントデータから loess (local regression) モデルを用いて期待される分散を推定し、標準化残差を算出することで、発現レベルに依存しない HVG (highly variable gene) を選択した。本手法の性能評価には、ヒト膵島細胞のほか、HEK293T 細胞株から得られたデータも使用した。
アンカー同定アルゴリズム: 2つのデータセット間で共通の生物学的状態を示す細胞ペア（アンカー）を同定するため、まず diagonalized CCA を用いて両データを共通の低次元空間に射影した。射影された正準相関ベクトルに対して L2 (L2 normalization) を施した後、MNN 探索を行い、相互に最近傍となる細胞ペアを抽出した。
アンカーのフィルタリングとスコアリング: 誤って同定されたアンカー（異なる細胞型間のペアなど）を排除するため、元の高次元データ空間において、アンカーを構成する細胞の近傍に相手方の細胞が位置するかを検証するフィルタリングステップを導入した。さらに、SNN (shared nearest neighbor) グラフの概念を応用し、各アンカーペアの周囲における相互近傍の重複度を算出してアンカーの信頼性を 0 から 1 の範囲でスコアリングした。
複数データセットの統合と情報転送:
複数データセットを統合する際は、データセット間のアンカー数に基づく類似度行列を作成し、hierarchical clustering (階層的クラスタリング) によりガイドツリー（guide tree）を構築した。このツリーに沿って、アンカーのスコアと距離に基づく重み付け補正ベクトルを適用し、再帰的にデータを補正した。ラベル転送においては、アンカーを介した重み付き多数決分類器を構築した。特徴量補完では、参照データの連続値（遺伝子発現やタンパク質発現）にアンカー重み行列を乗算することで、クエリデータの値を予測した。予測精度の評価には、Pearson correlation (ピアソン相関係数) を用いた。また、将来的な応用として、CRISPR-Cas9 を用いた遺伝子ノックアウトスクリーニングデータへの適用性も考慮した設計となっている。
検証用データセット:
本手法の検証には、ヒト膵島細胞データセット、10x Genomics プラットフォームで得られたヒト骨髄単核細胞データ、および C57BL/6J マウスの脳皮質から得られた scRNA-seq、scATAC-seq、STARmap (spatially-resolved transcript amplicon readout mapping)、osmFISH (cyclic single-molecule fluorescence in situ hybridization) データを使用した。マッピングやアライメントの検証には、Dobin et al. Bioinformatics 2013 や Langmead et al. GenomeBiol 2009 などの標準的ツールで処理されたデータを使用した。また、scRNA-seq データの取得には Smart-Seq2 (Single-cell RNA-seq Smart-Seq2) 技術で得られたプロファイルも参照データとして用いた。