ggplot2 を使って MA プロットを作成

2024年1月25日 00:19

遺伝子発現の変動を表示する方法として、ボルケーノプロットを紹介しましたが、もう1つよく使われるプロットとして、MAプロットがあります。もともと、マイクロアレイのデータを表示するのに用いられていましたが、RNA-seq のデータにも用いられます。

MA プロットに用いるデータ

遺伝子発現の変動をスコア化したデータは、logCPM, logFC, p-value などの値を含みます。このうち、MA プロットに用いるデータは、 logCPM と logFC です。（一方、ボルケーノプロットは、 logFC と p-value）

> de_result <- read_tsv("./results/edgeR.DE_results")
> de_result
# A tibble: 27,882 × 5
   GeneID      logFC logCPM   PValue      FDR
   <chr>       <dbl>  <dbl>    <dbl>    <dbl>
 1 MTND1P23     9.06   7.65 3.01e-51 8.39e-47
 2 CD177       -9.98   8.26 1.74e-32 2.43e-28
 3 CDH3         7.27   5.80 1.44e-31 1.34e-27
 4 CPNE7        7.21   4.76 9.51e-25 6.63e-21
 5 RP11-750B16  6.26   5.83 1.27e-24 7.08e-21
 6 ETV4         5.67   5.82 7.52e-24 3.50e-20
 7 MT2A        -4.63   7.07 7.08e-23 2.82e-19
 8 CLDN2        8.57   6.92 2.16e-22 7.54e-19
 9 SLC6A19     -8.82   5.15 3.64e-22 1.13e-18
10 BEST4       -7.08   5.66 3.12e-19 8.70e-16

logCPM

logCPM は、リードをカウントした値 CPM (count per million) の平均値を log2 変換したものです。（明記されないことが多いですが、一般的に底は 2 です。）平均値 = Average なので、 MA プロットの A です。x軸に用います。

logFC

logFC は、log2-fold change です。fold change (ratio) を log2 変換したものです。これも底は2です。複数のサンプルからなる2群を比較しているのであれば、それぞれの群の平均値どうしを比較して、fold-change を算出します。y軸 (M) に用います。（なぜ、こちらがMか、はっきりした説明を知りません。）

MAプロットのコード

散布図の x 軸に logCPM を、 y軸に logFC を指定します。 geom_point() を用いて、コードは下記のようになります。

g <- ggplot(de_result, aes(logCPM, logFC))
g + geom_point()

MA プロットは、下図のようになります。

MAプロットの広がりで、変動の大きさの概要を掴むことができます。logFC > 5 のものも多く、今回用いたデータは、変動が大きいデータと言えそうです。（変動が大きすぎると、あまりきれいなデータでないことが多いです。RNAのクオリティ的に。）