マイクロアレイデータの解析例 2.2 (stem cell 関連遺伝子のヒートマップ)

iPS 細胞のデータを正規化したら、何はともあれ、 stem cell 関連遺伝子の発現はどうなっているのかが気になるところだと思います。

このようなときは、アノテーションを元に、 stem cell 関連遺伝子だけを抽出して、クラスタリング処理し、さらにヒートマップで表示するとよいでしょう。(この段階では、ratio や p-value を算出してないので、変動していないものも含めています。)

stem cell 関連遺伝子のヒートマップ

作成したヒートマップの結果を示します。BioGPSのデータを用いたヒートマップでは、 ES細胞に特徴的な遺伝子発現のパターンが見えていましたが、今回のデータではどうでしょうか?

stem cell 関連遺伝子のヒートマップ。

とりあえず、大きく分けると3つのクラスターに分かれているように見えます(あくまで見えているというだけで、以下、決定的な事実というわけではありませんので、注意してください)。

一番左のクラスターには、HepG2 や MCF7 などのがん細胞が含まれています。どうやら、stem cell 関連遺伝子に関しては、 iPS細胞がん細胞で似ていないようです。また、iPS細胞では、がん細胞ほど共通のパターンが見られません。いろいろなバリエーションのiPS細胞がありそうです。

上部のツリーの分かれ方を見ると、iPS 細胞は、さらに大きく2つのグループに分かれるようです。この2つのグループで、はっきり分かれているように見える遺伝子としては、上から3分の1ほどのクラスターに含まれる PSCA, GATA2, CDX2, MSX2, HOXB4 あたりでしょうか。

ヒートマップから読み取れる情報

ヒートマップの色付けは、 median からの距離にしていますので、 色が黒いところ低いというわけではありません。特に、一番下の行のクラスターに含まれている POU5F1 (OCT3) や NANOG は、黒く見えていますが、左側のがん細胞のデータが明るい緑なので、がん細胞よりiPSにおける発現は高いことを意味します。

また、1つのサンプルの遺伝子を1つ取って見ると、それが MCF7 に近いほどであったり、であったりします。個別に見ることには、あまり意味がないかもしれません。例えば、hiPS201B2_p67 の MSX2 や CDX2 の色は、HepG2と同程度の色ですが、これだけを持って何か結論づけることはできません。

MSX2やCDX2の発現が低く見える部分。

繰り返しになりますが、色ではなく、実際の値を確認するほうがよいです。シグナル値が1桁程度のものでも、このような色付け方法では、差が強調されすぎる場合があります。

がん細胞だけではなく、BioGPSのように、いろいろな分化後の組織が混ざっていると、色のパターンが出やすかったかもしれません。(iPS細胞はがん細胞とは違う、ということを確認するという意味ではこれだけでよかったのかもしれませんが。)

 

マイクロアレイデータの解析例 2.1 (iPS細胞のデータ、正規化)

この解析例の紹介では、遺伝子発現のデータはどういうものか実感できるように、いろいろな種類のデータを取り上げたいと思います。

単純なタイムコースのデータに続き、今回は、iPS細胞のデータを取り扱います。

iPS細胞のマイクロアレイデータを取得

例のように GEO から、 iPS 細胞のマイクロアレイデータを取得します。GSE42445 のデータです。様々なバリエーションのヒトのiPS細胞 (hiPS) と幹細胞 (hES)、がん細胞 (MCF7など) のデータです。70サンプルあります。

マイクロアレイデータの正規化とボックスプロット

ここでは、rawデータを取得し、 シグナル値を取り出して、 quantile 法で正規化しました。(正規化後のデータ: Series Matrix File を取得して、そのまま用いても良いでしょう。)

正規化前 (raw) のデータをボックスプロットで表示すると下図のようになります。特に大きく外れたサンプルはないようです。右端のがん細胞のデータが若干低いようですが、よくある程度だと思います。

iPS細胞のデータのボックスプロット。正規化前。

正規化後のデータをボックスプロットで表示すると下図のようになります。データの分布がそろっていることが確認できます。

iPS細胞のデータのボックスプロット。正規化後。

参考

  • Koyanagi-Aoi M, Ohnuki M, Takahashi K, Okita K et al. Differentiation-defective phenotypes revealed by large-scale analyses of human pluripotent stem cells. Proc Natl Acad Sci U S A 2013 Dec 17;110(51):20569-74. PMID: 24259714
 

特定の組織で発現の高い遺伝子を探す(ES細胞だけで発現の高いものは?)

ヒートマップとクラスタリングを応用することで、特定の組織で発現の高い遺伝子を探すことも可能です。

例えば、ES細胞だけで発現の高い遺伝子を探してみましょう。

公開データ (GEO, BioGPS) を利用して、いろいろな組織のデータを取得

まず、データとしては、様々な組織のマイクロアレイデータが必要になります。個々の研究室レベルで揃えようと思うと大変ですが、幸い近年では公開データがあるので、うまく利用しましょう。

いろいろな組織の遺伝子発現を比較したければ、 BioGPS のデータが使えます。心臓、肝臓、肺などの種々の組織、約180サンプルぶんのマイクロアレイデータがあります。

クラスタリングしてヒートマップで表示

データを取得したら、正規化して、クラスタリングを行い、結果をヒートマップで表示します。ここでは、RMA と quantile アルゴリズムによって正規化し、階層的クラスタリングを行いました。

