RNA-Seqボルケーノプロット：カットオフ、コード、そして落とし穴

ボルケーノプロットは、1回の差次的発現解析のすべての遺伝子を1つのパネルに描画します。x軸はlog2 fold changeです。y軸は調整後p-valueの−log10です。中心から遠く、かつ上方にある遺伝子が候補遺伝子です。

このプロットは10行のコードで作成できます。正しく読み取るにはもう少し注意が必要です。多くの人がコピーする2本の破線、|log2FC| ≥ 1 と padj < 0.05 は、表示上の習慣であってルールではありません。そしてedgeRユーザーガイドは、この組み合わせで遺伝子を選ぶと「一般にFDRの制御を破壊する」と述べています（セクション2.14）。

この記事では、各軸の意味、どのカットオフを引くべきかとその理由、雲状の各領域の読み方、そして動作するRコードとPythonコードを解説します。図だけが必要な場合は、ボルケーノプロット作成ツールが結果テーブルを受け取り、ラベル付きのプロットを返します。

シミュレーションデータ（12,000遺伝子、5%が真の効果を持つ）から作成したボルケーノプロット。matplotlibで作成。データを生成したコードは後述のPythonセクションにあり、完全なプロットスクリプトはこの画像と同じリポジトリ内にあります。

2つの軸、そしてなぜ両方とも対数なのか

両方の軸が対数なのは同じ理由です。生のスケールでは、注目すべき部分が隠れてしまうからです。

生のスケールでのfold changeは偏っています。2倍になる遺伝子は2.0、半分になる遺伝子は0.5となるため、上方は無限大まで広がる一方、下方は1未満に押し込められます。log2はこれを修正します。2倍は+1、半分は−1となり、プロットの両側は鏡像になります。

p-valueは逆の問題があります。p-valueはゼロ付近に密集し、そこに興味深い値がすべてあります。−log10を取るとその密集が広がり、マイナス符号によって最も強いエビデンスが上部に来ます。

この表から2つの事実が導かれます。y = 4 の点は、y = 2 の点よりpadjが100分の1です。そしてy軸に上限はないため、1つの遺伝子の極小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ではp.adjustでmethod = "BH"を指定して実行します。Benjamini and Hochberg (1995)に基づきます。

主要な2つのツールは、すでに補正済みの列を提供しています。

• DESeq2（vignette）：log2FoldChange、pvalue、padj。results()はpAdjustMethod = "BH"で実行されるため、padj列はデフォルトでBHです。
• edgeR（ユーザーガイド）：topTags()はlogFC、logCPM、検定統計量の列（FまたはLR）、PValue、FDRを返します。必要な列はFDRです。

ここでDESeq2に関する2つの注意点があります。まず、results()はデフォルトでalpha = 0.1を使用し、独立フィルタリングのステップもその値に合わせて調整されます。破線を0.05に引く場合は、フィルタリングを描いた線に合わせるため、results(dds, alpha = 0.05)を呼び出してください。

次に、一部の遺伝子はpadj = NAで返されます。これは独立フィルタリングと外れ値除去が機能しているためで、vignetteの「Why are some p values set to NA?」セクションに説明があります。プロット前にこれらの行を除外してください。NAを1や0に変換してはいけません。

カットオフ線をどこに引くか

標準的なカットオフの組み合わせはなく、よく使われるものがあるだけです。一般的な数値の由来は以下のとおりです。

習慣ではなく、シグナルの量で選んでください。6,000遺伝子が通過するなら、そのリストは役に立たないので、より高いfold changeの基準が役立ちます。9遺伝子しか通過しないなら、fold changeの線を外し、padjで順位付けしてください。

edgeRガイドには、x軸カットオフに関する最良の一言ルールがあります。「その遺伝子に確実に興味がないfold changeの値」と読み、遺伝子が興味深くなる閾値とは読まないでください。

知っておくべきカットオフの罠

p値でフィルタリングしてからfold changeでフィルタリングするのは、有効な検定ではありません。edgeRガイドは率直に述べています。そのような組み合わせは「アドホックであり、fold change閾値に対する差次的発現の検定として意味のあるp値を与えない」と。また「低発現だが変動の大きい遺伝子を優先し」「一般にFDRの制御を破壊する」と警告しています。

fold changeの基準が主張に関わるなら、事後的にフィルタリングするのではなく、その基準に対して検定してください。

# DESeq2: |log2FC| が 1 より大きいかを検定（0 より大きいかではない）
res <- results(dds, lfcThreshold = 1, altHypothesis = "greaterAbs", alpha = 0.05)

# edgeR: TREAT法による同様の考え方
fit <- glmQLFit(y, design)
tr  <- glmTreat(fit, coef = 2, lfc = 1)

これでpadjは大きさのカットオフをすでに考慮し、垂直線は実行した検定を表します。どちらのオプションも文書化されています。DESeq2のlog2 fold changeの閾値に関する検定と、glmTreatの背後にあるTREAT論文を参照してください。

プロットの読み方：領域ごとに

ラベルを読む前に雲全体を読んでください。その形状は、個々の遺伝子よりも実験について多くを語ります。

上部中央は人々が誤読する領域です。これらの遺伝子は通常、高カウントの遺伝子で、標準誤差が小さいため、15%の変化でもp値のカットオフを通過します。変化は本物ですが、15%の変化が重要かどうかは統計学の問題ではなく生物学の問題です。1つの警告サイン：ハウスキーピング遺伝子やリボソーム遺伝子がそこにある場合、生物学ではなく、正規化またはライブラリ構成の問題を疑ってください。

