皆様、ご無沙汰をしております。アメリエフ解析チームのkutsuyでございます。
今回はタイトルのままなのですが、GO解析結果の可視化について記事を書いております。
まずはGO解析について簡単に補足致しますと、例えばRNAseqデータの発現比較解析を実施して、二群間で有意に発現が変動している遺伝子セットがわかったとします（例えば100個とか）。
これらの遺伝子セットにどのような機能がエンリッチ（濃縮）されているかを検出するのがエンリッチメント解析です。エンリッチメント解析にもいくつか種類がありますが、Gene Ontology データベースを使用するのが通称GO解析になります。

GO解析の結果は基本的には表形式で出力されます。
このデータの特徴を論文などで読者にわかりやすく伝えるには、意味のある可視化が効果的です。以下のサイトにきれいな可視化の図が紹介されていましたので、今回引用させていただきます。

アメリエフの受託解析では生物学的解釈を助けるGO解析まで基本対応とし、受託解析のサービスを提供しております。あえて「データ解析が得意なアメリエフ」だからこそ持っているノウハウも一挙に公開します。コマンド一つ一つ丁寧に説明してみますので、皆様もRを開いてポチっと練習してみてください。
それではいってみましょう。

動作環境

R version 4.3.0
- DOSE v3.26.2
- enrichplot v1.20.3

データの準備

まずは Bar Plot ですね。コマンド先頭の「>」の記号ですが、これはRのコンソールでの入力を示すプロンプトです。実際にコマンドを実行するときは重ねて入力しないようにご注意ください。「>」で始まる行以外は標準出力を掲載しておりますので、入力する必要はありません。

# DOSEパッケージの読み込み
> library(DOSE)
DOSE v3.26.2  For help: https://yulab-smu.top/biomedical-knowledge-mining-book/

If you use DOSE in published research, please cite:
Guangchuang Yu, Li-Gen Wang, Guang-Rong Yan, Qing-Yu He. DOSE: an R/Bioconductor package for Disease Ontology Semantic and Enrichment analysis. Bioinformatics 2015, 31(4):608-609

DOSEパッケージは疾患系のテストデータとなっております。
上記のコマンドを実行すると、「DOSEパッケージが見当たりません」を意味するエラーが出る場合があります。そのような場合はまだDOSEパッケージがインストールされていない状態になっていますので、以下のコマンドを実行してDOSEパッケージをインストールしてください。

ちなみに以下では、DOSEパッケージをインストールするために必要なBiocManagerがインストールされていなければBiocManagerを先にインストールし、BiocManagerを使ってDOSEパッケージをインストールしております（長い…）。

# DOSEパッケージのインストール
> if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

> BiocManager::install("DOSE")

次に、DOSEパッケージに含まれるgeneListというサンプルデータを読み込みます。

> data(geneList)

特にエラーも出ずに次のコマンド待機状態に移行すれば成功しています。つまり、geneListという名前のオブジェクトがインポートされました。ここでgeneListデータの中身を確認しておきましょう。解析を進めるうえで「自分が何をしているのか」を理解して進めることは非常に重要です。くどいようですが、解析を進めるうえで「自分が何をしているのか」を理解して進めることは非常に重要です。

> class(geneList)
[1] "numeric"

まず、データ構造を確認しておきましょう。 class関数はオブジェクトのクラス（データ型）を確認するための関数です。 "numeric"と出力されましたので、geneListオブジェクトは数値型のベクトルであることがわかります。実際にどのようなデータなのか調べてみましょう。

> head(geneList)
    4312     8318    10874    55143    55388      991 
4.572613 4.514594 4.418218 4.144075 3.876258 3.677857

head関数は、ベクトルやデータフレームなどのオブジェクトの先頭部分を表示する関数です。データの中身をざっと確認したい時に使用します。私もよく使います。 head関数のデフォルトは6なので、先頭6個のデータが表示されます。出力が2行に分かれて表示されました。実際のデータ行は2列目になります。 1行目は2列目のデータの名前に相当するもので、今回は遺伝子IDに該当します。遺伝子IDについては説明が多くなるので割愛しますが、またどこかでお伝えできれば…と思います。

さらにgeneListオブジェクトのベクトルの長さを確認しておきましょう。

> length(geneList)
[1] 12495

先ほど先頭6個のデータを確認しましたが、同様のデータが全部で12495個存在することがわかりました。もしすべてのデータを確認する場合は次のようにオブジェクト名のみを入力し実行してください。出力が膨大になるので掲載は割愛させていただきます。

> geneList

データの抽出

ここまでで、テストデータがどのようなデータであるかを確認してきました。それではようやく次のコマンドに進みたいと思います。

> de <- names(geneList)[abs(geneList) > 2]

こちらのコマンドについても分解して説明してみます。
まず、abs(geneList)はgenelistオブジェクトのスコア（先ほど表示した2行目）の絶対値を返します。現状では発現変動が更新している場合は正の値に、減衰している場合は負の値になっているので、絶対値に変換することにより（正負問わず）発現が大きい遺伝子を抽出する意図があります。
次に、[abs(geneList) > 2]の部分でスコアの絶対値が2より大きいか否かの論理値を取得します。論理値とはTRUEあるいはFALSEで表されますので、1個目の遺伝子のスコアが条件を満たすか、満たす場合はTRUE、満たさない場合はFALSEの値を出力します。これを12495個分判定します。
さらに、names(geneList)はgeneListオブジェクトの名前部分（先ほど表示した1行目）を返します。よって、names(geneList)[abs(geneList) > 2]は「geneListオブジェクトのスコアの絶対値が2より大きい遺伝子の名前を抽出する」というコマンドです。この出力を代入演算子（<-）で新たな名前のオブジェクトdeに代入しています。

