【RNA-seq】発現変動解析の結果をグラフで可視化する方法（MA/火山/ヒートマップ）【R】

今回は、発現変動解析（DEG解析）の結果を視覚的に理解するための「可視化手法」 に焦点を当てて解説します。RNA-seq解析で得られる膨大な結果をわかりやすく表現するには、MAプロット、火山プロット、ヒートマップといったグラフの活用が不可欠です。この記事では、DESeq2による解析結果をもとに、RとBioconductorのパッケージを用いた実践的な可視化手法を、コードと解説付きで紹介します。

「結果をどう見せるか？」で、発見の質が大きく変わります。これを機に、可視化にも強いRNA-seq解析者を目指しましょう。

動作検証済み環境
macOS Monterey (12.4), Quad-Core Intel Core i7, Memory 32GB

発現変動解析結果の可視化について

RNA-seq実験から得られた発現変動解析の結果を効果的に可視化することは、データの理解と生物学的解釈において非常に重要です。適切な可視化によって、膨大な遺伝子発現データから有益なパターンや重要な遺伝子を迅速に識別することができます。この章では、RNA-seq解析における主要な可視化手法について説明します：

• MAプロット：平均発現量と発現変化の関係を視覚化し、発現変動パターンの全体像を把握
• 火山プロット：発現変化量（log2FoldChange）と統計的有意性（p値）を同時に表示し、生物学的・統計的に重要な遺伝子を特定
• ヒートマップ：複数サンプル間の発現パターンを色の強度で表現し、遺伝子クラスタリングやサンプルの類似性を視覚化

これらの可視化手法を通じて、実験条件間の発現変動の特徴を理解し、さらなる機能解析のための候補遺伝子を効率的に選別することができます。ここではR言語とBioconductorパッケージを活用した実践的な可視化の手法と解釈について詳しく説明します。

解析準備

GEOからカウントデータを取得

今回の解析に使うカウントデータを取得してください。

# GEOqueryパッケージを読み込む
library(GEOquery)
# GEOからカウントデータを取得
gse <- getGEO("GSE141413", GSEMatrix = FALSE)
# カウントデータを含む補足ファイルをダウンロード
supp_dir <- "data/geo_suppl"
dir.create(supp_dir, recursive = TRUE, showWarnings = FALSE)
count_file <- getGEOSuppFiles("GSE141413", baseDir = supp_dir, makeDirectory = FALSE)
# ダウンロードしたカウントデータを読み込む
counts_path <- file.path(supp_dir, "GSE141413_gene_count.txt.gz")
counts_data <- read.delim(counts_path, header = TRUE)

カウントデータの前処理

取得したカウントデータを解析に適した形式に整理します。このデータセットには発現量データと遺伝子アノテーション情報が含まれています。

# ステップ1: カウントデータの確認
dim(counts_data)
head(counts_data[, 1:8])
# ステップ2: データの構造確認
colnames(counts_data)
str(counts_data)
# ステップ3: 発現データとアノテーション情報の分離
count_cols <- c("Control", "NC_LM", "st2_LM", "IFNG", "st2_IFNG")
annotation_cols <- c("gene_name", "gene_chr", "gene_start", "gene_end", "gene_strand", "gene_length", "gene_biotype", "gene_description", "tf_family")
# サブセットの作成
counts <- counts_data[, count_cols]
gene_annotation <- counts_data[, annotation_cols]
# ステップ4: 行名の設定
rownames(counts) <- rownames(counts_data)
rownames(gene_annotation) <- rownames(counts_data)
# ステップ5: サンプル情報の作成
sample_info <- data.frame( sample_id = colnames(counts), condition = c("Control", "NC_LM", "st2_LM", "IFNG", "st2_IFNG"), group = c("Control", rep("LM", 2), rep("IFNg", 2)), treatment = c("None", "NC", "st2", "None", "st2"), row.names = colnames(counts)
)
# ステップ6: フィルタリングとデータ保存
# 空でない遺伝子のみを保持
counts <- counts[rowSums(counts) > 0, ]
gene_annotation <- gene_annotation[rownames(counts), ]
# 元の条件が複雑すぎてエラーになるため、シンプルな条件に変換
simplified_condition <- ifelse(grepl("Control", sample_info$condition), "Control", "Treatment")
sample_info$simple_condition <- factor(simplified_condition, levels = c("Control", "Treatment"))

