RでVolcano Plot

　Volcano Plotは発現変動遺伝子の解析などで用いられる図になります。発現変動遺伝子はその発現比を用いれば何倍になったかは分かりますが、その妥当性(p値)までは分かりません。逆にp値だけではどのくらいの比率で変わっているかの情報は得られません。それを同時にプロットしたものがVolcano Plotになります。以下参考記事

ggplot2 (R言語)でvolcano plotを描く RNA-seq解析 - Qiita

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("TCC")

library(TCC)

#シミュレーションデータを作成しcsvファイルに書き出し。
#10000の遺伝子とその発現パターンになっている。
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)
#rep関数で1をparam_Aだけ並べる
data.cl <- c(rep(1, param_A), rep(2, param_B))
data.cl
tcc <- new("TCC", data, data.cl) 
 
#多段階正規化手法を使用して正規化ファクターを計算
tcc <- calcNormFactors(tcc,iteration=param1)

#p値の計算
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)

生成されるm.valueはfold changeを表します。

#ここから図の作成準備
library(ggplot2)
library(dplyr)

countdata <- read.csv("TCC_simulation_data.csv")
head(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が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())

出力結果が以下のようになります。

実際にやる場合、変動前後のカウントデータ(csvファイル)のみもっている状態だと仮定して最初の方のコードも書きます。

library(TCC)
library(ggplot2)
library(dplyr)

count2 <- read.csv("count.csv" ,header = TRUE,row.names = 1)
head(count2)

#標本数によってここの数字を変える
param_A <- 3
param_B <- 3
param1 <- 3
param_DEmethod <- "edger"
param_FDR <- 0.05

data2 <- count2

data.cl_2 <- c(rep(1, param_A), rep(2, param_B))
data.cl_2
tcc <- new("TCC", data2, data.cl_2)  

あとは上と同じです。こんな感じのデータを用意すればＯＫです。