確認のためにdeオブジェクトのデータを確認しておきましょう。

> class(de)
[1] "character"

オブジェクトdeのデータ構造がcharacterなので文字列のベクトルになっています。

> head(de)
[1] "4312"  "8318"  "10874" "55143" "55388" "991"

それぞれの要素がダブルクオーテーションで囲まれていることに注意してください。各要素は数値（numeric）ではなく文字列（character）です。以降の解析で実際にエンリッチメント解析を実施していくわけですが、関数によって受け付けるデータ構造がベクトル（numeric）、ベクトル（character）、データフレーム、マトリクスなど、厳密に決められていることがほとんどですので、データ構造の確認と理解は非常に重要です。

> length(de)
[1] 207

オブジェクトdeの要素数は207でした。geneListの全データ12495個のうち、先ほどの条件（絶対値が2より大きい）を満たす遺伝子が207個であることがわかります。この207個の遺伝子は発現変動のスコアが大きいので、今回の実験で特徴的な遺伝子セットになります。この遺伝子セットにどのような機能がエンリッチされているかを見たい、というモチベーションになります。

エンリッチメント解析

で、次のコマンドですね。長くなっていますがもう少しお付き合いください。

> edo <- enrichDGN(de)

ここで先ほど用意した207個の遺伝子名が記載されたベクトルに対し、DOSEパッケージに含まれるenrichDGN関数を用いてエンリッチメント解析を実施しています。特にエラーなく次のRの待機状態に進めばOKです。詳細は割愛しますが、DisGeNET（DGN）データベースに基づいて解析を実施しています。

ここで新たに作成したedoオブジェクトを確認してみましょう。

> edo
#
# over-representation test
#
#...@organism    Homo sapiens 
#...@ontology    DisGeNET 
#...@keytype     ENTREZID 
#...@gene    chr [1:207] "4312" "8318" "10874" "55143" "55388" "991" "6280" "2305" ...
#...pvalues adjusted by 'BH' with cutoff <0.05 
#...247 enriched terms found
'data.frame':   247 obs. of  9 variables:
 $ ID         : chr  "C0007124" "C3642347" "C1609538" "C1827293" ...
 $ Description: chr  "Noninfiltrating Intraductal Carcinoma" "Basal-Like Breast Carcinoma" "Latent Tuberculosis" "Carcinoma of urinary bladder, invasive" ...
 $ GeneRatio  : chr  "24/204" "16/204" "14/204" "16/204" ...
 $ BgRatio    : chr  "486/21671" "245/21671" "183/21671" "287/21671" ...
 $ pvalue     : num  3.81e-11 1.56e-09 2.21e-09 1.49e-08 2.05e-08 ...
 $ p.adjust   : num  1.57e-07 3.04e-06 3.04e-06 1.54e-05 1.69e-05 ...
 $ qvalue     : num  1.27e-07 2.46e-06 2.46e-06 1.25e-05 1.37e-05 ...
 $ geneID     : chr  "4312/6280/6279/7153/6278/9582/55872/3669/4751/5080/2146/1111/3887/6790/9370/125/8839/2066/3169/9547/10647/2625/7021/5241" "2305/1062/4605/9833/7368/11065/10232/55765/5163/2146/2568/3620/6790/3169/2625/5241" "8685/3627/3669/6373/820/3002/50852/3902/3620/6355/10578/9370/3117/23158" "2305/6279/55165/220134/2146/9212/4321/6790/891/8544/9122/23158/652/2625/5241/10551" ...
 $ Count      : int  24 16 14 16 13 17 10 7 19 19 ...
#...Citation
 T Wu, E Hu, S Xu, M Chen, P Guo, Z Dai, T Feng, L Zhou, W Tang, L Zhan, X Fu, S Liu, X Bo, and G Yu.
 clusterProfiler 4.0: A universal enrichment tool for interpreting omics data.
 The Innovation. 2021, 2(3):100141

このオブジェクトにはエンリッチメント解析の結果が格納されています。こちらも一つ一つ説明するとかなり長くなりますので、「enrichDGN」で検索して公式ドキュメントなどを確認してみてください。とても理解が深まると思います。そしてついに最後のbar plotの可視化です。

bar plot の出力

> library(enrichplot)
> barplot(edo, showCategory=20)

まずは、library(enrichplot)でenrichplotパッケージを読み込みます。 DOSEパッケージと同様に、「enrichplotパッケージがありません」のようなエラーが出る場合は、以下のコマンドを試していただき再度上記のコマンドを実行してみてください。

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

> BiocManager::install("enrichplot")

実際にbarplotを描画するコマンドが、barplot(edo, showCategory=20) になります。私の手元ではこちらで以下の図が無事に出力されました。

p.adjust（補正p値）が小さい順に、つまり有意差が大きい機能の順にタームが上から並んでいます。また、bar plot がp.adjustの数値によって赤から青に色合いがグラデーションとなっており、各タームに含まれる遺伝子数が横軸に表現されています。二群間でどのような機能がエンリッチされているかとてもわかりやすいですね。

まとめ

今回はテストデータの確認とbar plotの描画を実施してみました。おそらく研究されている読者の方は自分のデータでも同様の解析を実施したいことが多くあると思います。データを丁寧に確認しながら進めることで自力でも解析が出来るようになると思いますのでぜひトライしてみてください。また、かっこいい可視化の方法があれはぜひ共有いただけると嬉しいです。

では次回は続きをやっていきますのでぜひご覧ください。