下部の両端はその鏡像です。一方の群で4カウント、もう一方で30カウントの遺伝子は、巨大なlog2 fold changeとほとんどエビデンスのない結果を与えます。これらの点はx軸を広げ、ストーリーを刺激的にしますが、再現性はほとんどありません。縮小（shrinkage）されたfold changeが表示を修正します。lfcShrink(dds, coef = 2, type = "apeglm")は、測定が不十分な遺伝子をゼロに向かって引き寄せ、良く測定された遺伝子はそのままにします。

2つの形状チェックで読み取りを締めくくります。

• 対称性。 偏った雲、例えば900が上方制御で40が下方制御、は本当の活性化である可能性があります。また、正規化が吸収しきれなかった構成シフトに由来する可能性もあります。その文を書く前にMAプロットとサイズファクターを確認してください。
• Vの幅。 x = 0 の周りのギャップと、両側で点が上昇していくのは予想通りです。大きな変化ほどp値のカットオフを通過しやすいからです。代わりにx = 0 にしっかりした垂直スパイクがある場合は、通常、生のp値をプロットしたことを意味します。

DESeq2の結果からRでプロットする

これはDESeqDataSet（ddsという名前）で動作し、ggplot2、ggrepel、dplyrが必要です。

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

lfc_cut  <- 1
padj_cut <- 0.05

# alpha を描く線と同じにしないと、フィルタリングが 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 は 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から来た場合は、2つの列名を変更すれば、残りのスクリプトはそのまま動作します。

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

コネクタと網掛けのカットオフボックス付きの完全にスタイルされたバージョンについては、EnhancedVolcano vignetteを読む価値があります。

pandas DataFrameからPythonでプロットする

DESeq2のテーブルをCSVにエクスポートしてから、これを実行します。pandas、numpy、matplotlibのみを使用します。

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

# 最初の列は遺伝子名。DESeq2 の log2FoldChange と 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 は 0 にアンダーフローし、-log10 が無限大になるため、先にクリップする
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に名前を変更してください。

上の図の背後にあるコード

サンプルの図は実データではなくシミュレーションデータを使用しています。モデルは意図的に単純です。ベース平均は4桁にわたって分布し、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)        # 低カウント遺伝子は標準誤差が大きくなる
lfc = true_lfc + rng.normal(0, se)
pval = np.array([erfc(abs(s) / np.sqrt(2)) for s in lfc / se])

# Benjamini-Hochberg。R の 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                # 独立フィルタリングの代用

ボルケーノプロットを誤解を招くものにする5つの間違い

• y軸に生のp-valueを使う。 何百もの遺伝子が偶然に線を越えます。padjまたはFDRを使用し、軸ラベルにどちらかを明記してください。
• fold changeの線を生物学的意味として扱う。 転写因子の2倍の変化と構造タンパク質の2倍の変化は、比較可能な主張ではありません。その線はフィルターであって、エビデンスではありません。
• 外れ値に軸を引き伸ばさせる。 log2FC = 12 の1つの遺伝子が他のすべてを平らにします。coord_cartesian()またはset_xlim()でビューをクリップし、クリップされた点は端に三角形でマークして、削除しないでください。
• fold changeだけで色付けする。 x = 1 を超えたすべての遺伝子を赤にすると、ノイズの多い低カウント遺伝子がヒットとしてマークされます。色付けには両方の条件が必要です。fold changeの線を超えていることかつpadjの線を超えていること。
• 破線をまったく引かない。 破線がないと、読者はカットオフをプロットから読み取れません。両方を引き、数値をキャプションまたは凡例に記載してください。

FAQ

ボルケーノプロットはMAプロットと同じですか？

いいえ、異なり、答える質問も異なります。MAプロットはx軸に平均発現量、y軸にlog2 fold changeを置くため、大きな変化が低カウント遺伝子に由来するかどうかがわかります。ボルケーノプロットは発現レベルを落とし、代わりにエビデンスを示します。作業中は両方を作成し、論文にはボルケーノプロットを掲載してください。

padjが正確に0の場合はどうすればよいですか？

それは浮動小数点のアンダーフローであり、p-valueがゼロという意味ではありません。-log10(0)は無限大になり、点が消えるか軸が壊れます。上記のPythonコードのように、対数を取る前にpadjを最小の正のdouble値にクリップし、テキストではその遺伝子をpadj < 1e-300 と報告してください。

何個の遺伝子にラベルを付けるべきですか？

10〜20個が読みやすく、それ以上は散らかります。fold changeの線も通過する遺伝子の中でpadjで順位付けし、固定数をラベル付けしてください。Rではggrepel、Pythonではオフセット付きのannotate呼び出しを使用して、テキストが点の上に乗らないようにしてください。アブストラクトで言及する遺伝子は、トップ10に入っていなくてもラベルを付けてください。

ボルケーノプロットをそのまま論文に載せられますか？

はい、ベクターアートまたはジャーナルが要求する解像度でエクスポートされ、基になる解析が再現可能であれば可能です。出版社はデータ図をイラストとは異なる扱いをし、そのルールは投稿前に重要です。AI生成図に関するジャーナルポリシーのガイドを参照してください。論文のサマリーパネルには、グラフィカルアブストラクト作成ツールが模式図部分をカバーし、ボルケーノプロットは自分の数値で描くプロットとして残ります。