マイクロアレイデータの解析例 1.4 (発現変動遺伝子、機能解析)

正規化後のシグナル値から、ratio と Z-score を算出して、発現変動を判定します。その結果、発現変動している遺伝子のリスト(例えば、エクセルファイルのようなテーブル)が得られます。ここでは、 sample2/control2 = 24hr どうしを比較したデータを使用しています。

発現変動遺伝子のテーブル。
発現変動遺伝子のテーブル。

発現変動遺伝子

この例では、下のMAプロットに示されたような部分が変動していると判定されています。判定条件としては、 ratio > 2 かつ Z-score > 2 (または ratio < 0.5 かつ Z-score < -2)という少し厳しめの条件にしています。増加しているもの、減少しているもの合わせて、1032個を変動していると判定しています。

発現変動遺伝子を色づけしたMAプロット。
発現変動遺伝子を色づけしたMAプロット。

まずは、変動している遺伝子の機能を調べてみましょう。機能を調べるには、データベースの DAVID が便利です。

DAVIDによる機能解析

DAVIDに変動している遺伝子の ID のリストをアップロードすると、そのリストに含まれる遺伝子にどのような機能の遺伝子が多いのか(膜タンパクが多いのか、転写因子が多いのか、アポトーシスに関連する遺伝子が多いのか)、確認することができます。操作方法は、リンク先を参照してください。

DAVID の Functional Annotation Clustering を実行すると、下記のようなテーブルが得られます。

DAVID の Functional Annotation Clustering の結果。
DAVID の Functional Annotation Clustering の結果。

左上の表示を確認すると、1032個の変動遺伝子のうち、558個が認識されたことがわかります。(ここでは、EntrezGeneID をアップロードしています。)変動しているもののうち、アノテーションがついていないもの、機能が未知のものが多かったようです。

テーブルの最初 (Annotation Cluster 1) から、signal peptide, glycoprotein 、つまり、シグナルを伝達する分子、糖タンパクの遺伝子が多かったことが分かります。だいたい、この機能のグループ(クラスター)は、どんなデータでも上位に来ることが多いです。これらの用語は、漠然とした分類のアノテーションなので、あまり、情報としての価値はないかもしれません。

アノテーションのクラスターとして、次は、ホルモンのようですが、遺伝子数はそれほど多くありません。

入手元のデータのサマリーによると、変異の結果、アポトーシスが活性化され、ケラチノサイトの分化が抑制とあります。

TRPV3 dysfunction may increase apoptotic activity, inhibit keratinocyte differentiation and disturb the intricate balance between proliferation and differentiation state of keratinocytes in the skin.

ブラウザの検索機能で、”apoptosis” を検索して見ると、テーブルの下の方に見つかりました。

apoptosis を含む Annotation Cluster.
apoptosis を含む Annotation Cluster.

確かに変動しているようですが、数は少ないようです。(1032個のうち、20個程度)Enrichment Score も 0.37 と低いです。

“keratinocyte” を検索すると、 “keratinocyte differentiation” も見つかりますが、これも数は少ないです。

"keratinocyte differentiation" を含む Annotation Cluster.
“keratinocyte differentiation” を含む Annotation Cluster.

24hr の比較では、あまり、変異体の影響が出ていないように思われます。ほかに 16hr, 40hr のタイムポイントもあるので確認してみる必要がありそうです。

 

マイクロアレイデータの解析例 1.3 (散布図、正規化の影響)

ボックスプロットを確認したら、次は散布図 (scatter plot) も確認してみましょう。

散布図

例として、 control2 と sample2 を比較した場合の散布図を示します。 ratio > 2 のプローブ(=遺伝子)を、 ratio < 0.5 の遺伝子をに色づけしています。正規化後データを用いています。

正規化後データの散布図。
正規化後データの散布図。

散布図の広がり方から、平均的なデータのように見えます。(がんのサンプルや、変動が大きいデータでは、もっと点が全体に散らばって見えます。)

正規化の影響

rawとして、正規化前のデータの散布図も作成しました。

raw データの散布図。シグナル値の分布に偏りが見られる。
raw データの散布図。シグナル値の分布に偏りが見られる。

よく見ると、rawデータと、正規化後のデータでは、(ratio で判定した場合)変動ありと判定される遺伝子に異なる部分があることがわかります。特にシグナル値の高い部分です。

散布図の左下から右上に引かれた赤線は、 y=x を意味しています。raw データの散布図は、集団が y=x より下に膨らんで見えています。そのため、rawデータをそのまま用いると、sample2で減少した遺伝子が多く見つかり、増加した遺伝子は少なく見つかることになります。

一方、正規化後データの散布図では、点の中心が y=x 上に載っていることが分かります。(=ほとんどの遺伝子が変動していない。偏りがない。)

