volcano plotとは、RNA-seq解析をした際に、発現比較解析の結果を表すグラフです。
横軸に発現比(log2)、縦軸にはp値(-log10)をとります。

R言語でvolcano plotを描く機会があったので、備忘録としてまとめたいと思います。

(2021年5月9日追記) ggplot2で書く方法はこちらから

本稿では、発現比較解析を行い発現比とp値はすでに算出している前提で、グラフを作ります。
なお、R studio (ver. 1.3.1093)、R (ver. 4.0.3)を用いました。

使用したスクリプトは以下の通りとなります。発現比はfold changeとしています。

# R. User (shibanosuke) 27 January, 2021
#This script is for making volcano plot

#Import dataset
library(readxl)
RNA-seqdata <- read_excel("RNA-seqdata.xlsx")
View(RNA-seqdata)  

#Fold change(log2)の代入
m <- RNA-seqdata$Fold.Change

#p-valueの代入
p <- RNA-seqdata$p-value

#p-valueを-log10で対数変換し、代入
logp <- -log10(p)

#p-valueの閾値設定の数値を代入
pval <- -log10(0.05)

#volcano plotの作成
#ポイントの形は20(黒丸)、大きさを0.5とした
#発現が減少したものを青、増加したものを赤とした
#閾値はp値が0.05未満、fold change(log2)が1以上 or -1以下
plot(m, logp, pch = 20, cex = 0.5, main = 'vorcano plot', xlab = 'log2(Fold Change)', ylab = '-log10(p-value)', col = ifelse(logp <= pval,'black',ifelse(m >= 1,'red',ifelse(m <= -1,'blue','black'))))

#境界部分に、darkgreenの点線を引いた
#p値が0.05、fold change(log2)が1 or -1
abline(v = 1, col= 'darkgreen', lty=2)
abline(v = -1, col= 'darkgreen', lty=2)
abline(h = pval, col= 'darkgreen', lty=2)

excelファイルからのimportは、R studioのimport機能を用いました。
excelファイルをimport後、fold change値とp値をそれぞれ代入し、p値は対数変換(-log10)します。

volcano plotの作図はplot関数を用い、閾値はp値が0.05未満、fold change(log2)が1以上もしくは-1以下としました。また、発現量が増加したものは赤色、減少したものは青色、それ以外は黒色でプロットすることで、分かりやすくしました。
閾値の境界部分を分かりやすいように、abline関数で点線を引き、volcano plot完成です。

最後に、今回使用したスクリプトで作図したvolcano plotを紹介します。
使用データは私が取ったデータを一部改変したものとなります。
volcano plotを作図することで、発現変動遺伝子の分布が分かりやすく可視化されたと思います。