RNA-Seq 볼케이노 플롯: 컷오프, 코드, 함정

볼케이노 플롯은 하나의 차등 발현 검정에서 나온 모든 유전자를 한 패널에 표시합니다. x축은 log2 배수 변화이고, y축은 조정된 p-value의 −log10입니다. 중심에서 멀리 떨어져 있고 위쪽에 있는 유전자가 후보입니다.

이 플롯은 코드 열 줄이면 만들 수 있습니다. 하지만 올바르게 읽는 데는 더 많은 주의가 필요합니다. 대부분의 사람들이 복사해서 쓰는 두 개의 점선, |log2FC| ≥ 1 및 padj < 0.05는 규칙이 아니라 표시 관례일 뿐입니다. 그리고 edgeR 사용자 가이드는 이 조합으로 유전자를 고르면 "일반적으로 FDR의 조절을 파괴"한다고 말합니다(섹션 2.14).

이 글에서는 각 축이 의미하는 바, 어떤 컷오프를 그려야 하는지와 그 이유, 구름의 각 영역을 읽는 방법, 그리고 동작하는 R 및 Python 코드를 다룹니다. 그림만 필요하다면, 저희 볼케이노 플롯 생성기가 결과 테이블을 받아 레이블이 표시된 플롯을 반환합니다.

시뮬레이션 데이터(12,000개 유전자, 5%는 실제 효과 보유)로 만든 볼케이노 플롯. matplotlib로 제작했으며, 데이터를 생성한 코드는 아래 Python 섹션에 있고, 전체 플롯 스크립트는 저장소에서 이 이미지 옆에 있습니다.

두 축, 그리고 둘 다 로그를 취하는 이유

두 축 모두 로그를 취하는 이유는 같습니다. 원래 스케일은 관심 있는 부분을 숨기기 때문입니다.

원래 스케일의 배수 변화는 한쪽으로 치우칩니다. 두 배가 되는 유전자는 2.0, 절반이 되는 유전자는 0.5가 되어, 상향은 무한대까지 공간이 있는 반면 하향은 1 아래에 갇힙니다. log2는 이를 해결합니다. 두 배는 +1, 절반은 −1이 되어 플롯의 양쪽이 거울처럼 대칭이 됩니다.

p-value는 반대 문제입니다. 흥미로운 값들이 모두 0 근처에 몰려 있습니다. −log10을 취하면 그 뭉침이 펼쳐지고, 마이너스 부호 덕분에 가장 강한 증거가 위쪽에 오게 됩니다.

이 표에서 두 가지 사실을 알 수 있습니다. y = 4에 있는 점은 y = 2에 있는 점보다 padj가 100배 더 작습니다. 그리고 y축에는 상한이 없으므로 p-value가 아주 작은 유전자 하나가 플롯 전체를 늘어뜨릴 수 있습니다.

y축에는 원래 p-value가 아닌 padj를 넣으세요

y축에 원래 p-value를 사용하면 수백 개의 가짜 발견을 얻게 됩니다. 벌크 RNA-seq은 모든 유전자를 검정하므로, 실제 효과가 없는 15,000개 유전자를 검정해도 p < 0.05인 유전자가 약 750개 나옵니다.

해결책은 오탐지율(FDR) 보정입니다. Benjamini-Hochberg 방법은 p-value를 재조정하여 0.05 컷오프가 "전체 검정의 5%"가 아니라 "내가 유지하는 유전자의 약 5%가 거짓"임을 의미하게 합니다. R에서는 Benjamini와 Hochberg(1995) 이후 p.adjust에서 method = "BH"로 수행합니다.

두 주요 도구는 이미 보정된 열을 제공합니다.

• DESeq2 (비네트): log2FoldChange, pvalue, padj. padj 열은 기본적으로 BH이며, results()가 pAdjustMethod = "BH"로 실행되기 때문입니다.
• edgeR (사용자 가이드): topTags()는 logFC, logCPM, 검정 통계량 열(F 또는 LR), PValue, FDR을 반환합니다. 원하는 열은 FDR입니다.

여기서 사람들을 곤란하게 하는 DESeq2 세부 사항이 두 가지 있습니다. 첫째, results()는 기본적으로 alpha = 0.1을 사용하며, 독립 필터링 단계도 그 값에 맞춰져 있습니다. 점선이 0.05에 있다면 results(dds, alpha = 0.05)를 호출하여 필터링이 그린 선과 일치하게 하세요.

