大阪大学医学部 Python会 (情報医科学研究会)

Now is better than never.

Gene Set Enrichment Analysis (GSEA)入門

2018-11-27(Tue) - Posted by 西田 in 技術ブログ    tag:Bioinformatics

Contents

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

    KEGG map

    終わりに

    最も簡単なGSEAを体験していただく例を示しました。

    追ってより進んだアンサンブル手法なども紹介できればと思っています。