マイクロアレイデータの構成
以前の投稿で述べたように、マイクロアレイデータの表示方法には、下記のようなものがあります。
- テーブル形式で表示（エクセル）
- ヒストグラム
- ボックスプロット
- 散布図
- MAプロット
まず、エクセルで開いてみましょう（参考）。マイクロアレイデータは、主に下記の項目から構成されています。（Agilent 社のマイクロアレイを例として用いています。）
- ProbeName (ProbeID): マイクロアレイ上のプローブのIDです。重複しません。例：A_23_P26810, A_33_P3315764 など。
- GeneSymbol: プローブのコードする遺伝子（以下、単に遺伝子）の公式な略称 (Official Gene Symbol) です。例：p53 の場合は、TP53 が GeneSymbol です。他の遺伝子と重複しないように考慮された略称です。
- Description: 遺伝子の機能などの説明文です。例：Homo sapiens tumor protein p53 (TP53), transcript variant 1, mRNA
- サンプル名（シグナル値）: サンプル名で示された列に、シグナル値を表示することが多いです。例：WT, KO など。任意の名称です。
- フラグ：Absent (A), Marginal (M), Present (P) のフラグです。そのプローブによって、対象の遺伝子が検出できたかどうか判定したものです。コールとも呼ばれます。（メーカーによっては、フラグがない場合もあります。）
- GO: 遺伝子に付けられた Gene Ontology (GO) の用語です。複数の用語を羅列することで、遺伝子の機能を表現します。例：GO:0000060(protein import into nucleus, translocation)|GO:0000075(cell cycle checkpoint) など。
- EntrezGeneID (EntrezID): 遺伝子の NCBI における ID です。例：7157
- その他：GenbankAccession, EnsemblID など、その他のデータベースにおける遺伝子の ID があります。
マイクロアレイデータは、1行に1プローブの情報が記録されています。これが、プローブの数だけ、縦に並んだ構造となっています。例えば、4万プローブを搭載したマイクロアレイの結果は、4万行のデータとなります。1行目にヘッダーと呼ばれる項目名の行があります。なお、ヘッダーより前の行に各種のパラメーターや情報が記載されている場合があります。列の並び（横の並び）は、自由なため、上記の項目の順に並んでいるとは限りません。また、行の並びも特に決まった順序がありません。ProbeID の順や、GeneSymbol の順（アルファベット順）、ratio の大きい順などの順で並べ替えられて（ソートされて）いることがあります。
マイクロアレイデータ
このうち、マイクロアレイデータの解析において、最低限必要なのは、「ProbeID」、「シグナル値」の2つの項目です。これ以外の情報はすべて、いわゆる、「アノテーション」と呼ばれる情報です。ProbeIDさえあれば、それに関連づけされた（annotateされた）アノテーションは、すべてデータベースから取得可能です。また、ヘッダーより前の情報も通常、解析に用いることはありません。
マイクロアレイデータ（最小構成）
アノテーションの情報が含まれていると、難しそうに見えますが、マイクロアレイ解析に用いる純粋なデータは、行数が多いものの、数字が並んでいるだけの単純な構造となっています。
マイクロアレイと聞いて、最初に思い浮かぶイメージはどのようなものでしょうか？多くの人が、赤と緑の画像をイメージするでしょう。しかし、この画像は誤解されたイメージかもしれません。
おそらく、イメージは、下の2パターンに別れると思います。
マイクロアレイのイメージ
左側をイメージされた方、正解です。これは、マイクロアレイの2色法のイメージです。比較するサンプルを赤と緑の2色でラベルするので、そう呼ばれています（two color, two channelとも）。一方のサンプルで高い発現を示す遺伝子は、赤または緑に見え、両方のサンプルで発現している遺伝子は、赤＋緑で黄に見え、両方のサンプルで発現していない遺伝子は黒に見えるという原理です。
しかし、これは古いイメージといえるでしょう。ここ数年の間にマイクロアレイの環境は変わり、近年では1色法が主流です。実際に NCBI のマイクロアレイのデータベースである Gene Expression Omnibus (GEO) に登録されているマイクロアレイデータのほとんどが、 Affymetrix 社の GeneChip® であり、 この GeneChip® は、1色法のマイクロアレイです。また、Agilent 社のマイクロアレイも2色法と1色法の両方に対応しています。
1色でどうやって比較するのか？と疑問を持たれる方もいるかもしれません。端的に言って、マイクロアレイ解析とは画像処理ではありません。解析に用いるのは数値となったデータであり、単純に数字の大小を比較しているに過ぎないのです。考え方としては、1色でも濃淡で比較できると考えていただいてもよいかもしれません。（1色法、2色法それぞれにデメリットとメリットがあります。）
マイクロアレイ != ヒートマップ
右の画像をイメージされたかたも、ある意味正解です。正確にはヒートマップですが、マイクロアレイデータの表示方法として、論文によく登場するため、このイメージが強いのかもしれません。
しかし、「1色法なのに、どうして、赤と緑の色がついているのだろう？」と思われた方もいるのではないでしょうか。その場合は、本サイトの「MeV の使い方 3. t-検定（続き）」などをご覧ください。シグナル値の色づけは自由です。コントラストが強いため、赤と緑を選ばれることが多いのでしょう。
色とシグナル値の大きさの関係
また、色の意味が異なることにご注意ください。赤＝発現が高い、緑＝発現が低い、黒＝発現していない、ではありません。「黒＝発現していない」というのは、上記に述べた2色法の画像での話です。ヒートマップにおいて、「黒＝中くらいの値（高くも低くもない）」を意味します。また、ヒートマップの色には「黄（オレンジ）」がないことにもご注意ください。
マイクロアレイの画像とヒートマップは、まったく異なるものです。マイクロアレイの画像を整理して並べ替える（クラスタリング処理する）と、ヒートマップが得られるというものではありません。
これらの2つの用語は、同じような意味で使われます。しかし、厳密には少し違った意味を持っています。「ヒートマップ」のほうが、広い意味を持っていると言えます。「クラスタリング図」は、いわば、クラスタリング処理を施されたヒートマップです。よって、どちらか分からない場合は、ヒートマップと呼ぶほうが無難でしょう。
クラスタリング処理は、簡単にいうと、並べ替えです。ヒートマップ中の個々の遺伝子を、発現パターンの似ている順に並べる作業です。また、並べ方には、方向があります。サンプルを似ている順に並べるか、遺伝子を似ている順に並べるか、というものです。
例：3つのサンプルA, B, C があり、それぞれ、gene1, gene2, gene3 のシグナル値を測定したとします。
|
サンプルA |
サンプルB |
サンプルC |
| gene1 |
30 |
20 |
35 |
| gene2 |
10 |
50 |
15 |
| gene3 |
10 |
20 |
15 |
クラスタリング図（＝クラスタリング処理したヒートマップ）
クラスタリング：サンプルを似ている順に並べ替え
クラスタリング：遺伝子を似ている順に並べ替え
両方向というパターンもあります。
クラスタリング：サンプルを似ている順に並べ替え＋遺伝子を似ている順に並べ替え
ratio を用いた発現変動遺伝子の判定の問題点
前回、ratio を用いた発現変動遺伝子の判定基準を述べました。また、下記の3つの問題点を挙げました。今回は、それらの問題点への対策を述べます。
- しきい値ギリギリの遺伝子。 (1.9倍に増加している場合など)
- シグナル値が低く、 ratio は高い遺伝子。 (10 から 100 に増加した場合など。)
- シグナル値が高く、 ratio が低い遺伝子。 (40,000 から 60,000 に増加した場合など。)
問題1. しきい値ギリギリの遺伝子（1.9倍に増加した遺伝子など）
前回、発現変動遺伝子を判定する条件として、2倍（0.5倍）という、しきい値を挙げました。このとき、1.9倍に増加した遺伝子は、発現変動遺伝子と判定されません。わずかの差で発現変動遺伝子とならないのは、惜しい気がします。この場合、どうしたらいいのでしょうか？
解決策1: 「あきらめる」
解決策の1つは、「あきらめる」です。仮に1.9倍もOKとした場合、1.8倍はどうするのか？1.8倍もOKなら、1.7倍は、、、とキリがありません。必ず、どこかで線引きが必要になります。1.9倍は、なにも発現変動しなかったと言っているわけではありません。いったん、2倍という基準を用いて、その後の解析を行って、それから次点の発現変動遺伝子の候補を探しても遅くはないでしょう。また、発現変動遺伝子の判定基準を低く設定すると、新たな問題が発生します。それは、発現変動遺伝子と判定される遺伝子の数が多すぎるという問題です。発現変動遺伝子の数が数千個となった場合は、その後の解析に支障が出ます。（詳しくは別で述べます。）
解決策2: 「ratio 以外の判定基準を使用」
別の解決策としては、「ratio 以外の判定基準を使用する」というものが考えられます。これは、次の問題点の解決策でもあるので、次で述べます。（p-value や Z-score などがあります。）
解決策3: 「他の実験結果をつける」
論文投稿時であれば、リアルタイムPCRなど「他の実験結果をつける」という方法もあります。特定の遺伝子の発現変動の有無に注目されている場合、その特定の遺伝子の ratio が、しきい値以上になっていなければ気になる方も多いと思います。前回、述べたように、判定基準の2倍という数字に根拠はありません。そのため、ほかの実験で確認できているのであれば、1.9倍であっても、1.5倍であってもレビュワーから問題にされることは少ないです。
解決策4: 「正規化のアルゴリズムを変更する」
別の解決策としては、「正規化のアルゴリズムを変更する」というものです。異なるアルゴリズムを用いて正規化を行えば、正規化のシグナル値は、若干、異なる値となります。うまくいけば、あるアルゴリズムで正規化したときに 1.9 倍の増加であった遺伝子が、別のアルゴリズムで正規化されたときに2倍の増加になるかもしれません。ただ、この方法だと、すべての遺伝子のシグナル値に影響があるため、逆のことも起こります。2倍の増加であったものが、1.9倍の増加になるかもしれません。注目しているすべての遺伝子に都合良く、正規化できるアルゴリズムが存在するとは限りません（遺伝子ごとに正規化方法を変えることはできません）。この方法は、いささか、恣意的な行為と考えられるため、奥の手と言えます。
結論1
いくつかの問題点はありますが、 ratio で判定すること自体がイケナイわけではありません。（問題があることを理解した上で判定するのであれば。）また、1度、決めたら判定基準に従うのが無難です。しきい値ギリギリの遺伝子は諦めるのをお勧めします。必ず、どこかで線引きが必要になります。さらに、最近の論文では、マイクロアレイデータだけでは認められず、ほかの実験（リアルタイムPCRなど）の結果も求められることが多いです。そのとき、ほかの実験で結果が出ていれば、ratio の大きさは、あまり問題にされません。いずれにしろ、ほとんどのマイクロアレイデータでは、相当数（数百から千個程度）の遺伝子が発現変動遺伝子と判定されますので、まずは、それらの明らかに変動している遺伝子を中心に結果を見ていくのがよいでしょう。（あえて、微妙な変動を示す遺伝子にこだわる理由があれば別ですが。）
つづきます。
今回も、引き続き、マイクロアレイ解析結果をLinuxコマンドを用いて簡単に操作する方法を紹介していきたいと思います。今回は、端末操作を楽にするコマンドを紹介したいと思います。
端末操作を煩わしく感じるのは、ディレクトリを行き来する際や、同じようなコマンドを何度も入力する時と思います。しかし、ちょっとした工夫で操作が楽しくなるので、活用頂けると幸いです。
使い方① （ディレクトリ移動）
pushd 記憶したいディレクトリ（例：. ~/work, /usr/bin 等）
popd 記憶されたディレクトリに戻る
サンプル＞
$ work # 現在作業中のディレクトリ
pushd . # pushdで、作業ディレクトリを記憶します。（. はカレント）
$ cd ../temp # 別のディレクトリに移動しました。
# 移動後、色々なコマンドを使って作業を行います。
popd # workディレクトリに戻りたい時、popdを使います。
$ work # 前回記憶した作業ディレクトリに戻ります。
使い方②（以前使用したコマンドを呼び出す）
コマンドプロンプトに対して、[Ctrl を押しながら、 r]
サンプル＞
$ work # 作業中に、以前のコマンドを編集したり、再利用したい場合
# [Ctrl + r]を押すと、コマンド候補検索用の表示に変わります。
(reverse-i-search)`':
# ヒントを入力すると、候補コマンドを表示します。
# 以下は、cuと入力した例です。コマンド群が表示されました。
(reverse-i-search)`cu': cut -f1-3 file1.txt | sort > check/cut_file1.txt
# 候補が表示された時、リターンキーで、候補が実行されます。
# カーソルを右や左に移動すると、編集が可能です。
# 編集が完了したら、リターンで実行します。
発現変動遺伝子
マイクロアレイデータの正規化の次に行う処理は、「発現変動遺伝子」の抽出です。「発現変動遺伝子」とは、具体的には、シグナル値をコントロールと比較して上がっている（または下がっている）遺伝子のことです。
例えば、ある遺伝子のプローブについて、コントロールとサンプルのシグナル値が下記のようであったとします。
| 遺伝子名 |
コントロール |
サンプル |
| A |
100 |
200 |
コントロールのシグナル値が100であるのに対して、サンプルのシグナル値は200です。直感的にサンプルのほうが高いことが分かります。どれくらい大きいかを議論するには、「差」を用いる方法と、「比」（何倍か）を用いる方法があります。マイクロアレイデータでは、「比」を用いることが多いです。この場合、比＝サンプル／コントロールを計算します。結果は、200/100 = 2 となります。つまり、遺伝子Aに関して、サンプルのシグナル値はコントロールの2倍に上がっていることになります。「比」は、ratio または fold-change と表記されます。
では、何倍に上がっていれば（下がっていれば）、発現変動遺伝子といえるのでしょうか？
発現変動遺伝子の定義
多くの研究者にとって、マイクロアレイ解析に期待するのは、「どの遺伝子の発現が上がっていて、どの遺伝子の発現が下がっているのか」ということだと思います。つまり、どれが発現変動遺伝子なのかということですが、単純そうに見えて、残念ながら、それに対する明確な定義は存在しません。しかしながら、多くの人に用いられている判断基準というものは存在します。それは下記のようなものです。
- 上がっている遺伝子（発現増加した遺伝子）： ratio が2以上。つまり、2倍以上増加した遺伝子。
- 下がっている遺伝子（発現減少した遺伝子）： ratio が0.5以下。つまり、半分以下に減少した遺伝子。
なぜ、”2″ 倍なのかということに対する答えはありません。（科学的または統計的な根拠のある数字ではありません。）とりあえず、2倍以上であれば、「上がっている」ということに対して、多くの人が賛同できるであろうという見込みの数字です。ある研究者が、1.5倍でも上がっていると主張しても、この点に関しては全く問題ありません。
* ratio (fold-change), logFC
論文中において、シグナル値の比は、ratio または fold-change と表記されます。fold-change は省略されて「FC」とも表記されます。log変換されたシグナル値の比が、 “logFC” です。
問題点
ratio を求めるだけであれば、マイクロアレイデータの計算は単純です。データが複雑に感じられるかもしれませんが、遺伝子Aのようなデータを遺伝子の数だけ（例えば4万個）集めたのがマイクロアレイデータです。ratio を求める割り算を4万回繰り返すだけです。
マイクロアレイ解析の難しさは、それ以外のところにあります。2倍以上という基準を用いたとき、1.9倍はどう考えたらよいのでしょうか？また、シグナル値が10から100になったとき、ratioは10ですが、一方で、シグナル値が1,000から10,000になったとき、ratioは同じ10です。同列に扱ってもよいものでしょうか？また、シグナル値が、40,000 から 60,000 になった遺伝子はどうでしょうか？ ratio としては 1.5 ですが、20,000増加しています。
つづきます。
エクセルで値を log2 変換する方法の一例です。LOG 関数を使います。
LOG 関数 LOG(数値,底)
セル A1 に変換したい値がある場合、セル B1 に =LOG(A1, 2) と入力します。
セル B1 に 2 を底としたセル A1 の値の対数が求められます。
底の指定がない場合 (=LOG(A1)) は 10 を底とした対数が求められます。
今回も、引き続き、マイクロアレイ解析結果をLinuxコマンドを用いて簡単に操作する方法を紹介していきたいと思います。今回は、行の整列に関する操作を扱いたいと思います。
データの昇順／降順の並べ替えは、データのチェック時に良く使用すると思いますが、
他の機能と組合わせて使うと効果なチェックを行ったり、形式を整えたりする事ができます。
使い方
sort data.txt または cat data.txt | sort
# タブ(\t)区切りのテキストファイルを、先頭列をキーにして昇順ソートします。
# 他にも、キーに使用するフィールドの指定や昇順／降順の指定などができます。
サンプル＞
$ cat sorttest.txt
orange 100
apple 150
banana 90
grape 300
melon 400
apple 150
# 単純なソート
$ cat sorttest.txt | sort
apple 150
apple 150
banana 90
grape 300
melon 400
orange 100
# 重複チェックと組み合わせ１（重複を除外）
$ cat sorttest.txt | sort | uniq
apple 150
banana 90
grape 300
melon 400
orange 100
# 重複チェックと組み合わせ２（重複のみ抽出）
$ cat sorttest.txt | sort | uniq -d
apple 150
例：t-検定の結果
- 左側の “Analysis Results” に “T Tests” の結果が表示される。
- “Significant Genes” に含まれる遺伝子群が、 t-検定の結果、有意となった遺伝子群。（標準の設定では、p-value < 0.01 で有意。0.01 または 0.05 を用いることが多い。）
- サンプルデータの場合、有意となる遺伝子はない。すべての有意でない遺伝子は、 Non-significant Genes に表示される。
例：t-検定の結果
3.(3) 検定結果の色づけ
- 現在の色づけは、他の遺伝子と比べている（図の赤枠）。遺伝子AとBが遺伝子Cより低いため、遺伝子AとBが緑になり、遺伝子Cが赤になっている。
- 他のサンプルと比べて、シグナルの高いところは赤、低いところは緑になって欲しい。図に青枠で示されたように、遺伝子Cの6サンプル間での差を見たい。
- 色づけのための調整が必要。調整の手段は、いろいろ考えられる。各遺伝子の中央値からの距離に変換してもよいし、各遺伝子の平均値からの距離（SD何個分か）に変換してもよい。ここでは、簡単な方法として、各遺伝子の中央値からの距離に変換する方法を紹介する。
推定結果の色づけ
3.(4) 遺伝子ごとの中央値からの距離を求める
1) “Adjust Data -> Gene/Row Adjustments -> Median Center Genes/Rows” を選択。 この処理は取り消すことができないので注意。元のシグナル値やほかの調整を行うには、データの読み込みから、すべての作業をやり直す必要がある。
2) 再度、“Display -> Set Color Scale Limits” を選択し、色づけのスケーリングを合わせる。
- Lower Limit = -2
- Midpoint Value = 0
- Upper Limit = 2
＊色の濃さは適宜、調節可能。中央値からの距離の場合、Midpoint Value は、0 以外を用いることはない。Lower = -1 , Upper Limit = 1 のように低く設定すると、色が濃くなり強調された状態となる。あまり、低い値を設定すると、どの遺伝子の差も大きいような誤解を与えるので注意。
遺伝子ごとの中央値からの距離
スケーリング後、再度、色づけをやり直した結果
3.(5) t-検定結果の保存
クラスタリング図を右クリックして出るメニューから、 “Save cluster…” を選択すると、画像ファイルとして結果を保存できる。
t-検定結果の保存
「エクセルで遺伝子名が曜日に自動変換される」の続きです。
1. Office ボタンをクリックして「開く」を選択
2. ファイルを選択して「開く」をクリック。
3. 「カンマやタブなどの区切り文字によってフィールドごとに区切られたデータ」を選択して「次へ」をクリック。
4. 区切り文字の「タブ」を選択して「次へ」をクリック。
5. 「GeneSymbol」の列を選択して、列のデータ形式の「文字列」を選択。「完了」をクリック。
6. GeneSymbol が正しく表示されました。このままエクセルブック形式で保存すると、次からはファイルを開いたときに GeneSymbol が曜日に自動変換されることはありません。
