Volcano Plot を綺麗に描きたい - アメリエフの技術ブログ

こんにちは、データ受託解析チームの tsuyuh です。

遺伝子の発現変動解析をしていると良く登場する Volcano Plot ですが、Excel や R で単に散布図を描くだけでは味気ないなと思ったことはありませんか？

今回は、簡単に格好いい Volcano Plot を描く方法をご紹介します。
以前の記事 (Seurat や Monocle のプロットの配色を調べる) に引き続き、"できるだけ作図にこだわりたい" シリーズです。

↑ 最終的にこんな図を描きます。

作業環境

R v4.1.1
Rパッケージ EnhancedVolcano v1.12.0
Rパッケージ gridExtra v2.3
Rパッケージ grid v4.1.1

↑ それぞれのドキュメントに従ってインストールしてください。

使用データ

まず、遺伝子の二群間比較解析の結果を用意します。
今回は、ヒト肺線維芽細胞のRNA-seqデータ SRP012607
(Trapnell C et al., "Differential analysis of gene regulation at transcript resolution with RNA-seq.", Nat Biotechnol, 2013 Jan;31(1):46-53)
を使って、Rパッケージ edgeR で二群間比較解析をおこないました。

その結果はこんな感じ↓

> head(final_res)
             logFC   logCPM        PValue           FDR
RASA2     2.071255 6.048938 9.881313e-324 2.490091e-321
OLFM2     2.856658 5.644147 6.818106e-322 1.686098e-319
MCC       1.875031 6.354475 2.636776e-317 6.395196e-315
PAG1      2.059575 5.926258 6.548451e-317 1.558284e-314
NIPSNAP1 -2.224537 5.383882 1.680666e-316 3.925288e-314
CSRP1    -2.346805 8.455595 7.743112e-316 1.775566e-313
......

スクリプト

(参考) bioconductor.org

以下のスクリプトはあくまで一例 (最近の私のお気に入り) です。
EnhancedVolcano には沢山のオプションが存在するため、自由自在にカスタマイズができます。

# 必要なパッケージの読み込み
> library(EnhancedVolcano)
> library(gridExtra)
> library(grid)

# PDFに書き込みたい
> pdf("VolcanoPlot.pdf", width = 15, height = 10)

# 左の図 (遺伝子名のラベルなし) を作る
> p1 <- EnhancedVolcano(final_res,
        lab = rownames(final_res),
        x = 'logFC',
        y = 'FDR',
        xlab = bquote(~Log[2]~ 'fold change'),
        ylab = bquote(~Log[10]~ 'FDR'),
        pCutoff = 5*10^(-10),
        FCcutoff = 1.0,
        pointSize = 2.0,
        labSize = 0,
        colAlpha = 3/5,
        legendLabSize = 12,
        legendIconSize = 4.0,
        legendPosition = 'top',
        drawConnectors = FALSE)

# 右の図 (遺伝子名のラベルあり) を作る
> p2 <- EnhancedVolcano(final_res,
        lab = rownames(final_res),
        x = 'logFC',
        y = 'FDR',
        xlab = bquote(~Log[2]~ 'fold change'),
        ylab = bquote(~Log[10]~ 'FDR'),
        pCutoff = 5*10^(-10),
        FCcutoff = 1.0,
        pointSize = 2.0,
        labSize = 4.0,
        colAlpha = 3/5,
        legendLabSize = 12,
        legendIconSize = 4.0,
        legendPosition = 'top',
        drawConnectors = TRUE,
        widthConnectors = 0.5,
        colConnectors = 'black')

# 2つの図を並べる
> grid.arrange(p1, p2,
            ncol = 2,
            nrow = 1,
            top = textGrob("ConditionA_vs_ConditionB",
            gp = gpar(fontsize = 32)))

> dev.off()

EnhancedVolcano() 関数のオプション説明
・lab：各プロットの名前 (遺伝子名) のベクトル
・x：入力データのうち x軸に使用する値の列名
・y：入力データのうち y軸に使用する値の列名
・xlab：x軸のラベル
・ylab：y軸のラベル
・pCutoff：y軸における色分けの基準値 ※今回は見やすくするために基準値を極端にしています
・FCcutoff：x軸における色分けの基準値
・pointSize：プロットの大きさ
・labSize：ラベルの大きさ
・colAlpha：プロットの不透明度
・legendLabSize：凡例の文字の大きさ
・legendIconSize：凡例の図形の大きさ
・legendPosition：凡例を書く位置
・drawConnectors：各プロットとそのラベル (遺伝子名) を線で結ぶか否か
・widthConnectors：各プロットとそのラベルを結ぶ線の太さ
・colConnectors：各プロットとそのラベルを結ぶ線の色