Volcano plot 기초부터 심화까지

그림 1. 완성 volcano plot — 오른쪽 위(코랄)=유의하게 증가, 왼쪽 위(딥그린)=유의하게 감소, 가운데 아래(회색)=유의하지 않음.

이 글에서 만들 것 — RNA-seq 차등발현(DEG) 결과를 한 장으로 보여 주는 volcano plot을 R의 ggplot2로 처음부터 그립니다. 패키지 설치 → 데이터 준비 → 색·임계선 → 유전자 라벨까지 복붙으로 따라 하면 발표·논문에 바로 쓸 그림이 나옵니다. 예제 데이터를 코드 안에서 만들기 때문에 별도 파일 없이 그대로 실행됩니다.

🔗 관련 글  ·  관련: NGS  ·  관련: 유전자 발현

1. Volcano plot이 뭔가요?

volcano plot(화산 그래프)은 수천~수만 개 유전자의 차등발현 결과를 한 번에 보여 주는 산점도입니다. 점 하나가 유전자 하나죠.

• 가로축 = log2 fold change: 두 조건 사이에서 발현이 몇 배 변했는지 (오른쪽=증가, 왼쪽=감소)
• 세로축 = −log10(adjusted p): 통계적으로 얼마나 믿을 만한지 (위로 갈수록 유의)

유의하게 변한 유전자들이 좌우 위쪽으로 솟아오른 모양이 화산 분출을 닮아 이런 이름이 붙었습니다. 앞선 NGS 글에서 읽어 낸 발현 데이터를, 유전자 발현 글에서 본 "어떤 유전자가 켜지고 꺼졌는지"로 해석하는 마지막 시각화 단계라고 보면 됩니다.

2. 사전 준비

R과 RStudio가 설치돼 있다고 가정합니다. 패키지 세 개만 준비하면 됩니다.

# 한 번만 실행 (이미 깔려 있으면 건너뛰기)
install.packages(c("ggplot2", "dplyr", "ggrepel"))

패키지	역할
ggplot2	그래프 그리기
dplyr	데이터 가공(컬럼 추가·필터)
ggrepel	유전자 라벨이 서로 겹치지 않게 밀어내기

3. 단계별 실행

STEP 1 — 데이터 준비

실제 분석에서는 DESeq2 결과를 그대로 씁니다.

# 실제 워크플로우라면 이렇게 시작합니다
# res <- as.data.frame(results(dds))   # DESeq2 결과 → 데이터프레임
# res$gene <- rownames(res)

이 글에서는 설치 없이 바로 따라 할 수 있도록 가상 데이터를 만듭니다. 모양은 실제 DESeq2 결과(log2FoldChange · pvalue · padj 컬럼)와 똑같습니다.

library(ggplot2)
library(dplyr)
library(ggrepel)

set.seed(42)                                  # 재현성 고정
n <- 2000
log2fc <- rnorm(n, mean = 0, sd = 1.2)        # 발현 변화량
pval   <- runif(n) * exp(-abs(log2fc))        # 변화가 클수록 작은 p (경향)

res <- data.frame(
  gene           = paste0("Gene", 1:n),
  log2FoldChange = log2fc,
  pvalue         = pval,
  padj           = p.adjust(pval, method = "BH")   # 다중검정 보정(BH)
)

head(res)

   gene log2FoldChange   pvalue     padj
1 Gene1        1.65574   0.0421   0.0751
2 Gene2       -0.67881   0.5123   0.6200
...

왜 padj인가? 유전자 수천 개를 한꺼번에 검정하면 우연히 "유의"해 보이는 게 쏟아집니다. 그래서 다중검정 보정을 거친 adjusted p(padj)를 기준으로 삼습니다.

STEP 2 — 결측치(NA) 제거

DESeq2의 padj에는 NA가 흔합니다(발현량이 너무 낮아 검정에서 빠진 유전자 등). 미리 걸러 줘야 그림이 깨지지 않습니다.

res <- res %>% filter(!is.na(padj))

STEP 3 — 유의성 분류

각 유전자를 증가(Up) / 감소(Down) / 유의하지 않음(NS)으로 나눕니다. 흔히 쓰는 기준은 |log2FC| > 1 그리고 padj < 0.05입니다.

fc_cut <- 1
p_cut  <- 0.05

res <- res %>%
  mutate(sig = case_when(
    log2FoldChange >  fc_cut & padj < p_cut ~ "Up",
    log2FoldChange < -fc_cut & padj < p_cut ~ "Down",
    TRUE                                    ~ "NS"
  ))

table(res$sig)     # 그룹별 개수 확인

