- 著者: Yoshiaki Yasumizu†, Masaki Hagiwara†, Yuto Umezu, Hiroaki Fuji, Keiko Iwaisako, Masataka Asagiri, Shinji Uemoto, Yamami Nakamura, Sophia Thul, Azumi Ueyama, Kazunori Yokoi, Atsushi Tanemura, Yohei Nose, Takuro Saito, Hisashi Wada, Mamoru Kakuda, Masaharu Kohara, Satoshi Nojima, Eiichi Morii, Yuichiro Doki, Shimon Sakaguchi, Naganari Ohkura (†joint first authors)
- Corresponding author: Naganari Ohkura (Department of Experimental Immunology, Immunology Frontier Research Center, Osaka University)
- 雑誌: NAR Cancer
- 発行年: 2024
- Epub日: 2024-04-15
- Article種別: Original Article
- PMID: 38751935
背景
DNAメチル化はゲノム全体にわたる安定した主要エピジェネティック修飾であり、遺伝子発現が細胞状態や環境条件によって変動しやすいのと対照的に、細胞の系統と特性を高再現性かつ安定的に反映する。近年、さまざまな細胞種にわたるDNAメチル化の網羅的プロファイリングにより、その細胞特異性と生物学的意義が解明されてきた。バルク組織やcfDNA (cell-free DNA) のメチル化プロファイルから細胞種構成比を数理的に推定するデコンボリューション解析は、腫瘍浸潤免疫細胞のプロファイリング・がん動態モニタリング・COVID-19や心筋梗塞などの病態解明において有望な臨床応用として期待を集めている。特に、免疫チェックポイント阻害薬 (ICI: immune checkpoint inhibitor) の治療効果予測において、PD-L1 単分子発現よりも腫瘍浸潤免疫細胞の組成プロファイリングが優れることが示されており (Havel et al. NatRevCancer 2019)、DNAメチル化ベースのデコンボリューションがその臨床実装を支える基盤技術として注目されてきた。さらに Chakravarthy らによるパンがん横断的なDNAメチル化デコンボリューション (Chakravarthy et al. NatCommun 2018) は腫瘍免疫環境の定量的評価を可能にしたが、活性化・ナイーブ・疲弊状態といった免疫細胞のより詳細なサブセット状態を識別する解像度は達成されていなかった。こうした精緻なサブセット分解が不可能であることがDNAメチル化ベースのTMEプロファイリングの技術的限界として残されており、解決が急務であった。
既存のDNAメチル化デコンボリューション手法には複数の重大な方法論的制約が存在する。CIBERSORT (Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts; Newman et al. NatMethods 2015)・ARIC (Accurate methylation-based cell Ratio estimation In Cancer)・EpiDISH (Epigenetic Dissection of Intra-Sample Heterogeneity) [RPC (Robust Partial Correlations)・CP (Constrained Projection)・QP (Quadratic Programming)] 等は、細胞種ごとのシグネチャーマトリクスを定義したうえで線形モデルまたはSVR (support vector regression) を適用する方式をとっており、高分散の大規模データセットにおける複数細胞種の同時推定への適応性が不足していた。また多くの既存手法はPBMC (peripheral blood mononuclear cell) デコンボリューション評価に重点が置かれており、腫瘍組織やcfDNAを対象とした評価には多様な組織・免疫細胞状態を包含した大規模リファレンスパネルの整備が手薄であった。さらに、WGBS (whole genome bisulfite sequencing)・RRBS (reduced representation bisulfite sequencing)・Illumina メチル化アレイ・ナノポアシーケンサーといった異なるプラットフォームを統一的にカバーするデコンボリューションモデルは存在せず、臨床実装における柔軟性に欠けるという問題があった。こうした課題を解消するため、本研究ではニューラルネットワークアーキテクチャとmixupデータ拡張を組み合わせた新規デコンボリューションツール MEnet (Methylation-based Estimation of cellular networks) の開発が着手された。既存手法に不足していたのは、異なる計測プラットフォームを統一的にカバーしながら詳細な免疫細胞サブセット状態まで識別できる、大規模コホートへの実用的な高速デコンボリューション手法であった。すなわち、WGBS・アレイ・ナノポアシーケンサーを横断して適用可能であり、かつ活性化・疲弊・ナイーブといった免疫細胞のサブセット状態を高解像度で識別できる統合的計算基盤は存在せず、この空白が解決されていなかったことが、この分野の中心的な未達課題として残されていた。
目的
多様な組織・腫瘍浸潤免疫細胞のDNAメチル化リファレンスデータ (945サンプル) でニューラルネットワークモデルを訓練し、血液・固形腫瘍・cfDNA・FFPE (formalin-fixed paraffin-embedded) 組織を含む幅広いサンプルタイプに汎用的に適用可能な高精度・高速な細胞デコンボリューションツール MEnet を開発・検証する。さらに難治性がんである肝内胆管がん (ICC: intrahepatic cholangiocarcinoma) 72例の腫瘍微小環境 (TME: tumor microenvironment) 解析に応用し、メチル化デコンボリューション由来の細胞組成プロファイルと患者の全生存・再発との関連を明らかにすることで、MEnet の臨床的有用性を実証する。
結果
モデル最適化とベンチマーク比較：
151サンプル (20カテゴリ、各カテゴリ最低5サンプル) を用いたモデル比較では、NN (neural network) + mixupモデルが精製細胞・混合細胞の双方で最高精度を達成した。精製細胞データでのMSEはNN + mixupが最低値を記録し、mixupなしNN・ランダムフォレスト・SVR・LightGBM (Light Gradient Boosting Machine)・NNLS (non-negative least squares) をすべて下回った (Fig. 1B)。特に混合細胞サンプルでは、NN + mixupのMSEが0.00730と他手法の0.018-0.084の誤差範囲を大幅に下回り、mixupデータ拡張によるoverfitting抑制の効果が現実の混合サンプルで特に顕著であることが示された。この精度優位性は6-fold交差検証の独立した各フォールドにおいても一貫して再現され (n=6 folds)、MEnetの汎化性能の堅牢さが確認された。一方、mixupなしNNは精製細胞での精度がわずかに良いものの混合細胞での精度が劣り、現実的な臨床サンプルには NN + mixupが優位であることが明確になった。既存デコンボリューション専用ツール (ARIC, CIBERSORT, RPC, CP/QP) と同一リファレンスで比較した場合も、MEnetは精製・混合細胞いずれの条件で上回る精度を示し (Fig. 1B)、B細胞・CD4+T・CD8+T・NK細胞を等割合で混合したサンプルでの細胞比率予測でも MEnetが最も小さな誤差を示した (Fig. 1C)。計算速度においてはMEnetの予測時間が既存手法比で781-fold以上高速であり (Supplementary Fig. S1C)、大規模臨床コホートへの実用的な適用に耐える水準が確認された。
MEnetのアーキテクチャと特徴領域同定：
945サンプルを包含するリファレンスデータから、ゲノム全体の1000 bpビンタイリングと2段階カイ二乗検定 (FDR上位2000・メチル化率差>40%) により合計87,726の特徴領域を同定した (Fig. 2C, D)。これらの特徴領域はプロモーターコア領域に最も集積しており (log10 Obs/Exp = 1.33、P < 0.001)、次いで5’ UTR (log10 Obs/Exp = 1.00)、エクソン (log10 Obs/Exp = 0.62) の順で有意な富化が確認された (Fig. 2E)。Bayesian最適化 (>3000 trial; EarlyStopping patience=1000) によって最終的な最適アーキテクチャ (隠れ層2層・ノード数1030) が確定した。シーケンス深度評価では、各条件につき n=5 independent replicates (独立サブサンプリングシード) を用いた検証において、RRBSは約2×10^8 bp相当で精度が飽和し、WGBSはその約10倍 (約2×10^9 bp) で飽和するが、WGBSでも3x以下のゲノムカバレッジで十分な精度を達成可能であることが示され (Fig. 3)、臨床適用における現実的なシーケンス深度要件が明確化された。
多様なサンプルタイプへの性能検証：
TCGA Illumina 450kデータ (17がん種) への適用では、MEnet由来腫瘍純度スコアがESTIMATE・ABSOLUTE (CNVベース)・IHC・CPEと複数がん種にわたって有意な正相関を示し (Fig. 4A)、MEnet由来免疫細胞スコアもLUMPと有意に相関した (Fig. 4B)。実臨床サンプル検証として、PBMC精製細胞DNAの人工混合12タイプ (GEO (Gene Expression Omnibus) accession GSE167998) を用いた独立した評価では、MEnet予測値が活性化・ナイーブ・疲弊を含む詳細な細胞状態分類においても実測値と高い一致を示した (Fig. 4D)。腫瘍組織 (卵巣・胃・皮膚・大腸がん) のフローサイトメトリーとの独立した比較 (n=11 腫瘍サンプル) では、リンパ球系統で強い相関が認められ、活性化CD8+T細胞など詳細なサブセットでも再現性が確認された (Fig. 4E)。骨髄系統では組織解離過程の細胞損傷による相関の弱化が観察され、今後の改善課題として認識された。FFPE組織への適用可能性の評価では、FFPE由来と新鮮組織由来のMEnet推定細胞比率間にPearson r = 0.75 (p < 0.001) の強い相関が確認された (皮膚がん n=6、Fig. 4F)。卵巣がん患者血液からのcfDNA (n=13) への適用では、MEnet推定卵巣由来比率とichorCNA由来CNA腫瘍分率の間にPearson r = 0.92の高い相関が認められ (Fig. 5)、血液単一検体でのリアルタイムがんモニタリングへの応用可能性が示唆された。
ICC腫瘍微小環境と患者予後の関連：
京都大学病院でICCに対して外科手術を施行した n=72例 (2004-2016年) の腫瘍WGBSデータにMEnetを適用した。同一患者の隣接正常肝とICC腫瘍部位のペア比較 (n=3 matched pairs) では、正常肝は肝細胞成分が主体だったのに対し、腫瘍部位では肝臓比率の低下と胆管・膵・腸管などの他消化器内胚葉組織の割合増加が観察され、ICCが胆管起源であることと整合的なパターンを示した (Fig. 6B)。腫瘍部位にのみ豊富な免疫細胞浸潤が観察され、MEnetによるTMEプロファイリングの有効性が確認された。CNA由来腫瘍分率とMEnet由来非腫瘍細胞比率の間には有意な負の相関が観察された (Spearman r = -0.61, p < 0.001、Fig. 6C)。72例×39細胞種の組成プロファイルをNMF (10成分) で次元削減し、Cox比例ハザード回帰で全生存との関連を検討した結果、NMF4 (高荷重: 活性化Tconv (conventional T cell; 従来型T細胞), 疲弊CTL (cytotoxic T lymphocyte), メモリーTreg) およびNMF7が高値の症例で良好な生存と関連し、NMF6 (高荷重: 心血管系細胞) 高値の症例で不良な生存と関連した (Fig. 6E)。カプランマイヤー解析でも、NMF4高値群は生存延長・NMF6高値群は生存短縮という有意な生存差が確認された (Fig. 6F)。ブートストラップ法 (100回反復) によるC-index評価では、MEnet由来NMFスコアを用いた全生存予測が従来のがんステージ分類を有意に上回り (p < 0.001、Fig. 6G)、再発予測においてもNMF4高値が低再発と関連したものの再発予測精度ではがんステージ分類が優位だった (p = 0.00401)。
考察/結論
本研究はDNAメチル化情報を活用したニューラルネットワークベースの細胞デコンボリューションツール MEnet を確立し、多様なサンプルタイプへの適用可能性と腫瘍微小環境解析における臨床的有効性を実証した。
既存手法との違いと優位性：これまでの研究で広く用いられてきた CIBERSORT・ARIC・EpiDISH 等の線形モデルや限定的な SVRベース手法と異なり、MEnetはニューラルネットワークの非線形性とmixupデータ拡張を組み合わせることで、大規模・高分散データセットにおける複数細胞種の同時推定精度を大幅に向上させた。混合細胞サンプルでの MSE 比較において既存手法が0.018-0.084の誤差を示すのに対して MEnetは0.00730を達成し、速度面でも 781-fold 以上の高速化を実現した点は既報と異なる特徴である。また、WGBS・RRBS・Illumina 450k/EPICアレイ・ナノポアシーケンサーを含む全ての主要プラットフォームを統一モデルでカバーする初の事例であり、プラットフォーム固有の従来モデルと対照的に臨床実装の柔軟性が格段に向上している。
本研究で新規に達成された知見：MEnetが新規に実現したのは、腫瘍浸潤免疫細胞の詳細なサブセット状態 (活性化・疲弊・ナイーブ) まで解像したDNAメチル化ベースのデコンボリューションである。疲弊CD8+T・活性化CD8+T・ナイーブTreg・エフェクターTregの4サブセットを独自WGBSデータとしてリファレンスに組み込んだ点は新規の試みであり、これまで報告されていない細胞状態識別の解像度を実現した。ICC の n=72例という大規模コホートにおいて、メチル化デコンボリューション由来TMEスコアが病理学的がんステージを超える全生存予測精度をもたらすことを本研究で初めて示したことは重要な知見である。NMF4が反映する免疫活性型TMEパターン (活性化T細胞と Treg の共存・高免疫原性) と良好な予後の関連は ICC の既知「炎症型」良好予後サブクラスと生物学的に整合し、手法の妥当性を支持する。
臨床的意義と橋渡し研究への展望：MEnetはフローサイトメトリーやmIHC (multiplex immunohistochemistry; 多重免疫組織化学) が実施困難な日常臨床サンプル (FFPE・少量cfDNA) に適用可能であり、bench-to-bedsideの実装が期待できる。cfDNAへの応用は、血液単一検体で腫瘍微小環境と免疫細胞動態をリアルタイムにモニタリングする液体生検の次世代的アプローチを提供する。ICI治療に伴う末梢血免疫細胞動態の観察を通じたirAE (immune-related adverse events; 免疫関連有害事象) の早期検知への臨床的有用性も今後期待される。モデルおよび訓練済み重みはGitHub・Zenodo・Figshareで公開されており、外部コホートへの即時適用が可能である。
残された課題と今後の展望：MEnetのナイーブTconv・メモリーTconv等一部サブセットの絶対値推定は実測値との乖離が認められ、臨床への直接応用には慎重な解釈が必要なlimitationが存在する。骨髄系細胞は組織解離過程での細胞損傷の影響を受けやすく、腫瘍組織フローサイトメトリーとの相関が弱まるという残された課題がある。ICC72例という比較的小規模なコホートでの予後解析は探索的であり、独立した大規模多施設コホートでの検証が今後の検討として不可欠である。シングルセルWGBSやscRNA-seq (single-cell RNA sequencing) 由来メチル化推定データとの統合によって希少細胞亜集団の解像度をさらに高める余地があり、多がん種前向きコホートでのICI応答予測・治療効果モニタリングへのMEnet適用がfuture researchとして期待される。
方法
リファレンスデータセット構築: ENCODE (Encyclopedia of DNA Elements)・Roadmap Epigenomics・Blueprint・ヒトメチロームアトラスの公開データベースに加え、大阪大学・京都大学の外科的切除腫瘍組織から独自に単離した腫瘍浸潤免疫細胞のWGBSデータを統合し、計945サンプルのゲノムワイドDNAメチル化データセットを構築した。細胞種分類は29 Majorグループ・39 Minorグループ。データ取得手法はWGBS・RRBS・Illumina 450k/EPIC (Infinium MethylationEPIC) アレイ・ナノポアシーケンサー (MinION, Oxford Nanopore Technologies; Rapid Sequencing Kit SQK-RAD004; R9.4.1フローセル) を包含。独自に取得した腫瘍浸潤リンパ球はFACSAria III (BD Biosciences) でソート: 疲弊CD8+T (CD3+CD8+CD45RA-PD1+)・エフェクターCD8+T (CD3+CD8+CD45RA-PD1-)・エフェクターTreg (CD3+CD4+CD25highCCR8+)・Treg (CD3+CD4+CD25high)。シーケンスデータはJGA (Japanese Genotype-phenotype Archive)/NBDC (National Bioscience Database Center) 制限公開 (Accession: JGAS000676)、人工混合検証データはGSE167998。
特徴領域選定とモデル設計: ゲノム全体を1000 bpビンにタイリングし、Majorグループのカイ二乗検定 (FDR (false discovery rate) 上位2000・メチル化率差>40%) とMinorグループ内カイ二乗検定の2段階選定により各細胞種の特徴領域を最大2000箇所ずつ同定、合計87,726領域を入力特徴とした。モデルはPyTorch (v1.6.0) 実装のMLP (multilayer perceptron; 多層パーセプトロン)。6-fold 交差検証とBayesian最適化 (Optuna CmaEsSampler; >3000 trial; 探索空間: 層数2-10・ノード数30-3000・ドロップアウト率0-0.3・活性化関数relu/tanh/leakyrelu・optimizer Adam/RMSprop/SGD (stochastic gradient descent)・学習率1e-5〜1e-1; EarlyStopping patience=1000) でハイパーパラメータを最適化し、最終アーキテクチャとして隠れ層2層・ノード数1030を確定。データ拡張にmixup (最大15サンプルの加重ランダム混合) を採用。計算は分散学習 (MySQL) + GPU (Intel Xeon Platinum 8180 @ 2.5GHz / NVIDIA Tesla V100 PCIe 32GB) で実施。
評価・統計解析: モデル比較はMSE (mean squared error) をscikit-learnで算出し、ARIC (Accurate methylation-based cell Ratio estimation In Cancer; v0.1)・CIBERSORT (Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts)・EpiDISH (v2.12.0) のRPC/CP/QP・NNLS (non-negative least squares) と同一リファレンスで比較。外部検証にはTCGA (The Cancer Genome Atlas) Illumina 450kデータ (17がん種) を使用し、ESTIMATE (Estimation of STromal and Immune cells in MAlignant Tumours)・ABSOLUTE (absolute quantification of somatic DNA alterations)・IHC (immunohistochemistry; 免疫組織化学)・CPE (cancer purity estimation) との腫瘍純度Pearson相関係数 (一側相関検定)、LUMP (Leukocytes Unmethylated in CpG loci Present) との免疫スコアPearson相関を評価。実臨床サンプル検証: フローサイトメトリー比較 (n=11 腫瘍サンプル)・FFPE vs 新鮮組織Pearson相関 (皮膚がん n=6)・cfDNA vs CNA (copy number alteration) Pearson相関 (卵巣がん n=13; CNA推定はichorCNA、hidden Markov model)。ICC72例 (京都大学病院2004-2016年手術例; 特性はSupplementary Table S3) のWGBSデータにMEnet適用後、NMF (non-negative matrix factorization; scikit-learn v0.24.1; 10成分; min-maxスケーリング) で次元削減。生存解析はCox比例ハザードモデル (lifelines v0.27.0; 共変量: NMF10成分+がんステージ整数値)、Kaplan-Meier曲線の群間差はlog-rank検定、NMF-臨床情報相関はSpearman相関 (R package corrplot) で可視化。全生存予測精度評価はブートストラップ法 (100回反復; 62例訓練/10例テスト) によるC-index (concordance index) 算出とMann-Whitney U検定でがんステージ分類と比較した。