ご無沙汰をしております、受託コンサルティングチームの kutsuy でございます。 今回はGSEA（Gene Set Enrichment Analysis）解析について実際にコマンドを打ちながら内容を理解していきたいと思います。 GSEA解析では、特定の機能を持つ遺伝子グループに着目し、二群間で個々の遺伝子の発現がどのように変化したかをグループ全体の傾向として捉えることができます。
動作環境
- R version 4.3.0
- airway v1.20.0
- DESeq2 v1.40.2
- clusterProfiler v4.8.3
- org.Hs.eg.db v3.17.0
- dplyr v1.1.4
今回は公開データとして R パッケージにもなっている airway データを使用したいと思います。 私はツールの動作確認などをする際によくこのairway データを使います。 ヒト気道平滑筋細胞に対し Dexamethasone（ステロイド）で処理したという実験系になっていて、 処理有りが 4 サンプル、処理無しが 4 サンプル、合計 8 サンプルのデータです。
airway のインストール
まず、R パッケージのインストールとライブラリの読み込みです。
> if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
> BiocManager::install("airway")
> library(airway)
データをオブジェクトとして取得し、内容を確認してみます。
> data("airway")
> airway
class: RangedSummarizedExperiment
dim: 63677 8
metadata(1): ''
assays(1): counts
rownames(63677): ENSG00000000003 ENSG00000000005 ... ENSG00000273492
ENSG00000273493
rowData names(10): gene_id gene_name ... seq_coord_system symbol
colnames(8): SRR1039508 SRR1039509 ... SRR1039520 SRR1039521
colData names(9): SampleName cell ... Sample BioSample
データ型は RangedSummarizedExperiment です（長い…）。 こちらの詳細については割愛しますが、発現カウントデータ、各遺伝子のゲノム情報、サンプル情報などが一体となっているオブジェクトです。 カウントデータは、遺伝子に相当する行が 63677、サンプルに相当する列が 8 であることがわかりました。
発現マトリクスをみてみましょう。assay 関数で取得できます。 データが大きいため、head 関数で上位 6 行のみを表示してみます。
> head(assay(airway))
SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517 SRR1039520 SRR1039521
ENSG00000000003 679 448 873 408 1138 1047 770 572
ENSG00000000005 0 0 0 0 0 0 0 0
ENSG00000000419 467 515 621 365 587 799 417 508
ENSG00000000457 260 211 263 164 245 331 233 229
ENSG00000000460 60 55 40 35 78 63 76 60
ENSG00000000938 0 0 2 0 1 0 0 0
行名は ENSG から始まる Ensembl Gene ID になっていますね。 列名は SRR から始まる Sequence Read Archive（SRA）のラン番号になっています。 SRAから生データ（fastq.gz）も取得できますので、最初から再解析も可能です。
二群間比較解析
次に、二群間比較解析を実施します。 今回は DESeq2 というツールを用いますが、edgeR など他のツールでも構いません。
> library(DESeq2) # 処理条件の順序設定 > airway$dex <- relevel(airway$dex, "untrt") # DESeq2による比較解析 > dds <- DESeqDataSet(airway, design = ~ dex) > dds <- DESeq(dds) estimating size factors estimating dispersions gene-wise dispersion estimates mean-dispersion relationship final dispersion estimates fitting model and testing > res <- results(dds)
結果をみてみましょう。オブジェクト res は DESeqResults というデータ型になっています。
> res
log2 fold change (MLE): dex trt vs untrt
Wald test p-value: dex trt vs untrt
DataFrame with 63677 rows and 6 columns
baseMean log2FoldChange lfcSE stat pvalue padj
<numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
ENSG00000000003 708.6022 -0.3788467 0.173141 -2.188080 0.0286637 0.139309
ENSG00000000005 0.0000 NA NA NA NA NA
ENSG00000000419 520.2979 0.2037605 0.100599 2.025479 0.0428182 0.183359
ENSG00000000457 237.1630 0.0340429 0.126279 0.269585 0.7874794 0.930569
ENSG00000000460 57.9326 -0.1171781 0.301237 -0.388990 0.6972833 0.895443
... ... ... ... ... ... ...
ENSG00000273489 0.275899 1.545726 3.52279 0.438779 0.660822 NA
ENSG00000273490 0.000000 NA NA NA NA NA
ENSG00000273491 0.000000 NA NA NA NA NA
ENSG00000273492 0.105978 -0.540733 3.53381 -0.153017 0.878385 NA
ENSG00000273493 0.106142 -0.540733 3.53381 -0.153017 0.878385 NA
全遺伝子IDに対する比較解析の結果が格納されています。
値がNAのセルもありますが、元々発現カウントに0が含まれる遺伝子などは計算が出来ない場合があります。
ここで必ず認識しておくべきポイントをお伝えしたいと思います。
処理条件の順序設定では、airway データセットの dex 列（dexamethasone 処理の有無）における因子の基準レベル（reference level）を "untrt" に設定しました。
この場合、DESeq2 は "untrt" を基準として "trt" との比較（trt vs untrt）を行います。
例えば、ENSG00000000003 では logFC が -0.3788467 です。
負の値になっているので「trt（処理有り）で発現が低下」していることを意味します。
この解釈が曖昧になっていると発現変化の解釈がすべて逆になりますのでくれぐれもご注意ください。
エンリッチメント解析の準備
それではいよいよエンリッチメント解析を実施しますが、 その前に、ORAとGSEAの違いについて確認しておきます。 ORA（Over-Representation Analysis）は「発現変動が有意な遺伝子IDのみ」を入力とします。 例えば、logFC > 1 や padj < 0.05 のように任意の閾値で事前にフィルタリングを実施します。 一方、GSEA（Gene Set Enrichment Analysis）は「全遺伝子の順位付きリスト」を入力とします。 例えば、logFC が大きい順（発現差が大きい順）にソートしたリストなどです。
それではまず事前準備をします。 解析実行のため、Ensembl Gene ID を Entrez ID に変換する必要があります。
> library(clusterProfiler)
> library(org.Hs.eg.db)
> library(dplyr)
# Entrez ID に変換
> gene.df <- bitr(rownames(res), fromType = "ENSEMBL",
toType = "ENTREZID", OrgDb = org.Hs.eg.db)
'select()' returned 1:many mapping between keys and columns
Warning message:
In bitr(rownames(res), fromType = "ENSEMBL", toType = "ENTREZID", :
44.96% of input gene IDs are fail to map...
# 結果と Entrez ID を結合
> res$ENSEMBL <- rownames(res)
> res2 <- inner_join(as.data.frame(res), gene.df, by = "ENSEMBL")
Warning message にも記載がありますが、遺伝子IDの変換で約半分くらいは変換できないようです。 念のため、行列数と head の結果を確認しておきます。
> dim(res2)
[1] 37362 8
> head(res2)
baseMean log2FoldChange lfcSE stat pvalue padj ENSEMBL ENTREZID
1 708.6021697 -0.37884667 0.1731411 -2.1880804 0.02866375 0.1393085 ENSG00000000003 7105
2 0.0000000 NA NA NA NA NA ENSG00000000005 64102
3 520.2979006 0.20376048 0.1005987 2.0254789 0.04281822 0.1833587 ENSG00000000419 8813
4 237.1630368 0.03404294 0.1262790 0.2695852 0.78747941 0.9305689 ENSG00000000457 57147
5 57.9326331 -0.11717811 0.3012365 -0.3889904 0.69728327 0.8954426 ENSG00000000460 55732
6 0.3180984 -1.72456894 3.4936334 -0.4936319 0.62156615 NA ENSG00000000938 2268
Entrez ID に変換された遺伝子ID 37362 個が最も右の列に挿入されていることがわかります。
これでエンリッチメント解析の準備が整いました。 次回はこのデータを使って実際に ORA と GSEA を実行していきたいと思います。
それではまた次回お会いしましょう、皆様ごきげんよう。
▼後編はこちら▼