このように散布図を確認すると、raw, 正規化前のデータに偏り(バイアス)がないか、また、正規化後のデータから偏りが解消されているのか確認できます。

スコアだけで判断して、結果を誤って解釈しないよう、散布図を必ずチェックするようにしましょう。

 

マイクロアレイデータの解析例 1.2 (ボックスプロット、正規化)

ダウンロードしたデータを用いて、ボックスプロットを作成してみましょう。マイクロアレイデータは、まず、ボックスプロットや散布図を書いて、シグナル値の状態を確認することをおすすめします。

シグナル値のばらつきが極端に大きい場合は、データのクオリティが良くない(サンプルのコンディションが悪い、RNAの分解が進んでいる)ことも考えられれます。

raw データのボックスプロット

raw データ(正規化前)の散布図を示します。便宜上、サンプルの名前を、293T_16hr_Control = control1, 293T_16hr_muTRPV3 = sample1 のように変更しています。サンプルによって、多少上下していることが確認できます。サンプルのクオリティは悪くないように見えます。

raw データのボックスプロット

使用した raw データの値は、下記よりダウンロードできます。

https://www.dropbox.com/s/igxolub38mrm0po/rawdata.txt.zip?dl=0

(Agilent の raw データのファイルからシグナル値だけを取り出す方法については割愛します。)

正規化後データのボックスプロット

続いて、正規化後データ (normalized data) のボックスプロットです。ここでは、正規化のアルゴリズムとして、 quantile 法* を用いています。データの分布がそろっていることが確認できます。

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

使用した正規化後データの値は、下記よりダウンロードできます。

https://www.dropbox.com/s/2m1poolbqz9vizh/normalized_data.txt.zip?dl=0

極端に分布の異なるサンプルは、正規化の際(アルゴリズムによっては)、他のサンプルのシグナル値にも影響を与えることもあります。物理的なサンプルの状況を確認して、RNAの分解などが疑われる場合は、そのサンプルを除外して正規化を行ったほうがよいでしょう。

* Bolstad et al. A comparison of normalization methods for high density oligonucleotide array data based on variance and bias. Bioinformatics (2003) vol. 19 (2) pp. 185-93.

 

マイクロアレイデータの解析例 1.1 (データのダウンロード)

GEO に公開されたデータを利用して、実際に解析を行ってみましょう。マイクロアレイ解析によって、どのような情報が得られるのかイメージがつかめると思います。

マイクロアレイデータ

まずは、 GSE53614 のデータを使ってみます。293T cell lines を用いた解析のようです。サンプル名から推測するに、3つのタイムポイント (16hr, 24hr, 40hr) があり、それぞれに、コントロールと対照サンプルがあるようです。また、使用されたアレイは、 Agilent SurePrint G3 Human GE 8x60Kv2 です。

GSE53614

マイクロアレイデータの取得

データを取得します。取得する方法は、いくつかあります。著者らが正規化したデータをそのまま利用するには、 “Series Matrix Files(s)” をダウンロードしましょう。正規化からやり直す場合は、その下の GSE53614_RAW.tar をダウンロードします。 (ここでは、正規化前のデータを用いました。)

マイクロアレイデータの取得

なお、 SOFT形式のデータをダウンロードすれば、アノテーションが付いた状態のファイルを取得できますが、計算時に読み込むときに苦労する場合があります。また、MINiML形式の場合も、XML形式のデータであるため、計算に用いない情報が含まれています。

Series Matrix File の始めの行(ヘッダー)は、サンプル情報です。 ID_REF から始まる行以降のデータが、正規化後のデータです。1列目はプローブのIDです。2列目以降に各サンプルのシグナル値があります。(6サンプルなので6列)プローブの ID が A_ で始まらないものは、バックグラウンドなどの値です。

サンプル情報がヘッダーに含まれています。
サンプル情報がヘッダーに含まれています。
 

GEO と NCBI アカウント

2012年12月より、GEO のアカウントが、 NCBI のアカウントに統合されたようです。現在、GEO のサイトからログインを試みると、下図のような注意書きが表示されます。

GEO アカウントから NCBI アカウントへ。
GEO アカウントから NCBI アカウントへ統合。

すでに、新しく GEO アカウントを作成することはできなくなっており、代わりに NCBI のアカウントを作成することを求められています。 NCBI のアカウントを持っている人は、GEO のアカウントとリンクさせる必要があるそうです。なお、既存の GEO アカウントによるログインは、2013年3月1日以降、行えなくなる模様です。

NCBI アカウントでは、NCBI を検索した結果を保存できます。これに加えて、GEO へのマイクロアレイデータの登録も行えるようになっています。NCBI アカウントでログイン後、GEOのサイトに移動すれば、これまで通り submit の画面に進めます。

NCBI にログインしているかどうかは、右上にユーザー名が表示されているかどうかで確認できます。ログインしていなければ、 “Not logged in” と表示されます。