- 著者: Zuguang Gu, Roland Eils, Matthias Schlesner
- Corresponding author: Matthias Schlesner (German Cancer Research Center (DKFZ), Heidelberg, Germany)
- 雑誌: Bioinformatics
- 発行年: 2016
- Epub日: 2016-05-20
- Article種別: Applications Note
- PMID: 27207943
背景
ゲノムデータに隠されたパターンを明らかにする上でヒートマップは基盤的な可視化手法であり、特に遺伝子発現解析 (Eisen et al. 1998) やメチル化プロファイリング (Sturm et al. 2012) において広く利用されてきた。ゲノムワイドな発現クラスタリングをヒートマップとして示した Eisen らの手法は現在も遺伝子発現解析の標準的視覚化として踏襲されており、メチル化解析においても同様の可視化パターンが確立されている。Gao et al. SciSignal 2013 が開発した cBioPortal の OncoPrint 機能はゲノム変異パターンの可視化に優れた成果をあげているが、こうした個別ツールはそれぞれ特定のデータ型に特化しており、複数の異なる情報源を統合的に表示する汎用フレームワークとしては機能しない。
次世代シーケンシング (next-generation sequencing) 技術の急速な普及により、同一患者コホートに対して変異・コピー数変化・転写産物量・エピゲノム修飾など多種類のゲノムデータが同時に生成される状況が標準化している。このような多次元データを統合解析する際、情報層間の相関関係を一つの視覚的フレームで把握できる可視化ツールの必要性はきわめて高い。しかし R プログラミング環境で利用可能な従来のヒートマップツール（基本 heatmap 関数、pheatmap、heatmapplus など）はいずれもアノテーショングラフィックの表示機能が限定的であり、複数の並列ヒートマップをプロットする機能を持たないという不足があった。Cancer et al. Nature 2014 のような大規模コホートゲノムプロファイリング研究が生み出す多種多様なデータ層を単一の可視化フレームに統合し、生物学的文脈情報と合わせてアノテートする柔軟な汎用ツールが手薄であることが研究上の gap in knowledge として認識されていた。高次元ゲノムデータのパターンと相関関係を効率的に明らかにするために、より機能豊富で拡張性の高いヒートマップフレームワークの開発が求められていた。
目的
本研究の目的は、高次元ゲノムデータにおけるパターンと情報源間の関係性を効率的に可視化するための R/Bioconductor パッケージ ComplexHeatmap を開発し、従来ツールが持つ制限を克服することである。具体的には以下の機能を実現することを目指した。(1) + 演算子を用いて異なる情報源のデータ行列を共通の行次元で連結し、複数の並列ヒートマップとして一括表示する機能。(2) 点グラフ・棒グラフ・ボックスプロットなどの定型アノテーションおよびユーザーが独自に設計したグラフィックをヒートマップの行・列に柔軟に付加できるアノテーション管理機能。(3) デンドログラム描画やユーザー定義距離関数を含む柔軟なクラスタリングサポート、k-means 分割やデータフレーム分類によるヒートマップ行スライス化機能、インタラクティブデバイス上での行列サブセット選択機能。加えて、4 つの実際のゲノムデータセットを用いてパッケージの有用性を実証することも目的に含まれる。
結果
ComplexHeatmap パッケージは n=4 の実データセットを通じ、多次元ゲノムデータの複雑なパターンと情報層間の相関関係を効果的に可視化することを実証した。
TCGA 肺腺癌 OncoPrint による多層ゲノム変化の統合表示:
TCGA 肺腺癌コホートにおける n=134 患者、n=38 選択遺伝子のゲノム変化（点変異・欠失・増幅などの体細胞変化）を ComplexHeatmap の拡張 OncoPrint として可視化した (Fig 1A)。従来の OncoPrint（cBioPortal 実装）では単一のゲノム変化マトリクスの表示に留まるが、本実装では ComplexHeatmap の並列ヒートマップ機能を活用して複数の情報層を同一フレームで提示した。n=38 遺伝子は患者間での増幅率に基づいて高増幅群と低増幅・変異主体群の 2 グループに行スライス化して分割し、視覚的に区別できる構成とした。ヒートマップの行には各遺伝子における全患者の変化総数を示す棒グラフアノテーションを、列には各患者における全遺伝子の変化総数を示す棒グラフアノテーションをそれぞれ追加し、バリアント (variant) 分布の定量的把握を可能にした。さらに OncoPrint の右側に第 2 の独立したヒートマップを並列配置し、各変異遺伝子の生物学的機能（Gene Ontology 機能カテゴリー）を色分けで表示した。2 枚のヒートマップが行を共有して並置されることで、ゲノム変化パターンと機能的カテゴリーの対応関係を一目で読み取ることができる。解析の結果、高度に増幅された遺伝子群は細胞応答プロセス (cellular response processes) に変異主体群の 2-fold 以上の頻度で濃縮される傾向が示された一方、変異を主体とする遺伝子群は発生 (development) および代謝 (metabolism) に関連するプロセスに濃縮されることが明らかになった。コピー数変化と点変異が異なる生物学的文脈で機能的役割を担う可能性を示唆するこの知見は、生物学的機能ヒートマップを OncoPrint と並列配置することで初めて可視化されたものであり、単一ヒートマップしか提供しない従来ツールでは 1 枚の図で完結させることができなかった。
マウス T 細胞の単細胞 RNA-seq による細胞不均一性と共発現構造の同時可視化:
マウス T 細胞の単細胞 RNA-seq (single-cell RNA sequencing) データを用いて、ComplexHeatmap が高次元発現データの細胞間不均一性と遺伝子共発現構造をどのように同時可視化できるかを実証した (Fig 1B)。n=721 個の選択遺伝子の発現プロファイルを階層的クラスタリングで解析した結果、T 細胞は明確な n=2 つのサブポピュレーションに分類された。列デンドログラム上で薄赤色にハイライトされた左側サブポピュレーションは、特定の細胞周期遺伝子サブセットの相対的な発現が高い一方でその絶対発現レベルが低いという特徴を持つ。対して右側サブポピュレーションはこの細胞周期遺伝子サブセットの相対発現が低い反面、その他の選択遺伝子（別の細胞周期遺伝子を含む）が高い絶対レベルで発現していた。特に注目すべき点として、リボ核タンパク質遺伝子群 (ribonucleoprotein complex genes) が両サブポピュレーションを通じて極めて強い共発現 (co-expression) パターンを示し、遺伝子対間の平均ペアワイズ Pearson 相関係数 r=0.83 を超える高相関クラスターを形成した (Fig 1B)。こうした転写コモジュールの存在は、サブポピュレーション分類を超えて多くの T 細胞に共通する転写調節機構の存在を示唆している。さらに、ヒートマップの右側には第 2 の並列ヒートマップとして n=721 遺伝子間のペアワイズ相関係数 (pairwise Pearson correlation) のマトリクスを配置し、相関値の階層的クラスタリングにより得られた行順序を全ての並列ヒートマップで共有した。この設計により、発現プロファイルの類似構造と遺伝子間相関構造が同一行次元上で直接比較可能となり、共発現モジュールの視覚的同定が容易になった。単細胞 RNA-seq データにおける細胞サブポピュレーション識別と遺伝子共発現ネットワーク解析を同一フレーム内で統合した本例は、ComplexHeatmap の並列ヒートマップ機能が単細胞トランスクリプトミクス解析に実用的な価値を持つことを示している。
多層エピゲノムデータの統合可視化:
未発表研究で発見されたパターンを基にランダム生成したデータを用いて、DNA メチル化・遺伝子発現・エンハンサー (enhancer)・遺伝子関連情報の 4 層間の相関を ComplexHeatmap で統合的に可視化した (Supplementary Figure S3)。ヒートマップの各行は 1 つの DMR (differentially methylated region; 差分的メチル化領域) または DMR に関連付けられたゲノムオブジェクトを表す。関連オブジェクトとしては、例えば関連 DMR のメチル化レベルと負相関する発現を持つ最近傍遺伝子、あるいは DMR と空間的に重複するエンハンサーが含まれる。複合ヒートマップの可視化により、高度にメチル化された DMR が遺伝子間領域 (intergenic region) および遺伝子内領域 (intragenic region) に濃縮され、エンハンサーとの重複が稀であることが明確に示された。対照的に、低メチル化 DMR は転写開始点 (TSS; transcription start site) およびエンハンサーに濃縮される傾向を示し、低メチル化が遺伝子活性調節領域に優先的に位置することを視覚的に捉えた。この可視化パターンは、DMR のメチル化状態とゲノム上の機能的文脈（遺伝子構造・調節エレメントとの関係）の間に系統的な関連性が存在することを直観的に提示するものであり、多次元並列ヒートマップが複数のエピゲノム情報層の相互関係を効率的かつ包括的に表現できることを実証した。
グリオブラストーマメチル化プロファイリングの再実装と新規アノテーション追加:
Sturm et al. (2012) によるグリオブラストーマのメチル化プロファイリング論文の Figure 1 を ComplexHeatmap で再実装し、パッケージが既存の複雑な可視化を忠実に再現しながら新たな情報層を容易に追加できることを示した (Supplementary Figure S4)。元の図に対して 2 つの新規ヒートマップを追加した。第 1 の追加ヒートマップは各 CpG (cytosine-phosphate-guanine) サイトと最も近傍の TSS との距離を可視化し、第 2 のヒートマップは CGI (CpG island; CpG アイランド) に対するゲノムアノテーション（CGI 本体・Shore・Shelf・オープンシー）を示した。ヒートマップは CGI アノテーションに従って行スライス化を行い、CGI Shelf およびオープンシーに位置する CpG サイトが高いメチル化レベルと TSS からの高い距離を持つことが示された。一方、CGI 本体に位置する CpG サイトはメチル化レベルが低く TSS 近傍に存在する傾向が確認された。この再実装例は、ComplexHeatmap の行スライス化機能と複数の並列アノテーションヒートマップ機能を組み合わせることで、既存の複雑な多層可視化をわずかなコードで拡張できる柔軟性を具体的に示した。
考察/結論
先行研究との違いと本パッケージの位置づけ:
従来の R ヒートマップツール（heatmap、pheatmap、heatmapplus）と異なり、ComplexHeatmap は複数の並列ヒートマップの連結と、ユーザー定義グラフィックを含む多様なアノテーションの柔軟な統合を単一フレームワークで実現している。既報のツールは単一の行列を視覚化することに特化しており、異なるオミクスデータ層を同一の行次元で整列させた並列ヒートマップとして一括表示する機能が欠如していた。これに対し、本パッケージの + 演算子によるヒートマップ連結は、行の対応関係を保ちながら任意個の並列ヒートマップを結合するという、これまでの研究では提供されていなかった直観的な API (Application Programming Interface) を提供している。また Gao et al. SciSignal 2013 の cBioPortal OncoPrint と比較すると、ComplexHeatmap はゲノム変異可視化に加えて任意の生物学的メタデータを行・列アノテーションとして自由に追加できる汎用性を持ち、特定のデータ型に依存しない可視化プラットフォームとして機能する点で対照的である。pheatmap や heatmapplus と比べた場合も、これらの相違は本パッケージが提供するモジュール式 API 設計に起因しており、ユーザーがヒートマップの構成要素を独立して定義・組み合わせる設計が差異化の核心にある。
新規性と技術的貢献:
ComplexHeatmap は本研究で初めて提供されたオブジェクト指向設計に基づくモジュール式ヒートマップパッケージであり、ヒートマップのカスタマイズ・並列配置・ユーザー定義アノテーション管理を統合した汎用ソリューションは新規のものである。+ 演算子による宣言的な連結記法、HeatmapAnnotation クラスの列/行アノテーション統一管理、ヒートマップグリッドのユーザーカスタマイズ（拡張 OncoPrint 等）、インタラクティブデバイス上での行列サブセット選択など、複数の新規な機能を組み合わせることで、研究者が目的に合わせた複雑な多次元可視化を短いコードで構築できる環境を初めて提供している。TCGA 肺腺癌コホート (n=134) を用いた OncoPrint 解析では、ゲノム変化パターンと生物学的機能の関連を並列ヒートマップで一括表示したことが新規な可視化アプローチとして示された。また、マウス T 細胞単細胞 RNA-seq 解析では n=721 遺伝子の発現構造と遺伝子間相関構造を同一フレームで統合して示した点が新規であり、こうした複合的可視化は従来のツールでは実現が困難であった。
臨床的意義と応用可能性:
ComplexHeatmap はゲノム医学研究のデータ解釈・プレゼンテーション全般において臨床的意義が高い。がんゲノムプロファイリング（変異・コピー数変化・転写産物量・エピゲノム修飾）と臨床情報（治療反応・予後・サブタイプ）を統合的に可視化することで、新規バイオマーカーの発見や疾患サブタイプの分類に直結した知見を効率的に提供できる。Cancer et al. Nature 2014 のような大規模コホート研究のゲノム多次元データを解釈する際に、この種の統合的ヒートマップは臨床現場における個別化医療戦略立案を支援するツールとして bench-to-bedside の橋渡しを担う。マウス T 細胞単細胞 RNA-seq データの可視化実証例は、腫瘍微小環境内の免疫細胞不均一性や細胞サブポピュレーション特異的な遺伝子発現プロファイルを解釈する免疫腫瘍学研究への応用可能性も示唆している。グリオブラストーマメチル化プロファイリングの再実装例は、既存の複雑な図を ComplexHeatmap で再現・拡張できることを実証しており、CpG メチル化と TSS 距離・CGI アノテーションの同時表示という臨床的含意のある多層エピゲノム可視化のテンプレートを提供している。
残された課題と今後の検討:
今後の検討として、非常に大規模な単細胞オミクスデータ（数百万細胞規模）に対するスケーラビリティの最適化が挙げられる。インタラクティブ機能が現状 X11 デバイスに限定されているため、ウェブブラウザベースのインタラクティブ可視化（Shiny アプリ統合など）への対応が将来の展開として期待される。時系列データを扱う動的ヒートマップや、ネットワーク解析結果との統合機能も潜在的な拡張方向として考えられる。Limitation として、R プログラミング環境への依存性と、複雑なユーザー定義アノテーション設計に一定のコーディングスキルが必要な点が挙げられる。また、本論文の実証例のうちメチル化相関データはランダム生成データを用いており、実際の測定データに基づく定量的性能評価が future research として残されている。これらの課題を克服することで、ComplexHeatmap はゲノム解析の可視化における更なる検討と応用の場を広げていくことが期待される。
方法
ComplexHeatmap パッケージはオブジェクト指向プログラミングパラダイム (object-oriented programming) を採用し、R 言語および Bioconductor フレームワーク上に実装された（Bioconductor パッケージ URL (Uniform Resource Locator): http://www.bioconductor.org/packages/devel/bioc/html/ComplexHeatmap.html）。主要クラスは 3 つである。
Heatmap クラスは単一ヒートマップの表現を担い、データ行列を処理してヒートマップ本体・行列名・タイトル・デンドログラム・列アノテーションの描画メソッドを提供する。HeatmapList クラスは複数の並列ヒートマップのリストを管理し、全体のグラフィカル設定調整とレイアウト生成を行う。+ 演算子をオーバーロードすることで、対応する行が同一オブジェクトを指す並列ヒートマップの連結（例：Heatmap(matrix1, ...) + Heatmap(matrix2, ...)）を宣言的に記述できる。HeatmapAnnotation クラスはヒートマップの列または行に整列するグラフィックを汎用的に表現し、定型タイプ（点・棒グラフ・ボックスプロット）とユーザー定義グラフィックの両方をサポートする。列アノテーションは top_annotation 引数で単一ヒートマップに付加し、行アノテーションは + 演算子でリストに連結する形式をとる。
統計手法として、データのグループ化には階層的クラスタリング (hierarchical clustering) が用いられ、dendextend パッケージ (Galili 2015) を利用したデンドログラム描画が可能である。ユーザー定義距離関数としてピアソン相関係数 (Pearson correlation) を含む任意の 2 ベクトル受入関数を指定でき、発現プロファイルの類似性計算に利用できる。データのサブグループ分割には k-means クラスタリングまたはデータフレームによるパーティショニングが利用可能である。
パッケージの実証には n=4 の実データセットが用いられた。(1) Cancer et al. Nature 2014 の TCGA 肺腺癌コホート：n=134 患者、n=38 選択遺伝子のゲノム変化を拡張 OncoPrint として可視化。(2) Buettner et al. (2015) のマウス T 細胞単細胞 RNA-seq：n=721 選択遺伝子の発現プロファイルによる細胞不均一性の解析。(3) 未発表研究のパターンに基づきランダム生成したメチル化・遺伝子発現・エンハンサー相関データ。(4) Sturm et al. (2012) のグリオブラストーマ (glioblastoma) メチル化プロファイリング Figure 1 の再実装と新規アノテーション層追加。各データセットに対してヒートマップ構成の違いを生かした可視化アプローチを設計し、従来ツールでは困難だった多層的解析を実証した。