ほくそ笑む

R言語と統計解析について

RNA-Seq とマイクロアレイでの変動遺伝子の比較方法について(1)

僕はバイオインフォマティクスの会社に勤めているんですけど、統計関係の仕事をしています*1
なので、たまにバイオの人から統計関係の相談を持ちかけられることがあります。
今回、ちょっとおもしろいことを相談されて、頭を悩ませたので、ここで紹介しようと思います。

RNA-Seq v.s. マイクロアレイ

今回相談されたことというのは、RNA-Seq の解析結果とマイクロアレイの解析結果で、同じような結果が出たのかどうか知りたいが、どうすればいいか?というものでした。
もう少し詳しく言うと、同じ組織に対して RNA-Seq とマイクロアレイで発現変動を調査したところ、RNA-Seq で変動がみられた遺伝子が 4000 個、マイクロアレイで変動がみられた遺伝子が 1000 個、両方で共通して変動がみられた遺伝子が 300 個だったそうです。
このとき、300 個というのは多いのか少ないのか、RNA-Seq の結果とマイクロアレイの結果は似ているのか似ていないのか、判定するにはどうすればいいのか、というのが相談の内容でした。
うーん、1000 個以上あって、300 個しか共通する部分がないんですよね? 全然似てないじゃんと思いながら、真剣にこの問題を考えることにしました。

問題を整理

ちょっと問題を整理します。RNA-Seq で調査できる遺伝子の総数は 30000 個、マイクロアレイで調査できる遺伝子の総数は 20000 個とします。このうち発現変動がみられた遺伝子の数は RNA-Seq で 4000、マイクロアレイで 1000、共通で 300 でした。これを図にすると下のようになります。

RNA-Seq では全遺伝子が調査できるので、マイクロアレイで調査できる遺伝子はすべて調査できます。
このとき、共通遺伝子 300 個というのは、多いのでしょうか? 少ないのでしょうか?

ランダムサンプリング

ランダムサンプリングによる手法を考えてみましょう。
まず、RNA-Seq の 30000 遺伝子からランダムに 4000 個取り出します。次に、マイクロアレイの 20000 遺伝子からランダムに 1000 個取り出します。こうして取り出した遺伝子の共通部分はどうなるでしょうか?

とりあえず 1000 回ほどサンプリングしてみた結果が下のグラフです。

縦軸が共通遺伝子数、横軸が頻度です。共通遺伝子数のピークは 130-140 個ぐらいで、最大でも 168 個でした。300 個にはまったく届きません。
つまり、この条件では、共通遺伝子が偶然 300 個以上になるのはかなりの低確率だということがわかりました。なるほどー、意外です。
したがって、実験で得られた共通遺伝子は偶然得られたものではなく、なんらかの共通の背景を持っていて、それが原因でこれほど多くの共通遺伝子を持っているのだということが結論づけられます。
まあつまり、この条件で 300 個の共通遺伝子は多いということです。

ところが

ところが、これで終わらないのが遺伝子解析の難しいところです。
上記の計算では、30000 遺伝子すべてが同じ確率で変動するという仮定を置いています。
しかし、遺伝子には、発現変動しやすいものとしにくいものがあると思います。
どんな実験条件でも、ほとんど発現変動しない遺伝子というのもあると思います。
極端な話、もし、遺伝子のうち 55% が発現変動しないと仮定すると、結果は次のように変わります。

だいたい半分くらいが共通遺伝子 300 個を超えています。
つまり、共通遺伝子が 300 個を超えるというのは頻繁に起こることになります。
これでは、RNA-Seq の結果とマイクロアレイの結果が似ているという結論は出せません。
次回はこの問題について考えたいと思います。

*1:統計の仕事「も」していますといったほうが正しいか