全遺伝子では結果を見るのが大変なので、正規化後、アノテーションを元に stem cell 関連遺伝子だけを抽出しています。また、MeV上で、シグナル値を log2変換して、中央値からの距離に変換しています。クラスタリングは、遺伝子方向、サンプル方向とも階層的クラスタリングで行っています。

結果は下図のようになります。遺伝子が縦に並んでいます(行)。組織(サンプル)が横に並んでいます(列)。

stem cell 関連遺伝子のクラスタリング結果。BioGPSのデータを利用。
stem cell 関連遺伝子のクラスタリング結果。BioGPSのデータを利用。

ES細胞だけで発現が高い遺伝子

左から3分の1ほどの列に、ES細胞のデータ(4列)が並んでいます。これらの4列だけで、赤くなっている遺伝子があれば、それらが、ES細胞だけで発現が高い遺伝子です。

探してみると、真ん中より少し上にある遺伝子のクラスターに、ちょうど4列だけ赤い部分があります。ハイライトすると下図の部分です。遺伝子名を確認すると、Pou5f1 = Oct3, Nanog など有名な遺伝子が該当しているのが分かります。(見やすいように横方向を縮小しています。)

ES細胞だけで発現の高い遺伝子をハイライト。
ES細胞だけで発現の高い遺伝子をハイライト。

なお、色付けの設定から、「黒=発現していない」ではない、ので注意してください。

 

ratio のヒートマップで変動パターンをチェックする

シグナル値ではなく、ratio でヒートマップを書いて、変動パターンにをチェックすることもあります。

ratio のヒートマップ

シグナル値ではなく、ratio を用いたヒートマップの例を示します。データとしては、解析例1の3つの比較において、いずれかの比較で変動のあった遺伝子(約2678個)だけを用いています。

左から、3つの比較(ratio1, ratio2, ratio3) を並べています。ratio1 = sample1/control1, ratio2=sample2/control2, ratio3=sample3/control3 となります。

色付けは、ratio < 0.5 を明るい、ratio > 2 が明るいになるようにしています。

また、図は全体が見えやすいようにリサイズしています。

ratio のヒートマップ。遺伝子方向にクラスタリングを適用。
ratio のヒートマップ。遺伝子方向にクラスタリングを適用。

シグナル値を用いた場合のヒートマップと異なり、ヒートマップにコントロールは表示されません。3列ともであれば、3つの比較において、共通に増加した遺伝子だと分かります。逆に3列ともであれば、共通に減少した遺伝子です。

変動パターンをチェックしやすい反面、計算していない比較の組み合わせ(例えば、sample2 と sample3 ではどちらが高いか)は読み取るのが難しい場合があります。

参考

 

クラスタリング結果の誤解しやすい部分

クラスタリングを解釈する際に、よく誤解しやすい部分があります。クラスタリングの結果というものは、それほど絶対的なものではありません。

クラスタリングの結果は、どのアルゴリズム(hierarchical, k-means, SOM など)を用いてクラスタリングするかで多少異なりますし、同じアルゴリズムを用いたとしても、パラメーターの設定によっても若干変化します。

また、クラスタリングするときに、全遺伝子を用いるか、変動している遺伝子だけを用いるかでも変わります。

サンプル間の距離の計算方法による違い(パラメーターの違い)

サンプル方向にクラスタリングを行い、似ている順位のツリーを書くためには、各サンプル間の距離を求めることになります。このサンプル間の距離を求める方法は、いくつかあります。例えば、単純に引き算して差を取るだけの場合や、2乗してから差を取る場合、相関係数を用いる場合などです。

MeV では、 Distance Metric Selection の部分で、この距離を求める手法を選択することになります。

距離を求める方法の選択。
距離を求める方法の選択。

このサンプル間の距離の求め方によって、最終的な結果のツリーの分類は多少異なります。解析例1の全遺伝子(5万プローブ)をサンプル方向のみでクラスタリングした結果を示します。

Euclidian Distance を選択してクラスタリングした結果。
Euclidian Distance を選択してクラスタリングした結果。

距離を求める手法に Euclidian Distance を選択した場合、control1 だけ外れて(似ていないように)見えます。

Pearson Correlation を選択してクラスタリングした結果。
Pearson Correlation を選択してクラスタリングした結果。

相関係数(Pearson Correlation)を用いて、クラスタリングした場合、control1だけ外れているようには見えません。一般的には、この相関係数を用いることが多いと思います。

クラスタリングに用いる遺伝子セットによる違い

上述の全遺伝子を用いてクラスタリングした結果に対して、変動している遺伝子だけをクラスタリングした結果を示します。 全遺伝子を対象とした場合と異なり、 sample2 と sample3 が近いように見えています。

いずれかの比較で変動していた遺伝子を用いてクラスタリングした結果。
いずれかの比較で変動していた遺伝子を用いてクラスタリングした結果。

このようにクラスタリングの結果は、必ずしも、絶対的なものではありません。よく論文で図示されていますが、あくまで1例であり、あるパターンに分かれる可能性を示したものにすぎません。(仮説の証明手段ではありません。)

特に全遺伝子でクラスタリングした場合、「sample1 と一番近いのは、 sample2 と sample3 のどちらか?」というような、結果の細部にあまりこだわらないほうがよいでしょう。前述のような計算方法の違いで多少前後する可能性があります。ある意味、当然のことなので、これらの結果に一喜一憂する必要はありません。

(その時、用いたデータセットとパラメーターを使えば、とりあえず、あるパターンに分かれて見えるという程度に理解した方が良いです。結果の解釈には生物学的な裏付けが求められます。)