【RNA-seq】GSEAで遺伝子セットエンリッチメント解析を行う方法【R】

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

GSEA解析とは？

GSEA（Gene Set Enrichment Analysis） は、発現変動解析で得られた遺伝子のリストを用いて、生物学的な機能や経路に関与する「遺伝子セット（gene sets）」が全体として有意に変動しているかを統計的に評価する手法です。

通常のエンリッチメント解析と異なり、GSEAは全遺伝子のランキング情報を使用するため、特に微細な発現変化を捉えることができ、発現変動傾向の全体的なシフトを検出するのに適しています。

単一の遺伝子ではなく、機能的なグループ単位（経路・ネットワーク）で解析したいときに有効で発現変動が小さくて個別には検出されない遺伝子群も、集合的に見ることで意味のある変動が見つかることがあります。

解析準備

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"))

これで準備完了です。

GSEA解析のやり方

こちらが実装コードになります。

# ランキングの作成
# log2FoldChangeとp値からrankingを計算（p値の符号付き-log10変換）
res_df$ranking <- sign(res_df$log2FoldChange) * -log10(res_df$pvalue)
ranking <- res_df$ranking
names(ranking) <- rownames(res_df)
ranking <- sort(ranking, decreasing = TRUE)
# GSEAの実行
# ENSEMBL IDからEntrez IDへの変換（GSEA用）
all_genes <- data.frame( ENSEMBL = rownames(res_df), stringsAsFactors = FALSE
)
all_entrez <- bitr(all_genes$ENSEMBL, fromType = "ENSEMBL", toType = "ENTREZID", OrgDb = org.Mm.eg.db)
# EntrezIDでのランキングの作成
ranking_entrez <- ranking[names(ranking) %in% all_entrez$ENSEMBL]
mapping <- all_entrez$ENTREZID
names(mapping) <- all_entrez$ENSEMBL
names(ranking_entrez) <- mapping[names(ranking_entrez)]
ranking_entrez <- sort(ranking_entrez, decreasing = TRUE)
# 重複するEntrez IDをチェックして除去
duplicated_ids <- names(ranking_entrez)[duplicated(names(ranking_entrez))]
if(length(duplicated_ids) > 0) { message(paste("重複するEntrez ID:", length(duplicated_ids), "個")) # 重複を削除（各Entrez IDの最初の出現のみを保持） ranking_entrez <- ranking_entrez[!duplicated(names(ranking_entrez))] # 代替方法: 重複するEntrez IDごとに最大値を保持する場合 # dup_list <- split(ranking_entrez, names(ranking_entrez)) # ranking_entrez <- sapply(dup_list, max)
}
# GSEA解析の実行
gsea_result <- gseKEGG(geneList = ranking_entrez, organism = "mmu", pvalueCutoff = 0.05, pAdjustMethod = "BH", nPerm = 1000, # パーミュテーション回数 minGSSize = 10, # 最小遺伝子セットサイズ maxGSSize = 500) # 最大遺伝子セットサイズ

GSEAの実行コードは、遺伝子のランキングの作成から始まり、KEGG経路に対するエンリッチメント解析まで行います。以下に詳細なステップを説明します：

1. 遺伝子ランキングの作成：
   • GSEAでは、すべての遺伝子を重要性または変動の大きさでランク付けする必要があります。
   • ここでは、log2FoldChangeの方向（符号）とp値の有意性（-log10変換）を組み合わせた複合スコアを使用しています。
   • sign(res_df$log2FoldChange) * -log10(res_df$pvalue)：
     • sign()関数は値の符号（正または負）を返します。
     • log10(p値)はp値を対数スケールに変換し、小さなp値（高い有意性）を大きな値に変換します。
     • 両者を掛け合わせることで、有意に上方制御された遺伝子は大きな正の値、有意に下方制御された遺伝子は大きな負の値を持つようになります。
   • 名前付きベクトルとして格納し、値の大きい順（降順）にソートします。
2. 遺伝子IDの変換：
   • KEGG経路のGSEAにはEntrez IDが必要なため、Ensembl IDからEntrez IDへの変換を行います。
   • すべての遺伝子のIDを含むデータフレームを作成し、bitr関数で変換します。
   • 変換後、名前付きベクトルを作成して、Ensembl IDをキーとしてEntrez IDを値とするマッピングを構築します。
   • このマッピングを使用して、ランキングベクトルの名前（Ensembl ID）をEntrez IDに置き換えます。
   • 最後に、変換されたランキングを再度降順にソートします。
3. 重複Entrez IDの処理：
   • 複数のEnsembl IDが同じEntrez IDにマッピングされることがあり、これが重複を引き起こします。
   • 重複するEntrez IDを検出し、最初に出現したものだけを保持します。
   • 代替方法として、各重複IDについて最大値を採用する方法もコメントとして示しています。
4. GSEAの実行：
   • gseKEGG関数を使用して、KEGGパスウェイに対するGSEA解析を実行します。
   • 主要パラメータ：
     • geneList：ランク付けされた遺伝子リスト（Entrez ID）
     • organism：生物種のKEGGコード（“mmu” = マウス）
     • pvalueCutoff：有意性の閾値（0.05）
     • pAdjustMethod：多重検定補正の方法（“BH” = Benjamini-Hochberg法）
     • nPerm：パーミュテーション（並べ替え）の回数（1000）- 統計的有意性を評価するために使用
     • minGSSize、maxGSSize：解析に含める遺伝子セットのサイズ範囲（10〜500遺伝子）

GSEAは、遺伝子セット（パスウェイなど）内の遺伝子がランキングの上位または下位に集中しているかを評価し、エンリッチメントスコア（ES）を計算します。正のESは遺伝子セットが上方制御遺伝子に濃縮されていることを、負のESは下方制御遺伝子に濃縮されていることを示します。

# GSEA結果の表示
head(gsea_result, 10)

> head(gsea_result) ID Description setSize enrichmentScore
mmu03010 mmu03010 Ribosome 144 -0.8345503
mmu00190 mmu00190 Oxidative phosphorylation 122 -0.7906797
mmu05171 mmu05171 Coronavirus disease - COVID-19 206 -0.6855225
mmu05012 mmu05012 Parkinson disease 241 -0.6315631
mmu05208 mmu05208 Chemical carcinogenesis - reactive oxygen species 200 -0.6257951
mmu04932 mmu04932 Non-alcoholic fatty liver disease 139 -0.6440056 NES pvalue p.adjust qvalue rank leading_edge
mmu03010 -2.946786 0.001988072 0.01164072 0.007214275 1551 tags=78%, list=11%, signal=70%
mmu00190 -2.721614 0.002036660 0.01164072 0.007214275 1996 tags=68%, list=15%, signal=59%
mmu05171 -2.534552 0.002016129 0.01164072 0.007214275 1204 tags=42%, list=9%, signal=39%
mmu05012 -2.367968 0.002044990 0.01164072 0.007214275 2200 tags=45%, list=16%, signal=38%
mmu05208 -2.307132 0.002024291 0.01164072 0.007214275 1770 tags=45%, list=13%, signal=40%
mmu04932 -2.248367 0.002036660 0.01164072 0.007214275 1770 tags=45%, list=13%, signal=40% core_enrichment
mmu03010 78523/27397/66448/78294/68463/68193/100039532/27395/64659/19989/108168114/19896/107732/121022/57294/20088/268449/68436/100042335/66242/100040416/20102/67281/75617/19942/19921/114641/670832/19946/20115/67945/108169197/27176/67891/267019/67025/56284/19941/19981/100503670/100043813/65019/19934/66475/20042/20068/57808/66258/66845/110954/20091/19935/26451/67671/22186/66163/66481/67427/27367/19933/19988/108169053/225215/19982/20084/56040/67681/66292/68537/20090/66489/54217/27370/20044/432502/20103/19899/67115/68052/94066/69527/16898/20055/19951/19944/20005/619547/20116/76846/108167871/20085/319195/94065/269261/20104/22121/14109/66407/66419/76808/27207/11837/270106/27050/100040298/68611/54127/19943/20054/67186/16785/26961/24030
mmu00190 100041273/67264/66290/73834/11974/114143/66414/12866/11946/17992/12859/66445/70316/68194/12867/11949/66142/54405/17993/68202/68197/66694/22273/11957/22272/67530/67130/66052/12864/230075/623286/12858/66925/67184/225887/67942/11958/67680/68198/20463/57423/12865/74776/66152/66416/66594/17991/72900/66218/11950/110323/27425/75406/66144/66576/11951/66091/595136/11972/76429/66108/66046/66495/68349/69875/71679/12857/104130/407785/66377/12868/66916/67273/67126/226646/12861/78330/28080/228033/68375/17995/66043/68342
mmu05171 19989/108168114/19896/57294/20088/268449/68436/100042335/100040416/20102/67281/75617/19942/19921/114641/670832/19946/20115/67945/108169197/27176/67891/267019/67025/19941/19981/100503670/100043813/65019/19934/66475/20042/20068/57808/110954/20091/26451/67671/22186/66481/67427/27367/19933/19988/108169053/225215/19982/20084/56040/20090/66489/54217/27370/20044/432502/20103/19899/67115/68052/16898/20055/19951/19944/20005/619547/54131/20116/76846/108167871/20085/319195/269261/20104/22121/14109/76808/27207/11837/270106/27050/100040298/54127/19943/20054/67186/16785/26961
mmu05012 102631912/11947/22333/22223/12028/19173/100041273/67264/22334/57296/19185/26446/140499/26445/66414/22213/12866/26444/23996/66413/19182/19172/17161/11946/20422/17992/22195/12859/14680/66445/70316/68194/78294/23997/68015/16593/19170/22166/12867/11949/11739/66142/54405/17993/68202/68197/66694/22273/11957/22272/67530/67130/66052/12864/68943/230075/623286/12858/12315/66925/67184/225887/67942/67680/68198/20463/12865/66152/66416/66594/17991/72900/66218/11950/110323/75406/66576/11951/66091/595136/67151/66108/22186/66046/66495/68349/69875/71679/12857/104130/407785/66377/12868/66916/67273/56551/67126/226646/12861/78330/28080/228033/20655/68375/17995/66043/57320/68342
mmu05208 11651/66414/12866/66447/57344/11946/17992/12859/66445/17970/70316/68194/12867/11949/11739/66142/54405/17993/26417/68202/68197/66694/14862/22273/11957/22272/67530/67130/66052/58810/12864/230075/623286/12858/66925/67184/225887/67942/67680/68198/20463/12865/66152/14865/66416/66594/15461/17991/72900/66218/11950/110323/13057/14863/75406/26396/66576/11951/66091/595136/66108/66046/12015/66495/68349/69875/71679/12857/104130/407785/66377/12868/66916/67273/14873/67126/226646/12861/78330/28080/12408/14866/228033/20655/211666/68375/14872/17995/66043/68342
mmu04932 11651/66414/12866/17992/12859/66445/70316/68194/12867/66142/54405/17993/68202/68197/66694/14102/22273/22272/67530/67130/66052/12864/230075/623286/12858/66925/67184/225887/67680/68198/20463/12865/66152/66416/66594/17991/72900/66218/110323/75406/12606/66576/66091/595136/66108/66046/66495/68349/69875/12857/104130/407785/66377/12868/66916/67273/226646/12861/78330/19082/68375/17995/68342

# エンリッチメントプロットの作成
if (!is.null(gsea_result) && nrow(gsea_result) > 0) { # 上位パスウェイのエンリッチメントプロット top_pathway <- gsea_result$ID[1] gseaplot2(gsea_result, geneSetID = top_pathway, title = gsea_result$Description[1])

 # 複数のパスウェイを一度に表示 if (nrow(gsea_result) >= 3) { multi_pathways <- gsea_result$ID[1:3] gseaplot2(gsea_result, geneSetID = multi_pathways, title = "Top 3 Enriched Pathways") }

# リッジプロットでの可視化
ridgeplot(gsea_result, showCategory = 20, label_format = 60) + theme(axis.text.y = element_text(size = 8), # Y軸ラベルのサイズ axis.text.x = element_text(size = 8)) # X軸ラベルのサイズ
# GSEAドットプロットでの可視化
dotplot(gsea_result, showCategory = 10, title = "GSEA Dotplot")

GSEA結果の可視化コードは、clusterProfilerパッケージが提供する複数のビジュアライゼーション関数を使用して、結果を様々な観点から可視化します：

1. 結果の数値表示：
   • head(gsea_result, 10)でGSEA結果の上位10パスウェイを表形式で表示します。主な情報には、パスウェイID、説明、遺伝子セットサイズ、エンリッチメントスコア（ES）、正規化されたES（NES）、p値、調整p値などが含まれます。
2. エラーハンドリング：
   • if (!is.null(gsea_result) && nrow(gsea_result) > 0)を使用して、有意な結果が存在する場合にのみ可視化を実行します。
3. エンリッチメントプロットの作成：
   • gseaplot2関数を使用して、特定のパスウェイのエンリッチメントプロットを生成します。
   • まず最も有意なパスウェイ（gsea_result$ID[1]）を選択し、個別のプロットを作成します。
   • タイトルにはパスウェイの説明（gsea_result$Description[1]）を使用します。
   • エンリッチメントプロットには主に3つの要素が含まれます：
     • 上部：ランキングにおける遺伝子セットのメンバー遺伝子の分布（黒いバーコード）
     • 中部：走行エンリッチメントスコア（Running Enrichment Score）の曲線
     • 下部：ランク付けされた遺伝子リストのスコア分布
4. 複数パスウェイの同時表示：
   • 十分な結果がある場合（nrow(gsea_result) >= 3）、上位3つのパスウェイを同時に1つの図に表示します。
   • gseaplot2関数でgeneSetIDパラメータにパスウェイIDのベクトルを渡し、複数の曲線を描画します。
5. 異なるビジュアライゼーション方法：
   • ridgeplot：パスウェイに関連する遺伝子のランキング分布を密度曲線（リッジライン）として表示します。異なるパスウェイにおける遺伝子の分布パターンを比較するのに役立ちます。
   • dotplot：パスウェイ（y軸）とNES（x軸）を点で表示します。点の大きさは遺伝子数、色はp値を反映します。

これらの可視化手法によって、GSEAの結果を多角的に解釈することができます。特にエンリッチメントプロットは、遺伝子セットのメンバー遺伝子がランキングのどの位置に分布しているかを詳細に確認でき、上方制御と下方制御のパターンを視覚的に理解するのに非常に役立ちます。

最後に

いかがだったでしょうか。ぜひGSEA解析挑戦してみてください！