差次的発現解析の実行

library(DESeq2)
# DESeqDataSetオブジェクトの作成（シンプルな条件を使用）
dds <- DESeqDataSetFromMatrix( countData = counts, colData = sample_info, design = ~ simple_condition # シンプル化した条件を使用
)
# 低発現遺伝子のフィルタリング
keep <- rowSums(counts(dds) >= 10) >= 2 # 少なくとも2サンプルで10リード以上
dds <- dds[keep, ]
# 差次的発現解析の実行
dds <- DESeq(dds)
# 結果の取得（TreatmentとControlの比較）
res <- results(dds, contrast = c("simple_condition", "Treatment", "Control"))

これで準備完了です。

発現変動解析結果の可視化のやり方

MA（Mean-Average）プロット

MAプロットは発現変動解析の結果を視覚化する基本的かつ重要なツールです。このプロットでは、x軸に平均発現量（Mean）、y軸にfold change（Average）をプロットし、すべての遺伝子の発現変動パターンを一度に把握することができます。

# DESeq2のMAプロット(ggplot2）
res_df <- as.data.frame(res)
res_df$significant <- ifelse(res_df$padj < 0.05 & !is.na(res_df$padj), "Yes", "No")
ggplot(res_df, aes(x = log10(baseMean), y = log2FoldChange, color = significant)) + geom_point(alpha = 0.5, size = 0.5) + scale_color_manual(values = c("No" = "gray", "Yes" = "red")) + geom_hline(yintercept = 0, linetype = "dashed") + theme_bw() + ggtitle("MA Plot") + xlab("log10(Mean Expression)") + ylab("log2FoldChange")

このMAプロット作成のコードでは、以下の重要なステップを実行しています：

1. データの準備：
   • DESeq2の結果オブジェクト（res）をデータフレームに変換して操作しやすくします。
   • significantという新しい列を作成し、調整p値（padj）が0.05未満の遺伝子を「有意」としてマークします。
   • ifelse関数を使用して、条件に基づいて「Yes」または「No」の値を割り当てます。
   • !is.na(res_df$padj)条件は、欠損値（NA）を持つ行を除外するために含まれています。
2. ggplot2によるプロット作成：
   • ggplot関数でプロットの基本構造を定義し、aes（aesthetics）でデータマッピングを指定します。
   • x = log10(baseMean)：x軸には平均発現量の対数値を指定します（可視化のため対数スケールを使用）。
   • y = log2FoldChange：y軸にはlog2スケールの発現変化量を指定します。
   • color = significant：点の色を統計的有意性に基づいて変化させます。
3. プロットのカスタマイズ：
   • geom_point：各遺伝子を点として表示します。alpha = 0.5で透明度を設定して密集した点も見やすくし、size = 0.5で点のサイズを調整します。
   • scale_color_manual：有意な遺伝子を赤色、非有意な遺伝子を灰色で表示するように色を手動設定します。
   • geom_hline(yintercept = 0)：発現変化なし（log2FoldChange = 0）を示す水平線を追加します。
   • theme_bw()：シンプルな白黒テーマを適用します。
   • ggtitle, xlab, ylab：プロットのタイトルと軸ラベルを設定します。

MAプロットを使用することで、以下の重要な特徴を視覚的に評価できます：

• 発現量に依存したバイアスの有無（点の分布が発現量に依存して偏っているか）

• 上方制御と下方制御の遺伝子数のバランス（点の分布が0を中心に対称か）

• 高発現遺伝子と低発現遺伝子における発現変動パターンの違い

  

火山プロット

火山プロット（Volcano Plot）は発現変動の大きさ（log2FoldChange）とその統計的有意性（p値）を同時に視覚化する強力なツールです。名前は、中央が窪んでその両側が山のように隆起する形状が火山に似ていることに由来します。

