ex802 - nibb-gitc/gitc2021mar-rnaseq GitHub Wiki
ex802: GO enrichment analysis using public annotation databases
演習802: 公開アノテーションデータを用いたGOエンリッチメント解析
ここでは、公開アノテーションデータを使ってGOエンリッチメント解析を行う例として、マウスのRNA-seq解析データ(GEO accession: GSE70520) を用いた解析を行ってみよう。このデータは、炎症性疾患や感染症に伴う貧血や血小板減少などの症状に、Toll like receptor(TLR)を介したシグナル伝達による食細胞の活性化が関与していることを示した研究の過程でとられており、TLR7 agonistありの条件となしの条件で発達させたマクロファージにおける遺伝子発現を比較したものである。
Step 0: データ
Rを使うので、ローカル環境で作業する。~/gitc/data/IU/mouse の下にある GSE70520_edgeR.tab を用いる。これは、上記のデータを edgeR を用いて解析し、topTagsで大きなnを指定して全遺伝子の結果を含むテーブルを作成し、それをwrite.tableでファイルに書き出したものである。これは、Rでread.tableでそのまま読み込むことができる。読み込んだら、headで中身を確認しよう。
> setwd("~/gitc/data/IU/mouse")
> etab <- read.table("GSE70520_edgeR.tab")
> head(etab)
logFC logCPM F PValue FDR
ENSMUSG00000004105 -6.908256 5.582768 1086.1199 4.181224e-09 3.059046e-05
ENSMUSG00000039616 -2.418207 5.760821 852.4535 9.923661e-09 3.059046e-05
ENSMUSG00000006800 -4.164396 6.332824 825.1360 1.114531e-08 3.059046e-05
Step 1: Fisher's exact testのための遺伝子リスト作成
興味のある遺伝子セットを抽出する。edgeRでは、topTagsのp.valueオプションで取り出せるが、それだけだとup regulateとdown regulateが区別できない。これは、logFCの符号で区別できる。ただし符号が＋のときが、up/downのどちらに対応するかは、カウントマトリクス上の並びに依存する。ここではその過程は省略しているが、edgeRで解析する際にそのカウントマトリクスを読み込んだはずであり、仮にそのときの変数をxとすると、以下のようにして調べられる。
> x["ENSMUSG00000004105",]
lib3731.15755752 lib3732.15757758 lib3733.15756756
ENSMUSG00000004105 26 23 16
lib3734.15756757 lib3735.15757759 lib3736.15755753
ENSMUSG00000004105 2774 2865 2176
前半３つがTLR7 agonistあり、後半３つがなしに相当しているので、agonistありを中心に見ると、logFCがマイナスの時にdown reguateされていることがわかる。そこで、logFC > 0 がup regulateの条件となるが、ここでは少し条件を厳しくしてlogFC > 1 としてみよう。これはfold change 2倍以上という条件に相当する。これに FDR < 0.01という条件もつけて、以下のようにして抽出する。
> upreg <- subset(etab, logFC > 1.0 & FDR < 0.01)
> head(upreg)
logFC logCPM F PValue FDR
ENSMUSG00000026399 7.764379 5.549768 569.3602 4.173203e-08 3.706083e-05
ENSMUSG00000020717 4.702431 6.315477 536.1777 5.164914e-08 3.742434e-05
ENSMUSG00000075254 5.553891 6.151117 532.1209 5.306030e-08 3.742434e-05
実際に必要なのは、遺伝子名のリストであり、以下のように作成する。全体のリストも同様にして作成する。
> upregGenes <- rownames(upreg)
> allGenes <- rownames(etab)
Step 2: アノテーションデータベースの準備
マウスのアノテーションデータとしては、org.Mm.eg.db パッケージを利用できる。
> library(org.Mm.eg.db)
ロードすると、パッケージ名と同じorg.Mm.eg.dbという名前のオブジェクトを介してデータベースにアクセスできる。これは、AnnotationDbiパッケージで定義された共通インターフェイスを介して利用できるデータベースの一つであり、key（ここでは遺伝子ID）とcolumn（ここでは各遺伝子の属性）の２つを指定して、select関数を使って情報を取り出す。keyとしては、デフォルトではEntrezIDがとられるが、keytypeを指定することで他の属性に変更できる。以下はデータベース検索例（後のステップには必要ない）。
> head(keys(org.Mm.eg.db))
[1] "11287" "11298" "11302" "11303" "11304" "11305"
> columns(org.Mm.eg.db)
[1] "ACCNUM" "ALIAS" "ENSEMBL" "ENSEMBLPROT" "ENSEMBLTRANS"
[6] "ENTREZID" "ENZYME" "EVIDENCE" "EVIDENCEALL" "GENENAME"
[11] "GO" "GOALL" "IPI" "MGI" "ONTOLOGY"
[16] "ONTOLOGYALL" "PATH" "PFAM" "PMID" "PROSITE"
[21] "REFSEQ" "SYMBOL" "UNIGENE" "UNIPROT"
> select(org.Mm.eg.db, keys=c("11287","11298"), columns=c("ENSEMBL","SYMBOL"))
'select()' returned 1:1 mapping between keys and columns
ENTREZID ENSEMBL SYMBOL
1 11287 ENSMUSG00000030359 Pzp
2 11298 ENSMUSG00000020804 Aanat
....
> select(org.Mm.eg.db, keytype="ENSEMBL", keys="ENSMUSG00000026399", columns="GO")
'select()' returned 1:many mapping between keys and columns
ENSEMBL GO EVIDENCE ONTOLOGY
1 ENSMUSG00000026399 GO:0001618 ISO MF
2 ENSMUSG00000026399 GO:0002376 IEA BP
3 ENSMUSG00000026399 GO:0004857 ISO MF
Step 3: clusterProfilerによるFisher's exact test の実行
clusterProfilerをロードしenrichGO関数を実行する。引数にはターゲットの遺伝子リスト、全体の遺伝子リスト、アノテーションデータベース名、オントロジーの種類(BP/MF/CC)、データベースにアクセスする際のkeytype、(adjusted) p-valueの閾値をそれぞれ指定する。結果はエンリッチしているGO termがqvalueの順に並んでおり、headで上位のGO termを取り出すことができる。
> library(clusterProfiler)
> cprof.up.bp.fisher <- enrichGO(gene=upregGenes, universe=allGenes, OrgDb=org.Mm.eg.db, ont="BP", keyType="ENSEMBL", pvalueCutoff=0.01)
> head(cprof.up.bp.fisher)
ID Description GeneRatio BgRatio
GO:0007059 GO:0007059 chromosome segregation 36/490 270/11352
GO:0098813 GO:0098813 nuclear chromosome segregation 30/490 211/11352
GO:0046651 GO:0046651 lymphocyte proliferation 29/490 220/11352
GO:0032943 GO:0032943 mononuclear cell proliferation 29/490 222/11352
GO:0070661 GO:0070661 leukocyte proliferation 29/490 235/11352
GO:0000819 GO:0000819 sister chromatid segregation 23/490 160/11352
結果を可視化するには、enrichplotをロードし、さまざまな可視化関数にenrichment解析の結果オブジェクトを与える形で実行する。
> library(enrichplot)
> barplot(cprof.up.bp.fisher, showCategory=10)
> goplot(cprof.up.bp.fisher, showCategory=10)
エンリッチしているGO termにそれぞれどういう遺伝子が含まれているかという情報も、結果オブジェクトに含まれており、可視化関数 cnetplot を使って可視化できる。ただし、記録されているのは入力に使った遺伝子リスト中のID（ここではENSEMBL ID)なので、わかりやすい名前(SYMBOL)で表示するために、setReadable関数を使う。また、upsetplotはどのGO term群に含まれるかによって遺伝子を分類して頻度をグラフ表示したものである。
> cprof.up.bp.fisher <- setReadable(cprof.up.bp.fisher, 'org.Mm.eg.db', 'ENSEMBL')
> cnetplot(cprof.up.bp.fisher)
> upsetplot(cprof.up.bp.fisher)
cnetplotは描画に時間がかかる。途中で止めたい場合は左上のStopボタンを押すなどして止める。
（追加課題）オントロジーの種類がMF, CCの場合についても解析せよ。また、logFC < -1.0 & FDR < 0.01 の条件でdown regulateされた遺伝子を抽出し、同様に解析せよ。
Step 4: Gene Set Enrichment Analysisのための遺伝子ごとのスコアリスト作成
GSEAを実行するためには、発現解析結果を反映したなんらかのスコアを定義し、その順に遺伝子を並べ替えた遺伝子スコアリスト（ベクトル）を作成する必要がある。ここでは発現スコア(expScore)の絶対値はFDRの値とし、up/downを区別するためにlogFCの符号をかけたものと定義する。
> expScore <- sign(etab$logFC) * -log(etab$FDR)
> head(expScore)
[1] -10.39482 -10.39482 -10.39482 -10.39482 -10.39482 -10.39482
このスコアベクトルに遺伝子名を対応づける必要がある。これは、names関数を用いて各値に名前をつけることで行う。その後、スコアの値順でソートする。
> names(expScore) <- rownames(etab)
> expScore <- sort(expScore, decreasing=TRUE)
> head(expScore)
ENSMUSG00000026399 ENSMUSG00000020717 ENSMUSG00000075254 ENSMUSG00000004359 ENSMUSG00000038807 ENSMUSG00000025964
10.202950 10.193189 10.193189 10.193189 10.193189 9.948559
スコアが降順でソートされていることはclusterProfilerの gseGO 関数の仕様として決められているが、GSEA解析では一様分布からのずれを見るので、いずれにしてもup regulateとdown regulateは同時に検出される。
Step 5: clusterProfilerによるGene Set Enrichment Analysis の実行
gseGO関数で行う。引数には、geneListとして先ほど作成したスコアベクトル(expScore)を与える。このほかの引数は、概ねenrichGOと同様である。
> cprof.bp.gsea <- gseGO(geneList=expScore, OrgDb=org.Mm.eg.db, ont="BP", keyType="ENSEMBL", pvalueCutoff=0.01)
> head(cprof.bp.gsea)
可視化関数として、ridgeplotとgseplotがある。ridgeplotはshowCategoryで与えた数の上位GO termについて、スコアの分布と一様分布との差分を表示する。gseplotは、指定された順位のGO termについて、スコアの分布と一様分布との差分を表示する。
> ridgeplot(cprof.bp.gsea, showCategory=30)
> gseaplot(cprof.bp.gsea, geneSetID=1)
（追加課題）オントロジーの種類がMF, CCの場合についても解析せよ。
Optional: TopGOによるエンリッチメント解析の実行
TopGOを使う際も、必要な準備はclusterProfilerとほぼ同様である。Step 0, 2, 4 はこのまま実行すること。TopGOでは、Fisher testにかける遺伝子リストをStep 1のように明示的に作る代わりに、Step 4で作成した遺伝子ごとのスコアベクトルから遺伝子を抽出する関数を与えることで、Fisher's exact testとKolmogorov-Smirnov test (=GSEA)などを同時に実行できるという特徴がある。スコアベクトルはStep 4で定義した expScore （定義はsign(logFC) * -log(FDR)）とし、FDR<0.01を抽出条件とすると、抽出関数は以下のようになる。
> selfun <- function(x) {return (x > -log(0.01))}
これを用いて、TopGOを実行するために必要な情報を集めた"topGOdataオブジェクト"を作成する。
> library(topGO)
> topgoData <- new("topGOdata", ontology="BP", allGenes=expScore, geneSel=selfun, annot=annFUN.org, mapping="org.Mm.eg.db", ID="ENSEMBL")
これを用いて、statistic(fisher/ksなど)とalgorithm(classic/elim/weight01/weightなど)を指定して検定を実行する。
> resultFisher <- runTest(topgoData, algorithm="classic", statistic="fisher")
> resultWeightFisher <- runTest(topgoData, algorithm="weight", statistic="fisher")
> resultElimKS <- runTest(topgoData, algorithm="elim", statistic="ks")
複数の方法で行なった結果は、GenTable関数を使って一つにまとめて表示できる。2番目以降の引数 name=result の指定に従って、結果オブジェクトresultにおける各GO termのp-valueをnameという名前のカラム上に順に表示する。デフォルトでは最初に指定したカラムのp-value順にソートされるが、orderBy 引数で変更できる。
> GenTable(topgoData, weightFisher=resultWeightFisher, classicFisher=resultFisher, elimKS=resultElimKS)
GO.ID Term Annotated Significant Expected Rank in classicFisher weightFisher classicFisher elimKS
1 GO:0098813 nuclear chromosome segregation 211 55 16.93 2 1.5e-08 1.5e-15 0.979
2 GO:0051301 cell division 486 84 39.00 12 1.3e-07 7.8e-12 0.952
3 GO:0051383 kinetochore organization 21 10 1.69 41 1.6e-06 1.6e-06 0.989
4 GO:1902850 microtubule cytoskeleton organization in... 123 28 9.87 29 1.4e-05 3.3e-07 0.979
5 GO:0030261 chromosome condensation 31 11 2.49 79 1.6e-05 1.6e-05 0.950
可視化のための関数としてはshowSigOfNodesがあり、有意性ありとして抽出されたGO termを、GO グラフ階層上に表示することができる。この表示から、algorithm="classic" の出力が冗長であり、algorithm="weight"を用いることで冗長性が除かれ、より多様なGO termが抽出されていることが確認できる。
> showSigOfNodes(topgoData, score(resultFisher), firstSigNodes=6, useInfo='pval')
> showSigOfNodes(topgoData, score(resultWeightFisher), firstSigNodes=6, useInfo='pval')