Top of Case Study

DNAのメチル化状態は、遺伝子発現にどのくらい影響を与えるのか (DNAメチル化データ解析のケーススタディ)

Cancer Genome Atlas Stomach Adenocarcinoma (TCGA-STAD) のデータが GDC Data Portal からダウンロードして解析に使えます。さらに、Subio Platform を使えば、RNA-SeqとDNAメチル化アレイのデータを簡単にインポートできます。(インポートの操作は、"TCGAのRNA-Seqデータをインポートする" または "TCGAのDNAメチル化データをインポートして解析する" のチュートリアルをご覧ください。) 下記は、遺伝子発現とDNAメチル化という2種類のマルチオミクスデータを統合解析する実践的なケーススタディです。

このマルチオミクスのデータセットは、396の Illumina HumanMethylation450K BeadChip サンプル (2 Solid Tissue Normal, 394 Primary Tumor) と、406の RNA-Seq サンプル (32 Solid Tissue Normal, 374 Primary Tumor) で構成されます。このうち、336 人の患者さんは遺伝子発現とDNAメチル化の両方のデータがそろっていますので、遺伝子発現パターンにDNAメチル化パターンがどの程度影響を与えているかを見るのに十分なデータと言えるでしょう。

それでは、まずDNAメチル化のデータから見ていきましょう。Subio Platform を使ってTSSプロットを描画してみます。

TSS Plot - histogram

このマイクロアレイに搭載されているメチル化サイトは、TSSの近傍500bp以内に極端に多くあることがわかります。さらにTSSプロットを散布図にすると、もっと多くの情報を視覚化することができます。たとえば、正常サンプルのベータ値をTSSの周辺にプロットしてみましょう。

TSS_BetaValues_Normal

TSSから2000bp以上離れた領域では、ベータ値は安定した二値的分布(つまり、メチル化しているかしていないか)をしているようです。これに対し、TSSの近くでは中間的な値が多く、メチル化状態が不安定になっているように見えます。それでは、腫瘍サンプルではどうなっているか見てみましょう。

TSS_BetaValues_Tumor

ベータ値の分布パターンは似ていますね。正常と比べて若干違っているところもあるようですが、その差は上の二つの図を比べてもあまりよくわかりません。そこで、正常と腫瘍の差分をTSSプロットで視覚化してみました。

TSS_DiffOfBetaValues

今度はメチル化状態の変化がTSSの近傍で限定的に起こっていることがよくわかります。しかし、その振れ幅は小さいです。上の図は336の腫瘍サンプルの平均値で、平均値は0に収束する傾向があるかもしれません。そこで、腫瘍サンプルを一つ取って正常との差分を見てみました。

TSS_DiffOfBetaValuesOfOne

予想通り、変化が大きくなりました。しかしそれでもメチル化状態は0から1のように大きくシフトするわけではないようです。エピジェネティックの微妙な調整によって遺伝子発現レベルを大きく変えることができるのかもしれません。また、上の図からは別の情報が得られました。このデータでは-0.15から0.15の振幅はノイズのようです。一方、0.2を超える変化は何か意味がありそうです。そこで、336中333サンプル以上で振幅が0.15内に収まるものを意味のある変化をしないものとして除くと、測定されたメチル化サイトのうち1/3を除去できました。残りの2/3を以降の解析に用います。

DNAメチル化状態と、近傍遺伝子の発現パターンの相関係数を計算してみよう

それでは、336人の患者さんのデータで、遺伝子発現とメチル化状態のパターンがどの程度逆相関を示すのか見てみましょう。この解析を行う前に、RNA-Seqのデータからは測定値がノイズの遺伝子を除く処理を行っています。下のヒストグラムが遺伝子の発現パターンと、そのTSSの上流500bp以内CpG island内にあるメチル化サイトの平均値の変化パターンの相関係数の分布を表しています。全体的に左側(-1の方向、つまり逆相関側)に偏っていますが、ほとんどは弱い逆相関関係(-0.3付近)にすぎないことがわかります。ただ、わずかですが強い逆相関(-0.7付近)を示すものもあります。

Met-GX 0to500 Island

上図とは対照的に、microRNAと(予測された)ターゲット遺伝子の発現パターンの間で相関係数を計算して、その分布を図にしたものが下図ですが、これはほぼ正規分布の形をしています。つまりmiRNAはmRNAの量に影響しないということを表しています。このような結果になる理由として、次の3つが考えられるのではないでしょうか。

いずれにせよ二つの相関係数の分布を見比べることで、すべてではなく一部の遺伝子において、メチル化状態と発現量に負の相関関係があることが「単なる偶然ではなく」示唆されていると言えるでしょう。

Corr Bw Gx And Mi Rna

また、この偏った分布形は、上流1500~2500bpCpG islandの外にあるメチル化サイトの平均値との相関係数を計算してみた結果ではなくなっています。紺結果から、発現パターンに強い影響を持つメチル化サイトは、TSS近傍(500bp以内)に限られるということが推測されます。

Met-GX 1500to2500 not-Island

相関係数の分布パターンを、下のマトリックスで見比べてみてください。DNAのメチル化状態が遺伝子発現に対してもつ影響は二種類あるのではないかということが考えられます。まずは強い影響力を持つタイプで、これはCpG islandにあるかどうかは関係なく、TSSから500bp以内のごく近傍にあることが条件のようです。次に弱い影響力を持つタイプですが、これは上流1500bp以内のCpG islandにあるようです。

Correlation Coefficient Matrix

DNAのメチル化状態の強い影響を受ける(相関係数が-0.6以下で抽出)遺伝子は、わずか216個でした。これは、弱い影響を受ける遺伝子と比べてかなり少ないです。しかし、216個中53個がZNFファミリーの遺伝子というのは注目に値するのではないでしょうか。下表はGOエンリッチメント解析の結果で、この結果を見ても転写調節因子に極端に偏っていることがわかります。メチル化状態の影響を強く受ける遺伝子の数は少なくても、これらの特殊な転写因子を介してゲノムワイドな遺伝子発現の変化をもたらしているのかもしれません。

enrichedGO_genes_under_the_strong_effect

これら特殊な転写因子を用いて、ステージが1と2の患者さんを階層型クラスタリングにかけてみた結果が下のヒートマップになります。患者さんが二つのクラスターに分かれています。

clustering of TFs over starge 1-2 patients

初期ステージの患者さんの二つのクラスターの生存曲線解析をしてみると、下図のようになりました。特に2年目に生存率に違いが出てきているように見えます。これらの転写因子の発現レベルが、初期ステージの患者さんの予後に影響しているのかもしれません。対照的に、ステージ3と4の患者さんではこのような影響は見られませんでした。

survival curves of stage 1-2 patients

解析ツールについて

このケーススタディで使用した解析ツールは、このページの右側の "Product" セクションにリストアップされています。解析ツールの詳細はプラグインのページでご覧ください。 

Plug-ins

Subio Platformは、無料で使えるオミクスデータブラウザーおよびデータ共有基盤です。詳細はSubio Platformのページでご覧ください。

Subio Platform

解析ツールの操作になれていない方は、データ解析サービスをご注文いただくことで、Subioが代わって解析を行うよう依頼することができます。詳しくは下のリンクでご覧ください。

データ解析サービス