マイクロアレイデータの解析例 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 のタイムポイントもあるので確認してみる必要がありそうです。

 

Intensity-based Z-score (シグナル値を考慮した Z-score )

マイクロアレイにおいて、シグナル値の低いものは、信頼性(データの再現性)がよくない場合があります。これは、MAプロットや散布図の形状からも推測されます。(同一条件で比較したサンプルにおいても、シグナル値の低い部分に変動が見られがち。)

例えば、シグナル値が、1000から10000になったら、ratio=10ですが、10から100になっても同じ ratio=10 です。ratio だけで、発現変動遺伝子を判定すると、特にシグナル値の低い部分に、ratio が高い遺伝子が見つかりがちですが、再現性を考えた場合、少し不安です。

一方、 Z-score で判定した場合も、MA プロットから分かるように、シグナル値の高い部分は、Z-score が高くなりにくい傾向があります。

この2点の問題を解消するために、  Intensity-based Z-score (Intensity-dependent Z-score)* というものがあります。(以前から提案されているもので、珍しい手法ではありません。)

Intensity-based Z-score

Intensity-based Z-score では、文字通り、Intensity (=シグナル値の強さ) を考慮した Z-score です。マイクロアレイデータをシグナル値の大きさごとに一定の区間に分割して、その区間ごとに logFC の平均値と標準偏差を求めます。そして、それぞれの区間ごとに、Z-score を算出します。

そうすることで、シグナル値の低い部分は、大きい標準偏差 (SD) を用い、一方、シグナル値の高い部分は、小さい SD を用いて Z-score を判定することになります。したがって、シグナル値の低い部分は、通常の Z-score より判定条件が厳しくなり、シグナル値の高い部分は、判定条件が緩くなります。

MAプロットで確認すると、違いがよくわかります。ここでは、便宜的に、シグナル値の大きさによって8区間に分割して、Z-score を算出しています(図下)。(8つという区間に必然性はありません。もう少し細く区切ってもよいです。)また、Intensity-based Z-score を使用した図は、事前に使用した2サンプルのフラグが、共にA(absent: 未検出)のプローブを除外しています。

Z-score で色付けした MA プロット。
Z-score で色付けした MA プロット。
Intensity-based Z-score で色付けしたMAプロット。(フラグでカットしている。)
Intensity-based Z-score で色付けしたMAプロット。(フラグでカットしている。)

図形的には、変動している遺伝子を外側から数パーセント取るという形になります。単純に ratio で判定するよりも慎重な方法といえるでしょう。しかしながら、これも万能な手法というわけではありません。可能であれば、繰り返しサンプル (replicates) が、3サンプル以上 (n=3以上) あって、検定などの手法が使える方が望ましいことには変わりません。

参考

* Quackenbush. Microarray data normalization and transformation. Nat Genet (2002) vol. 32 Suppl pp. 496-501.

 

Z-score の利点と欠点

MAプロットで確認すると、 ratio と Z-score の関係が、分かりやすいと思います。解析例1のデータ (sample2/control2) においては、ratio=2 のラインより、2SDのラインが外側にあります。この場合は、 Z-score のほうが厳しい判定と言えます。(追記:この例の場合、2SD=約1.344 (logFCで)、ratio > 2.54 が Z-score > 2 のラインです。)

MAプロットに示した Z-score の関係。
MAプロットに示した Z-score の関係。

また、 logFC を Z-score 化して、発現変動遺伝子の判定を行う場合、利点と欠点があります。

Z-score を判定に用いた場合の利点

単純に ratio で判定する場合に比べて、下記のような利点があります。

  • 判定の基準で迷わない。
  • 発現変動が大きいデータの場合、カットオフに使える。
  • 発現変動が小さいデータでも、発現変動遺伝子を判定できる。

Z-score > 2 または、 Z-score < -2 ということに、統計学的な意味があるので、判定の基準をいくらぐらいにしようかということで迷う必要はありません。(2 以下にあまり意味がありません。)

また、発現変動が大きいデータの場合、ratio > 2, ratio < 0.5 などで判定すると、4000個、5000個の遺伝子が、発現変動ありとなることもあります。そのような大きいデータの場合、 Z-score で判定することで、発現変動遺伝子の数を1000個から2000個程度に減らせます。(発現変動遺伝子が多すぎても困ることがあります。例えば、DAVID で1度で解析できる遺伝子の数は、3000個です。)