Down   NS   Up
 158 1690  152      # (시드에 따라 값은 조금 달라집니다)

STEP 4 — 기본 산점도

먼저 색 없이 뼈대만 그려 봅니다.

ggplot(res, aes(x = log2FoldChange, y = -log10(padj))) +
  geom_point()

그림 2. 색·임계선 없는 기본 산점도 — 화산 윤곽은 보이지만 어떤 점이 의미 있는지 아직 알 수 없습니다.

STEP 5 — 색 + 임계선 추가

sig로 색을 입히고, 임계값 위치에 점선을 그립니다. 블로그 색(증가=코랄, 감소=딥그린)을 그대로 썼습니다.

cols <- c("Up" = "#E76F51", "Down" = "#2E7D6B", "NS" = "grey80")

ggplot(res, aes(log2FoldChange, -log10(padj), color = sig)) +
  geom_point(alpha = 0.7, size = 1.8) +
  scale_color_manual(values = cols) +
  geom_vline(xintercept = c(-fc_cut, fc_cut), linetype = "dashed", color = "grey50") +
  geom_hline(yintercept = -log10(p_cut),      linetype = "dashed", color = "grey50") +
  theme_minimal(base_size = 13)

그림 3. 색과 임계선을 더하면 오른쪽 위(코랄)·왼쪽 위(딥그린)의 의미 있는 유전자가 또렷해집니다.

STEP 6 — 라벨 + 마무리 (최종)

가장 유의한 유전자 10개에만 이름을 답니다. 라벨이 겹치지 않게 geom_text_repel을 씁니다. (결과 = 맨 위 그림 1)

top <- res %>% filter(sig != "NS") %>% arrange(padj) %>% head(10)

ggplot(res, aes(log2FoldChange, -log10(padj), color = sig)) +
  geom_point(alpha = 0.7, size = 1.8) +
  scale_color_manual(values = cols, name = NULL) +
  geom_vline(xintercept = c(-fc_cut, fc_cut), linetype = "dashed", color = "grey50") +
  geom_hline(yintercept = -log10(p_cut),      linetype = "dashed", color = "grey50") +
  geom_text_repel(data = top, aes(label = gene), size = 3.2,
                  max.overlaps = 20, show.legend = FALSE) +
  labs(x = expression(log[2]~"fold change"),
       y = expression(-log[10]~"adjusted p")) +
  theme_minimal(base_size = 13) +
  theme(legend.position = "top")

ggsave("volcano.png", width = 7, height = 6, dpi = 300)   # 저장

geom_text_repel이 라벨을 자동으로 밀어내 겹치지 않게 배치하고, ggsave로 고해상 PNG까지 저장됩니다. 이게 맨 위 그림 1의 완성본입니다.

4. 자주 발생하는 에러 + 해결법

증상	원인 / 해결
점이 일부 사라지고 경고 발생	padj의 NA 때문 → filter(!is.na(padj)) 먼저
could not find function "geom_text_repel"	ggrepel 미설치/미로드 → library(ggrepel)
위쪽 점이 무한대로 튐	padj == 0 → -log10(0)=Inf. 아주 작은 값으로 치환: res$padj[res$padj==0] <- 1e-300
색이 안 입혀짐	aes() 안에 color = sig를 빠뜨림
라벨이 너무 많아 지저분	head(10)처럼 상위 N개만 라벨

5. 확장 / 응용

① 패키지로 한 번에 — EnhancedVolcano
publication-ready 그림을 한 줄로 뽑고 싶다면:

# BiocManager::install("EnhancedVolcano")
library(EnhancedVolcano)
EnhancedVolcano(res, lab = res$gene,
                x = "log2FoldChange", y = "padj",
                pCutoff = 0.05, FCcutoff = 1)

② 인터랙티브로 — plotly
ggplotly()로 감싸면 마우스를 올렸을 때 유전자명이 뜨는 인터랙티브 그래프가 됩니다.

References

1. Love, M. I., Huber, W., & Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15(12), 550.
2. Wickham, H. (2016). ggplot2: Elegant Graphics for Data Analysis (2nd ed.). Springer.
3. Blighe, K., Rana, S., & Lewis, M. EnhancedVolcano: Publication-ready volcano plots. Bioconductor.
4. Slowikowski, K. ggrepel: Automatically Position Non-Overlapping Text Labels with ggplot2. CRAN.
5. Benjamini, Y., & Hochberg, Y. (1995). Controlling the false discovery rate. J. R. Stat. Soc. B, 57(1), 289–300.

Pipette & Pipeline · A bio portfolio journal — bridging research, data, and industry.