둘째, 일부 유전자는 padj = NA로 반환됩니다. 이는 비네트의 "일부 p-value가 NA로 설정되는 이유는 무엇인가?" 섹션에 따르면 독립 필터링과 이상치 제거가 제 역할을 하고 있다는 뜻입니다. 플롯을 그리기 전에 해당 행을 제거하세요. NA를 1이나 0으로 바꾸지 마세요.

컷오프 선은 어디에 그릴까

표준 컷오프 조합은 없으며, 흔히 쓰이는 조합만 있을 뿐입니다. 인기 있는 숫자가 실제로 어디서 나왔는지는 다음과 같습니다.

습관이 아니라 신호의 양에 따라 선택하세요. 6,000개 유전자가 통과하면 목록은 쓸모없으므로 더 높은 배수 변화 기준이 도움이 됩니다. 9개 유전자만 통과하면 배수 변화 선을 빼고 padj 순으로 순위를 매기세요.

edgeR 가이드에는 x축 컷오프에 대한 가장 좋은 한 줄 규칙이 있습니다. 이를 "유전자에 확실히 관심이 없는 배수 변화 미만"으로 읽어야지, 유전자가 흥미로워지는 변화 이상으로 읽으면 안 됩니다.

알아둘 만한 컷오프 함정

p로 필터링한 다음 배수 변화로 필터링하는 것은 유효한 검정이 아닙니다. edgeR 가이드는 이에 대해 솔직합니다. 그러한 조합은 "모두 임시방편이며, 배수 변화 임계값에 대한 차등 발현 검정에 의미 있는 p-value를 제공하지 않습니다". 또한 "발현량이 낮지만 변동성이 큰 유전자를 선호"하며 "일반적으로 FDR의 조절을 파괴"합니다.

배수 변화 기준이 주장에 중요하다면, 사후 필터링 대신 그 기준에 대해 검정하세요.

# DESeq2: test whether |log2FC| is above 1, not just above 0
res <- results(dds, lfcThreshold = 1, altHypothesis = "greaterAbs", alpha = 0.05)

# edgeR: same idea, via the TREAT method
fit <- glmQLFit(y, design)
tr  <- glmTreat(fit, coef = 2, lfc = 1)

이제 padj는 크기 컷오프를 이미 반영하며, 세로선은 실행한 검정을 설명합니다. 두 옵션 모두 문서화되어 있습니다. DESeq2의 임계값 이상의 log2 배수 변화 검정과 glmTreat의 기반이 되는 TREAT 논문을 참조하세요.

플롯을 영역별로 읽는 방법

레이블을 읽기 전에 구름부터 읽으세요. 형태는 어떤 단일 유전자보다 실험에 대해 더 많은 것을 말해줍니다.

위쪽 중간은 사람들이 잘못 읽는 영역입니다. 그 유전자들은 대개 표준 오차가 작은 고발현 유전자라서 15% 변화도 p 컷오프를 통과합니다. 변화는 실제이지만, 15% 변화가 중요한지는 통계 문제가 아니라 생물학 문제입니다. 한 가지 경고 신호는 하우스키핑 유전자나 리보솜 유전자가 그 위에 있다면 생물학이 아니라 정규화 또는 라이브러리 구성 문제를 의심해야 한다는 것입니다.

아래쪽 가장자리는 그 반대입니다. 한 그룹에서 4개, 다른 그룹에서 30개의 카운트를 가진 유전자는 엄청난 log2 배수 변화를 보이지만 증거는 거의 없습니다. 이런 점들은 x축을 넓게 만들고 이야기를 흥미진진하게 만들지만, 재현되는 경우는 드뭅니다. 수축된 배수 변화가 표시를 해결합니다. lfcShrink(dds, coef = 2, type = "apeglm")는 측정이 부정확한 유전자를 0 쪽으로 끌어당기고, 잘 측정된 유전자는 그대로 둡니다.

두 가지 형태 확인으로 읽기를 마무리합니다.

• 대칭성. 900개 상향, 40개 하향 같은 한쪽으로 치우친 구름은 실제 활성화일 수 있습니다. 또한 정규화가 흡수하지 못한 구성 변화에서 비롯될 수도 있습니다. 문장을 쓰기 전에 MA 플롯과 크기 인자를 확인하세요.
• V의 너비. x = 0 주변의 간격과 양쪽 측면에서 점들이 올라가는 것은 예상된 현상입니다. 변화가 클수록 p 컷오프를 더 쉽게 통과하기 때문입니다. 대신 x = 0에서 수직으로 솟은 뾰족한 형태는 보통 원래 p-value를 플로팅했음을 의미합니다.

DESeq2 결과로 R에서 플롯 그리기

