ex302: MA plot

MA plot は2グループの遺伝子発現を視覚化する便利な散布図である。マイクロアレイ解析でもRNAseq解析でも頻用される。

Definition of M & A

R, G ２つのサンプルを比較していると仮定する。

• M: log of the ratio = 発現量の比
  • M = log(intensityR / intensityG)
• A: intensity average of log intensity = 発現量の相乗平均
  • A = log(sqrt(intensityR * intensityG))

intensity はRNAseq の場合、TPMやRPKM/FPKMなどが使われる。

Make an MA Plot

R環境で、''arab2.txt'' のMA plotを書いてみよう。(ex301によってデータを変数datに読み込み済とする)

Q1) m1 vs m2 のMA-plot を書きなさい。

上記のMAの数式をそのままRコマンドに変換すると以下のようになる。

M <- log2(dat$m2 / dat$m1)
A <- log2(sqrt(dat$m1 * dat$m2))
plot(A, M)

ただしこのままではエラーが発生する。Rのバージョンによっては途中で止まってしまう。バージョンによっては、一応それらしいプロットは描かれれる。この式の何が問題だろうか。３つの問題がある。

• 0で割り算が出来ないにもかかわらず、分母が0のケースがある。
• log(0)は計算できない。
• かけ算による巨大な整数の桁溢れ（オーバーフロー）

これらのエラーに対応するために、式を変換し（早め早めに対数をとることと、+1を加えlog0になるのを防ぐ）、スクリプトを修正。

M <- log2(dat$m2 + 1) - log2(dat$m1 + 1)
A <- 1/2 * (log2(dat$m2 + 1) + log2(dat$m1 + 1))
plot(A, M)

コメント：edgeRにはより気の利いたMAplotを描画するコマンドが定義されている。ほかのRNAseq解析用パッケージやマイクロアレイ解析用パッケージにも同様のコマンドが用意されている場合が多い。ただ、MAプロットは自力で作製できるようにしておきたい。

plotのオプションを変更し、見栄えをよくする。

plot(A, M, pch=16, cex=0.4, ylim=c(-8,8), main="MA plot (m1 vs m2)")

発展：MA plotに、発現比が1, 2, 1/2 を示す線分を追加する。

(例)

abline(h=log2(2), col="red", lty=2)
abline(h=log2(1/2), col="red", lty=2)
abline(h=0, col="blue", lty=2)

（練習）上の練習ではm1 vs m2 の比較のMAプロットを描画した。次にm1 vs h1 の比較をしてみよう。m1 vs m2 はbiological replicateであり、m1 vs h1 はコントロール実験とバクテリア接種実験の比較なので、後者の方がバラツキが大きいはずである。MAプロットでその差を確認しよう。

（解答例）

par(mfrow=c(1,2))
plot(A, M, pch=16, cex=0.4, ylim=c(-8,8), main="MA plot (m1 vs m2)")

A.m1h1 <- 1/2 * (log2(dat$m1 +1) + log2(dat$h1 + 1))
M.m1h1 <- log2(dat$m1 +1) - log2(dat$h1 +1)
plot(A.m1h1, M.m1h1, pch=16, cex=0.4, ylim=c(-8,8), main="MA plot (m1 vs h1)")