逆に、サンプルによっては、たかだか1.5倍程度しか変動しないような、小さい発現変動のデータもあります。そのような場合、Z-score だけで判定すれば、(ratioの大きさに関係なく)変動している遺伝子の上位5%程度を取ることになるので、必ず、変動遺伝子を得られます。

Z-score を判定に用いた場合の欠点

利点の裏返しです。

  • 判定が厳しすぎるときがある。
  • シグナル値が高い部分で、|Z-score| > 2 に比較的なりにくい。

例えば、研究対象としている遺伝子があり、 ratio だと 2 倍を超えているのに、 Z-score だと、 1.7 というようなケースです。Z-score  が2以下だからといって、「変動していない」という証明になるわけではありませんが、印象が良くないかもしれません。

また、散布図やMAプロットの先が細いということからも見えることですが、シグナル値の高い部分は、変動が小さく見えがちです。シグナルの高い部分にも、シグナルの低い部分と同程度の標準偏差を求めるならば、高い部分は、|Z-score| > 2 になりにくいといえます。

しかしながら、そもそも、変動が大きすぎたり、変動が小さい場合は、まず、サンプル数(n数)を増やすことを考えるべきです。n=3以上にして、limma, SAM, t-test などの検定を行ったほうが、p-value も算出できますし、より再現性の高い遺伝子を変動ありと判定できます。

 

Z-score (発現変動遺伝子を判定するもう1つの方法)

ratio (logFC) 以外に発現変動遺伝子を判定する方法として用いられるものに、 Z-score があります。あまり聞きなれない用語かもしれませんが、偏差値というとどうでしょうか?

ある値が、その群の平均値から、標準偏差 (SD) の何個ぶん離れているかを求めたものが、 Z-score です。Z-score の考え方自体は、特殊なものではありません。Z-検定 (Z-test) という使われ方もあります。(Z-score 化されるのは、logFC に限りません。例えば、 cBioPortal では、シグナル値が Z-score 化されています。)

ここで用いるのは、 logFC を Z-score 化したものです。ある logFC  が、 logFC の平均値から、標準偏差の何個ぶん離れているか計算します。

ratio (logFC) のヒストグラム

ratio (logFC) のヒストグラムを示します。データは、解析例1のものを使っています。

ratio (logFC) のヒストグラム。
ratio (logFC) のヒストグラム。

ほとんどの遺伝子の logFC が0(つまりratio=1)、変動していないことが分かるかと思います。ratio で判定するなら、 ratio > 2 または、 ratio < 0.5 を満たす部分に含まれる遺伝子を発現変動遺伝子と判定していることになります。(ヒストグラムの右側と左側)。

標準偏差と Z-score

次に、同じヒストグラムに標準偏差 (SD) と Z-score の関係を示してみます。

標準偏差とZ-scoreを示したヒストグラム。
標準偏差とZ-scoreを示したヒストグラム。

解析例1のデータ (sample2/control2) の場合、計算してみると、logFC の平均値は、ほぼ0になります。SDは、約 0.672 です。また、分布関数の考え方から、平均値から標準偏差2個分離れていれば、稀な値 (p-value < 0.05) ということになります。

したがって、 Z-score で判定する場合は、下記のようになります。

  • 増加した遺伝子: Z-score > 2
  • 減少した遺伝子: Z-score < -2

絶対値の記号を用いれば、両方を合わせて、  |Z-score| > 2 とも書けます。

追記:2SDだと、2 x 0.672 = 1.344 なので、 logFC > 1.344 であれば、 Z-score > 2 です。ratio で言い換えると、ratio > 2.54 であれば、Z-score > 2 となります。

参考

https://www.khanacademy.org/math/probability/statistics-inferential/normal_distribution/v/ck12-org-normal-distribution-problems-z-score

 

GSEA 操作ガイド (4): 結果の表示と解釈

正しくパラメーターを設定できれば、GSEA を実行後に結果のレポートを確認できます。そのレポートの読み方を紹介します。

1. 実行状況の確認

Run を押した後、ウィンドウ左下の Status には、”Running” と表示されています。Error が出ずに終了すれば、この表示が、 “Success” に変わります。(選択したリストが多くなければ、数分後には終了します。)

GSEAの終了を確認。
GSEAの終了を確認。

2. レポートの表示

Success” の部分をクリックすると、ブラウザが起動し、解析結果のレポートが表示されます。 続きを読む GSEA 操作ガイド (4): 結果の表示と解釈