이 코드는 dds라는 DESeqDataSet에서 실행되며 ggplot2, ggrepel, dplyr이 필요합니다.

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

lfc_cut  <- 1
padj_cut <- 0.05

# keep alpha equal to the line you draw, or filtering optimises for 0.1
res <- results(dds, alpha = padj_cut)

df <- as.data.frame(res)
df$gene <- rownames(df)
df <- df[!is.na(df$padj), ]
df$padj <- pmax(df$padj, .Machine$double.xmin)  # padj can underflow to 0

df$class <- "Not significant"
df$class[df$padj < padj_cut & df$log2FoldChange >=  lfc_cut] <- "Up"
df$class[df$padj < padj_cut & df$log2FoldChange <= -lfc_cut] <- "Down"

top <- df %>%
  filter(class != "Not significant") %>%
  arrange(padj) %>%
  slice_head(n = 10)

ggplot(df, aes(x = log2FoldChange, y = -log10(padj), colour = class)) +
  geom_point(size = 1, alpha = 0.7) +
  geom_hline(yintercept = -log10(padj_cut), linetype = "dashed") +
  geom_vline(xintercept = c(-lfc_cut, lfc_cut), linetype = "dashed") +
  geom_text_repel(
    data = top, aes(label = gene),
    colour = "black", size = 3, max.overlaps = Inf, show.legend = FALSE
  ) +
  scale_colour_manual(values = c(
    "Up" = "#c23b3b", "Down" = "#2f6fb2", "Not significant" = "#b8bec9"
  )) +
  coord_cartesian(xlim = c(-5, 5)) +
  labs(x = "log2 fold change", y = "-log10 padj", colour = NULL) +
  theme_classic()

edgeR에서 가져온다면 두 열의 이름을 바꾸면 나머지 스크립트는 그대로 작동합니다.

df <- topTags(qlf, n = Inf)$table
df$log2FoldChange <- df$logFC
df$padj <- df$FDR

연결선과 음영 처리된 컷오프 상자가 있는 완전한 스타일 버전을 원한다면 EnhancedVolcano 비네트를 읽어볼 만합니다.

pandas DataFrame으로 Python에서 플롯 그리기

DESeq2 테이블을 CSV로 내보낸 후 이 코드를 실행하세요. pandas, numpy, matplotlib만 사용합니다.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

# first column is the gene name; needs DESeq2's log2FoldChange and padj
df = pd.read_csv("deseq2_results.csv", index_col=0)
df = df.dropna(subset=["padj", "log2FoldChange"]).copy()

lfc_cut, padj_cut = 1.0, 0.05

# padj can underflow to 0, which makes -log10 infinite; clip it first
df["y"] = -np.log10(df["padj"].clip(lower=np.finfo(float).tiny))
df["cls"] = "not significant"
df.loc[(df["padj"] < padj_cut) & (df["log2FoldChange"] >= lfc_cut), "cls"] = "up"
df.loc[(df["padj"] < padj_cut) & (df["log2FoldChange"] <= -lfc_cut), "cls"] = "down"

colours = {"not significant": "#b8bec9", "down": "#2f6fb2", "up": "#c23b3b"}
fig, ax = plt.subplots(figsize=(7, 5), dpi=200)
for name, grp in df.groupby("cls"):
    ax.scatter(grp["log2FoldChange"], grp["y"], s=6, c=colours[name],
               label=name, alpha=0.7, linewidths=0)

ax.axhline(-np.log10(padj_cut), ls="--", lw=0.9, c="#444444")
for v in (-lfc_cut, lfc_cut):
    ax.axvline(v, ls="--", lw=0.9, c="#444444")

top = df[df["cls"] != "not significant"].nsmallest(10, "padj")
for gene, row in top.iterrows():
    ax.annotate(gene, (row["log2FoldChange"], row["y"]), xytext=(4, 4),
                textcoords="offset points", fontsize=7)

ax.set_xlabel("log2 fold change")
ax.set_ylabel("-log10 padj")
ax.set_xlim(-5, 5)
ax.legend(frameon=False, markerscale=2)
fig.tight_layout()
fig.savefig("volcano.png")

Python에서 edgeR 출력을 사용하려면 먼저 logFC를 log2FoldChange로, FDR을 padj로 이름을 바꾸세요.

위 그림의 배후 코드

예시 그림은 실제 데이터가 아닌 시뮬레이션 데이터를 사용합니다. 모델은 의도적으로 단순합니다. 기본 평균이 네 자릿수에 걸쳐 분포하고, 5%의 유전자가 실제 효과를 가지며, 기본 평균이 낮아질수록 표준 오차가 커집니다. 마지막 부분이 먼 가장자리에 증거가 약한 점들을 만듭니다.

