RでGO enrichmentをする方法

RでedgeRを使って差発現遺伝子解析を行った後、GOエンリッチメント解析を行う方法はいくつかあります。ここでは、代表的な方法とその手順を説明します。

必要なパッケージのインストール

まず、必要なパッケージをインストールします。

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install(c("edgeR", "GO.db", "topGO", "org.Hs.eg.db", "AnnotationDbi")) # 例: ヒトの場合

edgeR: 差発現遺伝子解析に使用
GO.db: Gene Ontologyデータベース
topGO: GOエンリッチメント解析に使用
org.Hs.eg.db (例): ヒトの遺伝子IDとEntrez IDをマッピングするアノテーションパッケージ。対象生物に合わせて変更してください（例: org.Mm.eg.db (マウス), org.Rn.eg.db (ラット)）。
AnnotationDbi: アノテーションパッケージを扱うための共通インターフェース

GOエンリッチメント解析のステップ

差発現遺伝子リストの準備:
- edgeRの結果から、有意に発現変動した遺伝子のリストを抽出します。
- 一般的に、FDR（調整済みp値）が一定の閾値（例: 0.05）を下回る遺伝子を抽出します。
- edgeRの結果オブジェクト（例: et) から必要な情報を抽出します。
```
# edgeRの結果オブジェクト et があると仮定
de_genes <- rownames(et$table)[et$table$FDR < 0.05]

head(de_genes) # 抽出された遺伝子リストを確認
```
遺伝子IDの変換:
- 抽出された遺伝子ID（例: エンサンブルID、遺伝子シンボル）を、Entrez IDに変換します。
- Entrez IDはGOアノテーションに使用される標準的なIDです。
- AnnotationDbiパッケージの機能を使って変換を行います。
```
library(AnnotationDbi)
library(org.Hs.eg.db) # 例: ヒトの場合

gene_symbols <- de_genes # 仮に抽出された遺伝子が遺伝子シンボルだった場合
entrez_ids <- mapIds(org.Hs.eg.db, 
                    keys = gene_symbols,
                    keytype = "SYMBOL", # シンボルでマッピングする場合
                    column = "ENTREZID")

entrez_ids <- na.omit(entrez_ids) # NAを取り除く

head(entrez_ids) # 変換されたEntrez IDを確認
```
- 注意点:
  - keytype には、入力された遺伝子IDの種類を指定します。（例: “SYMBOL”, “ENSEMBL”）
  - column には、出力したいIDの種類を指定します。（ここでは、Entrez ID）
  - NA（マッピングできなかった）IDは na.omit() で取り除きます。

全遺伝子リストの準備:

edgeR解析に使用した全遺伝子のリストをEntrez IDに変換します。
GOエンリッチメント解析では、全遺伝子リストを背景（バックグラウンド）遺伝子として用います。

all_genes <- rownames(et$table)
all_entrez_ids <- mapIds(org.Hs.eg.db, 
                       keys = all_genes,
                       keytype = "SYMBOL", # シンボルでマッピングする場合
                       column = "ENTREZID")
all_entrez_ids <- na.omit(all_entrez_ids)

head(all_entrez_ids)

GOエンリッチメント解析:

topGOパッケージを使って、GOエンリッチメント解析を行います。
topGOでは、各GOタームに対して、有意にエンリッチされているかを統計的に評価します。

library(topGO)

# GOカテゴリー (BP: 生物学的プロセス, CC: 細胞構成要素, MF: 分子機能)
go_categories <- c("BP", "CC", "MF")

for(category in go_categories) {
    gene_list <- factor(as.integer(all_entrez_ids %in% entrez_ids))
    names(gene_list) <- all_entrez_ids
    
    go_data <- new("topGOdata",
                  ontology = category,
                  allGenes = gene_list,
                  geneSel = function(p) p == 1,
                  annot = annFUN.org,
                  mapping = "org.Hs.eg.db",
                  ID = "entrez")
    
    result_fisher <- runTest(go_data, algorithm = "classic", statistic = "fisher")
    
    go_results <- GenTable(go_data, 
                          classicFisher = result_fisher, 
                          topNodes = 20)
    print(paste("GO enrichment for category: ", category))
    print(go_results)
}

ontology には、解析したいGOカテゴリー（“BP”, “CC”, “MF”）を指定します。
allGenes には、全遺伝子リストの0/1のベクトルを与えます。（1は差発現遺伝子、0はそうでない）
geneSel は、差発現遺伝子を選択する関数を指定します。
annFUN.org は、アノテーション関数を指定します。
mapping には、アノテーションパッケージを指定します。
ID には、IDの種類（Entrez）を指定します。
runTest で、Fisher’s exact test を実行します。
GenTable で、結果を表形式で表示します。topNodes は、表示するGOタームの数です。

結果の解釈

GenTable の結果には、各GOタームに対するp値（classicFisher）が表示されます。
p値が小さいほど、そのGOタームが差発現遺伝子群で有意にエンリッチされていることを示します。
結果を基に、どのGOタームが差発現遺伝子の機能に関わっているかを解釈します。
必要に応じて、GO DAG（有向非巡回グラフ）を可視化して、結果をより深く理解します。

補足

