Gene Set Enrichment Analysis (GSEA)入門
GSEAとは
GSEAは発現差異解析の結果などで得られる遺伝子群がどういった機能のものかを明らかにするために用いられる解析手法です。機能表現として用いられる主なものとしてはKEGGのパスウェイ分類やGene Ontology(GO)があります。
clusterProfilerでGSEA入門
GSEA には様々なGSEAがあります。単に注目する遺伝子のリストをGSEAへの入力とするもの、発現プロファイルを入力とするもの、機能のデータベースやGSEA手法を複数組み合わせたりするもの (EGSEA) など様々です。[余談ですがGeneを代謝物(Metabolite)に置き換えたMSEAもあります。]
ここでは前述したものの内、最初の「単に注目する遺伝子のリストをGSEAへの入力とするもの」をbioconductorのパッケージ clusterProfiler を用いて体験してみましょう。
(clusterProfilerのインストールについては省略します。)
まず注目する遺伝子のリストを下記で用意します。
※編集注 : 以下はすべてRでの実行になります。
> data(geneList, package="DOSE")
> gene <- names(geneList)[abs(geneList) > 2]
> gene
[1] "4312" "8318" "10874" "55143" "55388" "991" "6280" "2305" "9493" "1062" "3868"
[12] "4605" "9833" "9133" "6279" "10403" "8685" "597" "7153" "23397" "6278" "79733"
[23] "259266" "1381" "3627" "27074" "6241" "55165" "9787" "7368" "11065" "55355" "9582"
[34] "220134" "55872" "51203" "3669" "83461" "22974" "10460" "10563" "4751" "6373" "8140"
[45] "79019" "820" "10635" "1844" "4283" "27299" "55839" "27338" "890" "9415" "983"
[56] "54821" "10232" "4085" "6362" "9837" "5080" "7850" "81930" "5918" "81620" "332"
[67] "55765" "79605" "3832" "6286" "5163" "2146" "3002" "50852" "7272" "2568" "64151"
[78] "51806" "366" "2842" "9212" "140578" "51659" "8715" "4902" "8208" "1111" "9319"
[89] "9055" "3833" "146909" "23475" "4321" "11182" "10112" "3902" "3620" "3887" "51514"
[100] "6790" "4521" "891" "57110" "8544" "1448" "24137" "6355" "10578" "4174" "9232"
[111] "643314" "1307" "776" "4129" "9370" "196740" "25924" "8857" "1602" "51161" "3708"
[122] "23090" "10742" "51760" "9122" "10699" "8416" "60598" "79148" "64799" "4629" "1556"
[133] "55096" "26289" "6038" "771" "51313" "23704" "3117" "80129" "23158" "125" "4958"
[144] "4857" "1311" "5105" "5174" "730" "2018" "81563" "2532" "1308" "4250" "23362"
[155] "2167" "51705" "2593" "652" "80736" "4036" "57502" "5507" "56521" "22885" "4137"
[166] "8483" "8839" "2066" "4693" "4148" "79083" "1101" "3158" "3169" "5346" "1408"
[177] "9547" "2922" "11283" "64499" "54829" "1524" "10234" "1580" "10647" "25893" "24141"
[188] "10351" "2330" "5304" "79846" "8614" "2625" "7021" "7802" "79689" "11122" "55351"
[199] "9" "4239" "5241" "10551" "10974" "79838" "79901" "57758" "4969"
geneはヒト遺伝子のENTREZのIDです。
この遺伝子リストを入力としてKEGGを対象としたGSEAを行ってみましょう。
library(clusterProfiler)
> kk <- enrichKEGG(gene = gene,
+ organism = 'hsa',
+ pvalueCutoff = 0.05)
> head(kk)
ID Description GeneRatio BgRatio pvalue p.adjust qvalue geneID Count
hsa04110 hsa04110 Cell cycle 11/90 123/7469 2.135411e-07 0.0000409999 4.023565e-05 8318/991/9133/890/983/4085/7272/1111/891/4174/9232 11
hsa04114 hsa04114 Oocyte meiosis 10/90 125/7469 2.216832e-06 0.0002128158 2.088489e-04 991/9133/983/4085/51806/6790/891/9232/3708/5241 10
hsa04218 hsa04218 Cellular senescence 10/90 160/7469 2.014286e-05 0.0012891427 1.265113e-03 2305/4605/9133/890/983/51806/1111/891/776/3708 10
hsa03320 hsa03320 PPAR signaling pathway 7/90 74/7469 2.724093e-05 0.0013075646 1.283191e-03 4312/9415/9370/5105/2167/3158/5346 7
hsa04914 hsa04914 Progesterone-mediated oocyte maturation 7/90 98/7469 1.657604e-04 0.0063651988 6.246549e-03 9133/890/983/4085/6790/891/5241 7
kk がGSEAの結果になります。
head(kk) でGSEAの結果を上から確からしいもの順にトップ5を出力しています。
もし遺伝子リストgeneが発現差異解析の結果であるなら、コントロールとターゲット間でCell cycle等に異変が起きているだろう、とこの結果から推測できるということになります。
次にkkの各列名の意味ですが
- IDはKEGG pathwayのID
- DescriptionはKEGG pathwayの名前
- GeneRatioは入力
geneの遺伝子の内、何個の遺伝子がそのKEGG pathwayにマップされるか - BgRatioは全KEGG pathwayにマップされる遺伝子の内、何個の遺伝子がそのKEGG pathwayにマップされるか
- p.adjust, qvalueは多重検定に伴う補正後のpvalueになり小さいほど確からしいことを示します。
入力geneに含まれる遺伝子がKEGG pathway上のどのgene productにmapされるかを可視化するには下記を実行します。
browseKEGG(kk, 'hsa04110')
終わりに
最も簡単なGSEAを体験していただく例を示しました。
追ってより進んだアンサンブル手法なども紹介できればと思っています。
- 前の記事 : 化合物データ処理ソフト RDKit
- 次の記事 : PythonとKerasでどうぶつしょうぎ
- 関連記事 :
