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

R言語のplot関数を用いてvolcano plotを描く方法を以前記事にしましたが、今回はシミュレーションデータを作成し、そこからggplot2を用いてvolcano plotを描く方法を備忘録としてまとめたいと思います。
前回の記事はこちらから

なお、筆者の用いたR studioはVersion 1.4.1106、Rはversion 4.0.5を用いました。ライブラリのバージョンは以下のスクリプト内に記載しています。

シミュレーションデータの作成

まず、volcano plotを描くためのシミュレーションデータを作成していきます。
RのTCC packageを用いてデータを作成しました。
TCC packageについては、Sunらの論文を参照してください(ページ最下部に記してます)。
シミュレーションデータの作成は以下のURLを参考にしました。

以下のスクリプトを用いて、シミュレーションデータをcsvファイルに書き出します。
今回は、biological replicateを3、2群間比較とし、iDEGES/edgeRで発現変動解析を行いました。

library(TCC)
packageVersion("TCC")
packageVersion("edgeR")

count <- simulateReadCounts(Ngene = 10000, PDEG = 0.20, DEG.assign = NULL,
                            DEG.foldchange = NULL, replicates = NULL, group = NULL,
                            fc.matrix = NULL)
in_f <- count$count

out_f <- "TCC_simulation_data.csv"

param_A <- 3
param_B <- 3
param1 <- 3
param_DEmethod <- "edger"
param_FDR <- 0.05

data <- as.matrix(in_f)
data.cl <- c(rep(1, param_A), rep(2, param_B))
data.cl
tcc <- new("TCC", data, data.cl)  

tcc <- calcNormFactors(tcc,iteration=param1)

tcc <- estimateDE(tcc, test.method=param_DEmethod, FDR=param_FDR)

result <- getResult(tcc, sort=FALSE)

tmp <- cbind(rownames(tcc$count), tcc$count, result)
write.csv(tmp, out_f, append=F, quote=T, row.names=F)

volcano plotの作成

さて、ここからvolcano plotを描いていきます。
まずは、必要なパッケージの読み込みとデータの読み込みを行います。

library(ggplot2)
library(dplyr)
packageVersion("ggplot2")
packageVersion("dplyr")

countdata <- read.csv("TCC_simulation_data.csv")
head(countdata)

次に、volcano plotを描いた際に、発現が上昇した遺伝子、低下した遺伝子それぞれのプロットを色付けするための条件を記述し、新たな変数としてcountdataに追加します。

countdata <- mutate(countdata, thcolor = ifelse(p.value >= 0.05 , 0,
                                   ifelse(m.value >= 1, 1,
                                          ifelse(m.value <= -1, 2,0))))

p値が0.05未満で、fold change(m.valueのこと)が1以上の時は1、fold changeが-1以下の時は2、それ以外は0を指定しています。

いよいよグラフの作成です。

ggplot() +
  geom_point(aes(x = countdata$m.value, y = -log10(countdata$p.value), color = as.factor(countdata$thcolor)), size = 0.1) +
  labs(x = "log2 (Fold Change)", y = "-log10 (p-value)") +
  scale_color_manual(values = c("black", "red", "blue")) +
  geom_hline(yintercept = -log10(0.05), size = 0.2, color = "dark green") + 
  geom_vline(xintercept = -1, size = 0.2, color = "dark green") +
  geom_vline(xintercept = 1, size = 0.2, color = "dark green") +
  theme(legend.position = "none", panel.grid = element_blank())

fold changeが1以上の遺伝子を赤、-1以下のものを青、それ以外を黒色でプロットしました。
また、閾値を分かりやすくするため、深緑色の線を引きました。
volcano plotはp値を対数変換する必要があるので、-log10()で変換しました。

以上のスクリプトを実行し描いたグラフは以下のようになります。

発現変動遺伝子が分かりやすく可視化されました。
plot関数よりもggplot2の方がこだわれば、よりきれいなグラフが描くことができるかもしれません。(少なくとも私自身はplot関数で描くのと、どちらが良いか分かりませんが。。。)

最後に、本記事で用いたスクリプトを実行した結果を以下に示します。

> library(TCC)
~読み込みが長いので割愛~