上記のコード例は、ヒトの場合を想定していますが、対象生物に合わせてアノテーションパッケージを変更してください。
topGO以外にも、clusterProfiler パッケージもよく使用されます。
必要に応じて、FDR補正やその他の統計的検定も適用してください。
GOエンリッチメント解析は、あくまで遺伝子機能の傾向を把握するためのツールであり、結果の解釈には生物学的知識が不可欠です。
全遺伝子リストの定義（例えば、RNA-seqデータの全遺伝子か、アノテーションされた遺伝子のみか）によって、結果が異なる場合があります。

ご自身の解析データや目的に合わせて、上記の手順を適宜調整してください。
もし、具体的なデータやコードについて質問があれば、お気軽にお尋ねください。

Mouseの場合

マウス（Mus musculus）の場合のGOエンリッチメント解析も、基本的な流れはヒトの場合と変わりませんが、いくつか変更が必要な点があります。以下に、マウスの場合の具体的な手順とコード例を示します。

必要なパッケージのインストール（マウス版）

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install(c("edgeR", "GO.db", "topGO", "org.Mm.eg.db", "AnnotationDbi")) # マウスの場合

org.Mm.eg.db: マウスの遺伝子IDとEntrez IDをマッピングするアノテーションパッケージ。

GOエンリッチメント解析のステップ（マウス版）

差発現遺伝子リストの準備:

ヒトの場合と同様に、edgeRの結果から、有意に発現変動した遺伝子のリストを抽出します。

# edgeRの結果オブジェクト et があると仮定
de_genes <- rownames(et$table)[et$table$FDR < 0.05]
head(de_genes) # 抽出された遺伝子リストを確認

遺伝子IDの変換:

抽出された遺伝子IDをEntrez IDに変換します。
アノテーションパッケージを org.Mm.eg.db に変更します。

library(AnnotationDbi)
library(org.Mm.eg.db) # マウスのアノテーションパッケージ

gene_symbols <- de_genes # 仮に抽出された遺伝子が遺伝子シンボルだった場合
entrez_ids <- mapIds(org.Mm.eg.db, 
                    keys = gene_symbols,
                    keytype = "SYMBOL", # シンボルでマッピングする場合
                    column = "ENTREZID")

entrez_ids <- na.omit(entrez_ids) # NAを取り除く
head(entrez_ids) # 変換されたEntrez IDを確認

全遺伝子リストの準備:

edgeR解析に使用した全遺伝子のリストをEntrez IDに変換します。
アノテーションパッケージを org.Mm.eg.db に変更します。

all_genes <- rownames(et$table)
all_entrez_ids <- mapIds(org.Mm.eg.db, 
                       keys = all_genes,
                       keytype = "SYMBOL", # シンボルでマッピングする場合
                       column = "ENTREZID")
all_entrez_ids <- na.omit(all_entrez_ids)

head(all_entrez_ids)

GOエンリッチメント解析:

topGOパッケージを使って、GOエンリッチメント解析を行います。
topGOdata の mapping に、マウスのアノテーションパッケージ org.Mm.eg.db を指定します。

library(topGO)

# GOカテゴリー (BP: 生物学的プロセス, CC: 細胞構成要素, MF: 分子機能)
go_categories <- c("BP", "CC", "MF")

for(category in go_categories) {
    gene_list <- factor(as.integer(all_entrez_ids %in% entrez_ids))
    names(gene_list) <- all_entrez_ids
    
    go_data <- new("topGOdata",
                  ontology = category,
                  allGenes = gene_list,
                  geneSel = function(p) p == 1,
                  annot = annFUN.org,
                  mapping = "org.Mm.eg.db",  # マウスのアノテーションパッケージ
                  ID = "entrez")
    
    result_fisher <- runTest(go_data, algorithm = "classic", statistic = "fisher")
    
    go_results <- GenTable(go_data, 
                          classicFisher = result_fisher, 
                          topNodes = 20)
    print(paste("GO enrichment for category: ", category))
    print(go_results)
}

変更点まとめ

アノテーションパッケージの変更:
- org.Hs.eg.db の代わりに org.Mm.eg.db を使用します。
その他:
- 上記以外は、基本的にはヒトの場合と同じです。
- keytype には、抽出された遺伝子IDの種類（例: “SYMBOL”, “ENSEMBL”）を指定する必要があります。
- 遺伝子シンボルではなく、エンサンブルIDを使っている場合も同様にkeytypeを設定します。

補足

他の生物種でも、同様に、対応するアノテーションパッケージを使用することで、GOエンリッチメント解析を行うことができます。
結果の解釈は、マウスの生物学的な知識に基づいて行う必要があります。
必要に応じて、FDR補正やその他の統計的検定も適用してください。
clusterProfiler パッケージなど、他のGOエンリッチメント解析ツールも利用できます。

上記コードは、マウスの遺伝子シンボルからEntrez IDに変換する場合の例です。エンサンブルIDから変換する場合は、 keytype を “ENSEMBL” に変更してください。

もし、具体的なデータやコードについて質問があれば、お気軽にお尋ねください。

このブログを検索/Search this blog

Translate

石川昌和 | Masakazu Ishikawa | いしかわまさかず

RでGO enrichmentをする方法

Mouseの場合

アーカイブ