ggplot2 (R 언어)로 volcano plot을 그리는 RNA-seq 분석

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 패키지를 사용하여 데이터를 만들었습니다.
TCC package에 대해서는, Sun 등의 논문을 참조해 주세요 (페이지의 최하부에 기재하고 있습니다).

다음 스크립트를 사용하여 시뮬레이션 데이터를 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())

참조

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 .