> packageVersion("TCC")
[1] ‘1.30.0’
> packageVersion("edgeR")
[1] ‘3.32.1’
> count <- simulateReadCounts(Ngene = 10000, PDEG = 0.20, DEG.assign = NULL,
+                             DEG.foldchange = NULL, replicates = NULL, group = NULL,
+                             fc.matrix = NULL)
TCC::INFO: Generating simulation data under NB distribution ...
TCC::INFO: (genesizes   :  10000 )
TCC::INFO: (replicates  :  3, 3 )
TCC::INFO: (PDEG        :  0.18, 0.02 )
> in_f <- count$count
> out_f <- "TCC_simulation_data.csv"
> param_A <- 3
> param_B <- 3
> param1 <- 3
> param_DEmethod <- "edger"
> param_FDR <- 0.05
> data <- as.matrix(in_f)
> data.cl <- c(rep(1, param_A), rep(2, param_B))
> data.cl
[1] 1 1 1 2 2 2
> tcc <- new("TCC", data, data.cl)  
> tcc <- calcNormFactors(tcc,iteration=param1)
TCC::INFO: Calculating normalization factors using DEGES
TCC::INFO: (iDEGES pipeline : tmm - [ edger - tmm ] X 3 )
TCC::INFO: Done.
> tcc <- estimateDE(tcc, test.method=param_DEmethod, FDR=param_FDR)
TCC::INFO: Identifying DE genes using edger ...
TCC::INFO: Done.
> result <- getResult(tcc, sort=FALSE)
> tmp <- cbind(rownames(tcc$count), tcc$count, result)
> write.csv(tmp, out_f, append=F, quote=T, row.names=F)
 警告メッセージ: 
 write.csv(tmp, out_f, append = F, quote = T, row.names = F) で: 
   'append' への変更の試みは無視されました 
> library(ggplot2)
> library(dplyr)

 次のパッケージを付け加えます: ‘dplyr’ 

 以下のオブジェクトは ‘package:baySeq’ からマスクされています: 

     groups 

 以下のオブジェクトは ‘package:Biobase’ からマスクされています: 

     combine 

 以下のオブジェクトは ‘package:matrixStats’ からマスクされています: 

     count 

 以下のオブジェクトは ‘package:GenomicRanges’ からマスクされています: 

     intersect, setdiff, union 

 以下のオブジェクトは ‘package:GenomeInfoDb’ からマスクされています: 

     intersect 

 以下のオブジェクトは ‘package:IRanges’ からマスクされています: 

     collapse, desc, intersect, setdiff, slice, union 

 以下のオブジェクトは ‘package:S4Vectors’ からマスクされています: 

     first, intersect, rename, setdiff, setequal, union 

 以下のオブジェクトは ‘package:BiocGenerics’ からマスクされています: 

     combine, intersect, setdiff, union 

 以下のオブジェクトは ‘package:stats’ からマスクされています: 

     filter, lag 

 以下のオブジェクトは ‘package:base’ からマスクされています: 

     intersect, setdiff, setequal, union 

>
> packageVersion("ggplot2")
[1] ‘3.3.3’
> packageVersion("dplyr")
[1] ‘1.0.5’
> countdata <- read.csv("TCC_simulation_data.csv")
> head(countdata)
  rownames.tcc.count. G1_rep1 G1_rep2 G1_rep3 G2_rep1 G2_rep2 G2_rep3 gene_id  a.value   m.value
1              gene_1      50      27      32      10      13      10  gene_1 4.320660 -1.688783
2              gene_2      23      55      28       8       4       0  gene_2 3.571005 -3.115365
3              gene_3     246     256     158     106      54      74  gene_3 7.032787 -1.463880
4              gene_4      87      47      68      14      21      10  gene_4 4.989496 -2.131754
5              gene_5     270     137     235      50      46     115  gene_5 6.938669 -1.569962
6              gene_6     177     151     116      50      44      47  gene_6 6.381648 -1.621381
       p.value     q.value rank estimatedDEG
1 0.0105380207 0.078700677 1339            0
2 0.0008475635 0.008552609  991            1
3 0.0008382452 0.008492859  987            1
4 0.0002320392 0.002726665  851            1
5 0.0013314925 0.012827481 1038            1
6 0.0001367345 0.001700678  804            1
> countdata <- mutate(countdata, thcolor = ifelse(p.value >= 0.05 , 0,
+                                    ifelse(m.value >= 1, 1,
+                                           ifelse(m.value <= -1, 2,0))))
>
> ggplot() +
+   geom_point(aes(x = countdata$m.value, y = -log10(countdata$p.value), color = as.factor(countdata$thcolor)), size = 0.1) +
+   labs(x = "log2 (Fold Change)", y = "-log10 (p-value)") +
+   scale_color_manual(values = c("black", "red", "blue")) +
+   geom_hline(yintercept = -log10(0.05), size = 0.2, color = "dark green") + 
+   geom_vline(xintercept = -1, size = 0.2, color = "dark green") +
+   geom_vline(xintercept = 1, size = 0.2, color = "dark green") +
+   theme(legend.position = "none", panel.grid = element_blank())

reference

Jianqiang Sun, Tomoaki Nishiyama, Kentaro Shimizu, Koji Kadota. TCC: an R package for comparing tag count data with robust normalization strategies. BMC Bioinformatics. 2013, 14(219) doi: 10.1186/1471-2105-14-219.