- 著者: Yuhan Hao, Stephanie Hao, Erica Andersen-Nissen, William M. Mauck III, Shiwei Zheng, Andrew Butler, Maddie J. Lee, Aaron J. Wilk, Charlotte Darby, Michael Zager, Paul Hoffman, Marlon Stoeckius, Efthymia Papalexi, Eleni P. Mimitou, Jaison Jain, Avi Srivastava, Tim Stuart, Lamar M. Fleming, Bertrand Yeung, Angela J. Rogers, Juliana M. McElrath, Catherine A. Blish, Raphael Gottardo, Peter Smibert, Rahul Satija
- Corresponding author: Rahul Satija (New York Genome Center / New York University); Peter Smibert (New York Genome Center)
- 雑誌: Cell
- 発行年: 2021
- Epub日: 2021-05-31
- Article種別: Original Article
- PMID: 34062119
背景
単一細胞RNAシーケンス (scRNA-seq) 技術は、不均一な組織における新規細胞種や細胞状態の発見に大きく貢献してきた。しかし、トランスクリプトームのみでは機能的に異なる免疫細胞集団を明確に分離できない場面が依然として多いことが課題として認識されている。例えば、CD4陽性T細胞とCD8陽性T細胞は、ショートリードscRNA-seqデータではRNA特徴量のみで完全に識別することが困難な場合がある。これは、T細胞のRNA含量が一般的に低く、RNase発現が高いという技術的制約に起因すると考えられる (Andreeff et al. 1978; Lu et al. 2018)。さらに、エフェクターT細胞、制御性T細胞、γδT細胞、MAIT (mucosal associated invariant T) 細胞など、機能的に多様なT細胞サブセットが転写産物プロファイルのみからでは区別しにくい事例も多く報告されている (Ding et al. 2020; Mereu et al. 2020)。Stubbington et al. Science 2017は単一細胞トランスクリプトミクスによる免疫系の包括的解析の可能性を示唆したが、RNA単独の限界は依然として残された課題であった。
一方、CITE-seq (Cellular Indexing of Transcriptomes and Epitopes by sequencing) に代表されるマルチモーダル単一細胞技術（転写産物、表面タンパク質、クロマチンアクセシビリティ、DNAメチル化、空間情報などの同時取得）が成熟するにつれて、複数のデータモダリティを統合して細胞状態を定義する計算手法の必要性が高まっていた。既存の統合手法（MOFA+ (multi-omics factor analysis plus)、totalVIなど）は、全細胞および全データセットに対して均一な統合戦略を仮定する傾向がある。しかし、実際には各モダリティのデータ品質や情報量は細胞ごとに大きく異なるため、既存の手法では柔軟な対応が不足していた。例えば、抗体パネルが標的細胞タイプを包括的に網羅している場合はタンパク質モダリティが優位となるが、重要なマーカーが含まれない場合や未知の細胞種ではRNAモダリティが優位となる。このような細胞特異的なモダリティ依存性を柔軟に扱える統合手法の開発が課題として残されていた。Stuart et al. Cell 2019による包括的統合手法の開発はscRNA-seq解析の標準化をもたらしたが、マルチモーダルデータへの適用には未解明な部分が多く残されていた。
マルチモーダル単一細胞データの統合における主な不足点として、(1) 各モダリティの相対的有用性を細胞単位で学習する手法の欠如、(2) 異なるデータ品質を動的に補正する機構の未確立、(3) 大規模参照アトラスの構築と新規データセットへの迅速なマッピング機能の不足が挙げられる。これらの課題を解決し、マルチモーダル単一細胞データから細胞状態を高解像度で定義する統合的フレームワークの開発が急務であった。
目的
本研究の目的は、複数モダリティの相対的有用性を細胞単位およびモダリティ単位で学習し、統合的に細胞状態を定義する教師なし学習フレームワークである「weighted-nearest neighbor (WNN)」解析を開発することである。さらに、このWNN解析を211,000個のヒト末梢血単核球 (PBMC) のCITE-seqデータに適用し、循環免疫系のマルチモーダル参照アトラスを構築することを目指す。このアトラスを用いて、新規リンパ系亜集団の同定、ワクチン応答やSARS-CoV-2 (severe acute respiratory syndrome coronavirus 2) 感染における免疫応答の解釈に応用する。これらの解析結果は、RツールキットSeurat v4としてオープンソース実装し、広く配布可能にすることも目的とする。
具体的には、以下の目標を設定した。(1) WNNアルゴリズムの理論的基礎と実装を確立する。(2) 228抗体パネルを用いた大規模PBMCアトラスを構築し、高解像度な細胞状態分類を実現する。(3) 新規リンパ系亜集団（例：CD8 T細胞におけるCD103/CD49aの相互排他的発現、NK細胞におけるCD16/CD38の二次元勾配）を同定し、その生物学的意義を検証する。(4) 参照ベースマッピング機能を用いて、COVID-19 (coronavirus disease 2019) 患者およびワクチン接種者の免疫応答を解析する。(5) WNNフレームワークが3モダリティ以上のデータにも拡張可能であることを実証する。これらの目的を達成することで、トランスクリプトームに限定されない、統一されたマルチモーダルな細胞アイデンティティの定義に貢献することを目指す。
結果
WNNによるモダリティ重みの細胞特異的学習とRNA-タンパク質データの相補的統合: 臍帯血CITE-seqデータ (n=8,617細胞、10マーカー抗体パネル) を用いた検証では、RNA k-NNグラフではCD8+ T細胞の近傍にCD4+ T細胞が混在し、CD8/CD4間の誤ったエッジが944本観察された。一方、タンパク質グラフでは誤エッジは12本に留まり、タンパク質データがCD8+ T細胞の識別に優れていることが示された。対照的に、cDC (conventional dendritic cell) はRNAデータで明確に識別できたが、タンパク質データでは他の細胞種と混在した。WNNはこれらのモダリティ間の情報量の差異を細胞特異的なモダリティ重みとして定量化し、CD8+ T細胞には高いタンパク質重みを、形質細胞には高いRNA重みを自動的に割り当てた。この細胞種特異的な重み配分は、既知の表現型知識と高い一致を示した (Figure 1F)。ボーンマローCITE-seqデータ (n=30,672細胞、25抗体パネル) においても、WNN統合によりT細胞亜集団の分離がRNA単独解析 (984誤エッジ) やタンパク質単独解析と比較して大幅に改善され (WNN: 322誤エッジ)、造血前駆細胞の識別も向上した (Figure 2A)。
WNNの堅牢性と他手法との比較: タンパク質データに段階的にガウスノイズを加えるシミュレーション実験では、ノイズレベルの増加に伴いタンパク質モダリティ重みが用量依存的に低下し、十分なノイズが加えられた後には重みが0に収束し、RNAデータへの依存が自動的に切り替わることが確認された (Figure 2C)。WNNをMOFA+およびtotalVIと比較したところ、25種類のタンパク質すべてにおいてWNNが同等以上の性能を示し、特にCD25 (制御性T細胞マーカー) とCD57 (エフェクターT細胞マーカー) の再構成精度で顕著な優位性を示した (Figure 2D)。計算速度においては、WNNが最大15倍高速であった (n=30,672細胞での全データセット解析) (Figure 2E)。ATAC+RNA multiomeデータ (n=11,351 PBMC) では、ATACデータがナイーブCD8+ T細胞とCD4+ T細胞の分離に優れ (ATAC KNN: 373誤エッジ vs RNA KNN: 984誤エッジ)、WNN統合により誤エッジが322本に削減された。MAIT細胞では、ATAC-seq (Assay for Transposase-Accessible Chromatin using sequencing) でアクセシブルなDNA領域がプロ炎症性転写因子RORγt (retinoid-related orphan receptor gamma t) のモチーフに富んでおり、転写産物でのRORγt上昇と一致した (Figure S3)。
211,000 PBMCマルチモーダル参照アトラスの構築と高解像度細胞分類: 228抗体CITE-seqにより取得したn=210,911細胞（健常成人8名、3時点採血）をWNNで解析し、level-2解像度で31の細胞タイプ、level-3解像度で57の細胞サブタイプを定義した (Figure 3D)。具体的には、CD4 T細胞（12クラスター：ナイーブ、TCM (central memory T cell)、TEM (effector memory T cell)、Treg (regulatory T cell) など）、CD8 T細胞（12クラスター：ナイーブ、TCM、TEM、TEM_KLRC2 (killer cell lectin like receptor C2)、MAITなど）、NK細胞（6クラスター：CD56bright、CD56dim、増殖性NKなど）、樹状細胞（cDC1、cDC2、pDC (plasmacytoid dendritic cell)、ASDC (AXL+ SIGLEC6+ dendritic cell)）、および各種B細胞サブセット（ナイーブ、メモリー、形質細胞、形質芽細胞）が識別された。scRNA-seq単独では分離不可能だった細胞集団の明確な識別が実現し、マルチモーダル統合の優位性が示された (Figure 3B, C)。各クラスターはRNAとタンパク質の両方で明確なバイオマーカーを示し、24サンプル全体で一貫性を保持していた (Figure 4A)。
新規リンパ系亜集団の同定と検証: CD8 T細胞において、インテグリンCD103とCD49aの相互排他的発現パターンを同定した (Figure 5A)。CD8+ CD103+ T細胞はインテグリンβ-7を高発現したが、CD8+ CD49a+ T細胞では発現が欠如していた (Figure 5B)。両集団はカノニカルな組織常在性マーカーCD69を発現しておらず (Figure S5)、組織常在化に向けて分化中の細胞と考えられた。両群間で差次的発現遺伝子モジュールを同定し、T細胞活性化、分化、シグナル応答、走化性に関連する遺伝子オントロジー項が濃縮された (Figure 5E, F)。フローサイトメトリーにより、独立した健常人PBMC (n=5 CMV陽性、n=3 CMV陰性) で両集団の存在を検証した (Figure 5C, D)。NK細胞では、CD16発現に相関する遺伝子モジュール（細胞傷害性関連）と、CD38発現に相関する第二の勾配（CMV関連適応性NK細胞シグネチャ）を同定した。CD38とCD16の相関は最小限 (r < 0.1) であり、NK細胞が二次元勾配上に分布することを示した (Figure 5K)。KLRC2 (NKG2C) はCD38勾配と負相関し、FCER1G (シグナルアダプター) は正相関した (Figure 5I, J)。
TCRクローン拡大とCMV関連免疫応答: 10x 5’ VDJキットによりTCR (T cell receptor) 配列を取得した細胞 (n=49,147) を解析し、31個の拡大クローン (10細胞以上) を同定した。これらのクローンはCD8_TEM_4、CD8_TEM_5、CD4 CTL (cytotoxic T lymphocyte) クラスターに限定されていた。各クローン内の細胞は極めて類似した分子プロファイルを示した (Figure 5H)。CMVペプチドプール刺激試験により、5名がCMV陽性、3名が陰性と判定され、拡大クローン細胞の91%がCMV陽性者に由来していた (Table S3)。MAIT細胞と不変性NKT細胞は予想通り半不変性TCRレパートリーを示し、複数の健常人で共通のTCR配列が観察された (Figure 5G)。
ワクチン接種後の自然免疫応答の解析: ワクチン接種前 (day 0) と接種後3日 (day 3) の比較では、CD14+古典的単球とCD16+非古典的単球で強い転写産物応答が観察された (Figure 6A)。両細胞型でタイプIインターフェロン応答遺伝子 (62遺伝子モジュール) が共通に上昇した (Figure 6B)。Siglec-1 (CD169) はタンパク質レベルでday 3サンプルのみでロバストに誘導された (p < 0.001) (Figure 6C)。樹状細胞サブセット間では、CD1C+ cDC2のみが顕著な応答を示し、CD141+ cDC1、ASDC、pDCでは差次的発現遺伝子が検出されなかった。リンパ球では全体的な比率変化は観察されなかったが、NK_proliferationクラスター（細胞周期遺伝子発現）がday 3とday 7で一貫して増加した。古典的単球と非古典的単球の比率はday 0からday 3で有意に変化し (p < 0.05, paired Wilcoxon test)、day 7でベースラインに復帰した (Figure 6F)。
COVID-19患者への参照ベースマッピングとMAIT細胞の枯渇: Wilk et al. (2020) のCOVID-19患者PBMCデータ (n=44,721細胞、患者7名 + 健常対照6名) を、構築した参照アトラスへsPCAベースマッピングにより射影した。教師なしクラスタリングではT細胞が6クラスターに分類されたが、教師ありマッピングにより12のlevel-2細胞タイプに細分化された (Figure 7C, D)。形質芽細胞頻度の増加（重症COVID-19）は教師なし解析でも検出されたが、MAIT細胞の枯渇は教師ありマッピングでのみ検出された。scRNA-seqベースのMAIT細胞頻度予測値とCyTOF (cytometry by time of flight) 測定値の相関は高かった (Pearson r = 0.911, n=13サンプル) (Figure 7E)。独立した検証コホート (n=16追加患者) のCyTOF解析でも、COVID-19患者におけるMAIT細胞の有意な枯渇が確認された (p < 0.05, unpaired Wilcoxon test) (Figure 7F)。このMAIT細胞枯渇は、抗原特異的なバリア組織への遊走を反映する可能性が示唆された。
考察/結論
先行研究との違い: 本研究は、Stuart et al. Cell 2019による包括的統合手法と異なり、WNN (weighted-nearest neighbor) は各モダリティのデータ品質と情報量を細胞単位で動的に学習する点が特徴である。これにより、抗体パネルの包括性が高い場面ではCITE-seqタンパク質データが主導し、未知細胞種が含まれる場面ではRNA転写産物が補完するという適応的統合が実現される。WNNは、ノイズの多いモダリティの影響を自動的に軽減することで、データ品質の変動に対する堅牢性も示した。
新規性: 本研究で初めて、228抗体パネルを用いた211,000個のPBMC (peripheral blood mononuclear cell) の大規模マルチモーダル参照アトラスを構築し、従来のscRNA-seq (single-cell RNA sequencing) では識別困難なリンパ系亜集団を高解像度で同定した。特に、CD8 T細胞におけるCD103/CD49aの相互排他的発現、NK細胞のCD16/CD38二次元勾配、CMV (cytomegalovirus) 関連クローン拡大の詳細な特性化は、これまで報告されていない新規知見である。これらの発見は、従来のフローサイトメトリーやCyTOF (cytometry by time of flight) では捉えられなかった細胞状態の連続体を明らかにした。
臨床応用: 本知見は、構築されたPBMC参照アトラスへ患者検体を射影することで、細胞タイプ比率と細胞状態シフトを標準化された解像度で迅速に評価できるため、臨床応用に直結する。ワクチン試験、免疫療法バイオマーカー研究、自己免疫疾患研究、感染症応答解析において広い有用性が期待される。肺がん・胸部腫瘍学の文脈では、腫瘍微小環境の免疫細胞サブセット解析への応用が想定され、WNNによる免疫機能サブセットの精密識別が腫瘍浸潤免疫細胞のプロファイリング精度向上に寄与する。特に、PD-1/PD-L1 (programmed cell death protein 1/programmed death-ligand 1) チェックポイント阻害薬の応答予測や、腫瘍関連マクロファージ・樹状細胞の機能的亜集団の同定に有用である。
残された課題: 今後の検討課題として、(1) 組織特異的アトラス（肺、腫瘍微小環境、脳など）の構築、(2) 参照に存在しない新規細胞状態が含まれるクエリデータでの対応、(3) 空間トランスクリプトミクスとの統合、(4) 多施設データでのバッチ効果が大きい場合のマッピング精度の検証が挙げられる。WNNフレームワークは3モダリティ以上への拡張が可能であり、今後の技術進展（クロマチン状態、DNAメチル化、空間位置の同時計測）に対応する基盤を提供する。Seurat v4は本論文公開後、scRNA-seq解析のデファクトスタンダードとして急速に採用が広がっており、学術・臨床研究の両領域で広く利用されている。
方法
WNNアルゴリズム: WNNアルゴリズムは、各細胞についてRNAとタンパク質それぞれのモダリティで独立にk=20の最近傍細胞を特定することから開始する。次に、各モダリティのクロスモダリティ予測精度を計算する。これは、自身のモダリティの近傍細胞の平均発現量から相手モダリティの発現量を予測し、実際の測定値と比較することで行われる。この予測精度に基づいて、各細胞に固有のモダリティ重みを計算する。モダリティ重みは非負であり、合計が1となるようにsoftmax変換される。最終的に、これらの重み付き距離（RNAとタンパク質の加重平均）を用いてマルチモーダル最近傍グラフを構築する。このWNNグラフは、下流のクラスタリング、UMAP (uniform manifold approximation and projection) 可視化、および細胞軌跡解析に利用される。アルゴリズムの詳細な数学的記述はSTAR Methodsに記載されている。
CITE-seqアトラス構築: 健常成人8名から末梢血単核球 (PBMC) サンプルを採取した。これらのサンプルは、HIVワクチン試験 (HVTN 087, NCT01578889) の参加者から、ワクチン接種前 (day 0)、接種後3日 (day 3)、7日 (day 7) の3時点で採血された。BioLegend TotalSeq-A抗体パネル（最大228抗体）を用いたCITE-seqにより、合計211,000個のPBMCを取得した。具体的には、10x Genomics 3’ v3キットを用いて161,764細胞（228抗体）を、10x Genomics 5’ VDJキットを用いて49,147細胞（54抗体 + TCR/BCR）を処理した。品質管理とダブルレット除去後、WNNクラスタリングを適用し、level-2解像度で31のサブセット、level-3解像度で57のサブセットの細胞状態を定義した。
追加データセットと参照ベースマッピング: WNN解析の汎用性を評価するため、SARS-CoV-2感染患者PBMCデータ (Wilk et al. 2020, n=7患者 + 6健常対照)、10x Genomics Multiome ATAC+RNAデータ (n=11,351 PBMC)、ASAP-seq (Assay for Single-cell ATAC and Protein by sequencing) データ（RNA + ATAC + タンパク質の3モダリティ同時計測、骨髄細胞）、およびSHARE-seq (Simultaneous High-resolution Assay for RNA Expression and Chromatin Accessibility) データ（マウス皮膚細胞, n=34,774）にWNNを適用した。参照ベースマッピング機能では、sPCA (supervised principal component analysis) を用いて、参照データセットのWNNグラフ構造を最大限保持する転写産物ベクトルを学習した。その後、クエリデータセットをこの参照アトラスに射影し、細胞タイプラベルとタンパク質発現量を転送した。
統計解析: 差次的発現解析にはWilcoxon rank-sum testを用い、p値が0.05未満を有意とした。モダリティ重みの比較にはpaired Wilcoxon testを使用した。CMV (cytomegalovirus) 感染状態の判定は、CMVペプチドプール刺激後の細胞内サイトカイン染色 (ICS) により、CD8+ T細胞のIFN-γ (interferon-gamma) 産生に基づいて行った。
実装: 開発されたWNNフレームワークは、Seurat v4としてGitHubおよびCRANでオープンソース公開された。また、Azimuthウェブアプリケーション (https://azimuth.hubmapconsortium.org/) を通じて、ユーザーは自身のデータセット（最大50,000細胞）を5分以内に参照アトラスにマッピングできる機能が提供された。