Rでvolcano plotを描いてみる #RStudio

参考

ggplot2で描いてみる

サンプルデータをダウンロードする

suppressPackageStartupMessages(library(ggplot2))
suppressPackageStartupMessages(library(ggrepel))

# Download the data we will use for plotting
download.file("https://raw.githubusercontent.com/biocorecrg/CRG_RIntroduction/master/de_df_for_volcano.rds", "de_df_for_volcano.rds", method = "curl")

# The RDS format is used to save a single R object to a file, and to restore it.
# Extract that object in the current session:
tmp <- readRDS("de_df_for_volcano.rds")

# Remove rows that contain NA values
de <- tmp[complete.cases(tmp), ]

データをみてみる

head(de)

gene_symbol	log2FoldChange	pvalue
Ndufa9	-0.8764536	3.385150e-34
Rpp21	-0.9764875	1.560439e-26
Amotl1	-0.7582692	2.164811e-24
Gm44064	-0.8485025	5.125048e-20
Hcfc1r1	-0.6021286	3.799504e-17
Rgs18	-0.8171239	1.055586e-16

発現変動遺伝子をラベルする

# set thresholds
padj <- 0.05
logFC <- 0.5
de$significant <- ifelse(de$pvalue < padj & abs(de$log2FoldChange) > logFC,
                            paste0("FDR<", padj) , "Not Sig")

volcano plot関数を定義してプロットする

# Define the "volcano_plot" function
volcano_plot <- function(x){
  ggplot(x, aes(x = log2FoldChange, y = -log10(pvalue))) +
    geom_point(size = 3, aes(color = significant)) +
    scale_color_manual(values = c("#F8766D", "grey")) +
    theme(legend.position = "none") +
    geom_text_repel(
      data = subset(x, pvalue < padj & abs(x$log2FoldChange) > logFC),
      aes(label = gene_symbol), size = 4, box.padding = unit(0.15, "lines"),
      point.padding = unit(0.15, "lines")) +
    labs(x = "Log2FoldChange", y = "-log10(padj)")
}
volcano_plot(de)

出力結果をみてみる

それぞれの遺伝子をプロットし、発現変動遺伝子のみをラベルすることができた