import numpy as np
from math import erfc

rng = np.random.default_rng(11)
n = 12000

base_mean = 10 ** rng.uniform(0, 4.2, n)
true_lfc = np.where(rng.random(n) < 0.05, rng.normal(0, 1.1, n), 0.0)
se = np.sqrt(0.035 + 9.0 / base_mean)        # low-count genes get a larger standard error
lfc = true_lfc + rng.normal(0, se)
pval = np.array([erfc(abs(s) / np.sqrt(2)) for s in lfc / se])

# Benjamini-Hochberg, same as R's p.adjust(method = "BH")
order = np.argsort(pval)
ranked = pval[order] * n / np.arange(1, n + 1)
padj = np.empty(n)
padj[order] = np.clip(np.minimum.accumulate(ranked[::-1])[::-1], 0, 1)
padj[base_mean < 1.5] = np.nan                # stand-in for independent filtering

볼케이노 플롯을 오해하게 만드는 다섯 가지 실수

• y축에 원래 p-value 사용. 수백 개의 유전자가 우연히 선을 넘습니다. padj 또는 FDR을 사용하고, 축 레이블에 어느 것인지 명시하세요.
• 배수 변화 선을 생물학적 의미로 취급. 전사 인자의 2배 변화와 구조 단백질의 2배 변화는 비교할 수 없는 주장입니다. 그 선은 필터이지 증거가 아닙니다.
• 이상치가 축을 늘어뜨리도록 방치. log2FC = 12인 유전자 하나가 다른 모든 것을 납작하게 만듭니다. coord_cartesian() 또는 set_xlim()으로 보기를 자르고, 잘린 점은 버리지 말고 가장자리에 삼각형으로 표시하세요.
• 배수 변화만으로 색칠. x = 1을 지난 모든 유전자를 빨간색으로 표시하면 노이즈가 많은 저발현 유전자가 발견으로 표시됩니다. 색칠에는 두 조건이 모두 필요합니다. 배수 변화 선을 넘었는지 그리고 padj 선을 넘었는지.
• 점선을 전혀 그리지 않음. 점선이 없으면 독자가 플롯에서 컷오프를 알 수 없습니다. 두 선을 모두 그리고, 숫자를 캡션이나 범례에 넣으세요.

자주 묻는 질문

볼케이노 플롯은 MA 플롯과 같은가요?

아니요, 두 플롯은 서로 다른 질문에 답합니다. MA 플롯은 x축에 평균 발현량, y축에 log2 배수 변화를 두어 큰 변화가 저발현 유전자에서 나오는지 보여줍니다. 볼케이노 플롯은 발현량을 빼고 대신 증거를 보여줍니다. 작업 중에는 둘 다 만들고, 출판할 때는 볼케이노 플롯을 사용하세요.

padj가 정확히 0이면 어떻게 하나요?

그것은 p-value가 0이라는 뜻이 아니라 부동소수점 언더플로입니다. -log10(0)은 무한대가 되어 점이 사라지거나 축이 깨집니다. 위 Python 코드처럼 로그를 취하기 전에 padj를 가장 작은 양의 double로 클리핑하고, 텍스트에서는 해당 유전자를 padj < 1e-300으로 보고하세요.

유전자 레이블은 몇 개나 붙여야 하나요?

10~20개는 읽을 만하고, 그 이상은 지저분합니다. 배수 변화 선도 통과하는 유전자 중에서 padj 순으로 순위를 매긴 다음 고정된 개수만큼 레이블을 붙이세요. R에서는 ggrepel을, Python에서는 위치를 조정한 annotate 호출을 사용하여 텍스트가 점 위에 놓이지 않게 하세요. 초록에서 언급하는 유전자는 상위 10개에 없더라도 레이블을 붙이세요.

볼케이노 플롯을 논문에 바로 넣을 수 있나요?

네, 벡터 아트로 내보내거나 저널이 요구하는 해상도로 내보내고, 기본 분석이 재현 가능하다면 가능합니다. 출판사는 데이터 그림과 일러스트레이션을 다르게 취급하며, 제출 전에 규칙이 중요합니다. AI 생성 그림에 대한 저널 정책에 대한 저희 가이드를 참조하세요. 논문의 요약 패널에는 그래픽 초록 생성기가 개략적인 부분을 담당하고, 볼케이노 플롯은 여러분 자신의 데이터로 만든 그림으로 남습니다.