# 基本的な火山プロット（ggplot2）
ggplot(res_df, aes(x = log2FoldChange, y = -log10(pvalue), color = significant)) + geom_point(alpha = 0.5, size = 0.5) + scale_color_manual(values = c("No" = "gray", "Yes" = "red")) + geom_vline(xintercept = c(-1, 1), linetype = "dashed") + geom_hline(yintercept = -log10(0.05), linetype = "dashed") + theme_bw() + ggtitle("Volcano Plot") + xlab("log2FoldChange") + ylab("-log10(p-value)")
# 遺伝子シンボル情報を結果に追加
# EnsemblIDをカラムとして追加
res_df$ensembl_id <- rownames(res_df)
# 遺伝子シンボルを取得（前章で定義した関数を使用）
res_df$gene_symbol <- ensembl_to_symbol(res_df$ensembl_id)
# EnhancedVolcanoパッケージを使用した高度な火山プロット
library(EnhancedVolcano)
EnhancedVolcano(res_df, lab = res_df$gene_symbol, # 遺伝子シンボルをラベルとして使用 x = "log2FoldChange", y = "pvalue", pCutoff = 0.05, FCcutoff = 1, pointSize = 0.8, # ポイントサイズを小さく labSize = 2.5, # ラベルサイズを小さく title = "Volcano Plot", subtitle = "Treatment vs Control", axisLabSize = 10, # 軸ラベルサイズ titleLabSize = 12, # タイトルサイズ subtitleLabSize = 10, # サブタイトルサイズ captionLabSize = 8, # キャプションサイズ legendLabSize = 10, # 凡例ラベルサイズ legendIconSize = 3.2) # 凡例アイコンサイズ

火山プロットの作成には、まず基本的なggplot2によるアプローチと、より高度な機能を持つEnhancedVolcanoパッケージによるアプローチの2つが示されています。

基本的な火山プロット（ggplot2）

1. データマッピング：
   • x = log2FoldChange：x軸には発現変化量（log2スケール）を指定。負の値は下方制御、正の値は上方制御を示します。
   • y = -log10(pvalue)：y軸にはp値の負の対数を指定。小さなp値（高い統計的有意性）は大きな-log10(p)値として表示され、プロットの上部に位置します。
   • color = significant：有意な遺伝子（padj < 0.05）を赤色でハイライトします。
2. 閾値線の追加：
   • geom_vline(xintercept = c(-1, 1))：生物学的に意味のある発現変化（2倍以上の変化、log2スケールでは1）を示す垂直線を追加します。
   • geom_hline(yintercept = -log10(0.05))：統計的有意性の閾値（p = 0.05）を示す水平線を追加します。
3. プロット設定：
   • 点の透明度・サイズ設定、テーマ選択、タイトルと軸ラベルの設定はMAプロットと同様です。

高度な火山プロット（EnhancedVolcano）

EnhancedVolcanoパッケージは、火山プロットを効果的に作成するための専用機能を提供します。

1. 遺伝子シンボルの準備：
   • Ensembl IDをrownamesから取得し、ensembl_to_symbol関数を使って遺伝子シンボルに変換します。
   • この遺伝子シンボルをプロット上の各点のラベルとして使用します。
2. EnhancedVolcano関数の主要パラメータ：
   • lab：点のラベル（遺伝子シンボル）を指定します。
   • x, y：x軸とy軸のデータ列を指定します。
   • pCutoff：p値の閾値（0.05）を設定します。
   • FCcutoff：Fold Changeの閾値（1、つまり2倍の変化）を設定します。
   • 各種サイズパラメータ：点、ラベル、軸、タイトルなど様々な要素のサイズを調整します。

EnhancedVolcanoは自動的に次の4つのカテゴリに遺伝子を分類して色分けします：

• 非有意かつFold Change閾値未満の遺伝子（通常はグレー）
• 有意だがFold Change閾値未満の遺伝子（通常は青）
• 有意でないがFold Change閾値以上の遺伝子（通常はピンク）
• 有意かつFold Change閾値以上の遺伝子（通常は赤）

また、統計的に最も有意な遺伝子には自動的にラベルが表示され、視覚的に重要な遺伝子を素早く特定することができます。

火山プロットはfold changeの大きさとp値の両方を考慮して重要な遺伝子を特定するのに非常に役立ちます。プロットの右上（上方制御で有意）と左上（下方制御で有意）に位置する遺伝子が、一般的に最も注目に値する候補となります。

ヒートマップ

ヒートマップは複数のサンプル間での遺伝子発現パターンを視覚化するのに非常に効果的なツールです。色の強度によって発現レベルを表現し、行（遺伝子）と列（サンプル）の両方をクラスタリングすることで、類似した発現パターンを持つ遺伝子グループやサンプルグループを識別することができます。

# 上位50のDEGのヒートマップ
library(pheatmap)
# 有意なDEGをフィルタリング（シンプル化した条件を使用）
sigGenes <- subset(res_df, padj < 0.05 & abs(log2FoldChange) > 1)
# 上位DEGの選択（遺伝子シンボルを使用）
top_genes <- sigGenes$ensembl_id[1:min(50, nrow(sigGenes))]
# vst変換したカウントデータからヒートマップを作成
vsd <- vst(dds, blind = FALSE)
vst_counts <- assay(vsd)
# アノテーション情報の作成（シンプル化した条件を使用）
annotation_col <- data.frame( Condition = sample_info$simple_condition, row.names = rownames(sample_info)
)
# 遺伝子行名のマッピング作成（ヒートマップのラベル用）
gene_labels <- ensembl_to_symbol(top_genes)
names(gene_labels) <- top_genes
# ヒートマップの作成
pheatmap(vst_counts[top_genes, ], scale = "row", # Z-scoreでスケーリング annotation_col = annotation_col, show_rownames = TRUE, labels_row = gene_labels, # 遺伝子シンボルをラベルとして使用 main = "Heatmap of Top DEGs", cluster_rows = TRUE, cluster_cols = TRUE, color = colorRampPalette(c("navy", "white", "firebrick3"))(50))

このヒートマップ作成コードでは、以下の重要なステップを実行しています：

1. 対象遺伝子の選択：
   • 統計的に有意（padj < 0.05）かつ生物学的に意味のある変化（|log2FoldChange| > 1）を示す遺伝子をフィルタリングします。
   • フィルタリングされた遺伝子の中から上位50遺伝子を選択します。min(50, nrow(sigGenes))を使用して、有意な遺伝子が50未満の場合にもエラーが発生しないように配慮しています。
2. 発現データの正規化：
   • vst関数（variance stabilizing transformation）を使用して、カウントデータを分散安定化変換します。
   • blind = FALSEオプションは、デザイン情報（群間の差）を考慮して変換することを指定します。
   • assay関数で変換済みのカウントデータを取得します。
3. アノテーション情報の準備：
   • annotation_colデータフレームを作成し、サンプルの条件情報を指定します。
   • 行名をsample_infoのサンプルIDに設定することで、ヒートマップの列とアノテーション情報が正しく対応するようにします。
4. 遺伝子ラベルの準備：
   • ensembl_to_symbol関数で選択した遺伝子のEnsembl IDを遺伝子シンボルに変換します。
   • names関数で、各シンボルにEnsembl IDを対応させた名前付きベクトルを作成します。
5. pheatmap関数の主要パラメータ：
   • 第一引数：VST変換済みの発現データから選択した上位遺伝子のみを抽出します。
   • scale = "row"：各行（遺伝子）ごとにZ-scoreでスケーリングします。これにより、遺伝子間の発現量の絶対値の違いではなく、発現パターンの違いに焦点を当てることができます。
   • annotation_col：サンプルのアノテーション（条件など）を色付きのバーで表示します。
   • show_rownames = TRUE：行名（遺伝子名）を表示します。
   • labels_row：行名として使用するカスタムラベル（遺伝子シンボル）を指定します。
   • cluster_rows, cluster_cols：行と列をそれぞれ階層的クラスタリングするかどうかを指定します。
   • color：ヒートマップの色スケールを指定します。colorRampPalette関数で、青（低発現）から白、そして赤（高発現）へのグラデーションを50段階で作成します。

ヒートマップを使用することで、以下の情報を視覚的に把握することができます：

• サンプル間の遺伝子発現パターンの違いと類似性
• 実験条件（Treatment vs Control）に関連した発現パターン
• 協調して発現が変化する遺伝子グループ（遺伝子クラスター）

発現変動解析の結果を視覚化する際には、目的に応じて適切な方法を選択することが重要です。MAプロットは全体的な発現変動パターンを把握するのに、火山プロットは統計的・生物学的に重要な遺伝子を特定するのに、ヒートマップは複数サンプル間の発現パターンや遺伝子間の関係を理解するのに、それぞれ適しています。

最後に

いかがだったでしょうか。

可視化は、発現変動解析の結果を直感的に理解し、生物学的な意味を引き出すための非常に強力なツールです。

ぜひMAプロット、火山プロット、ヒートマップを使いこなして、解析結果の説得力を高めていきましょう。