ex802 - nibb-gitc/gitc2025mar-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フォルダ（~/Desktop/gitc）に置いてあることを確認しよう。 ~/Desktop/gitc/data/IU/mouse の下にある GSE70520_edgeR.tab を用いる。Rの解析結果をプロジェクトごとに管理するため、前問同様に新たなディレクトリを作ってファイルをコピーしてから解析するとよいだろう（以下では、ディレクトリ名をex802とする）。
このデータは、上記のデータを edgeR を用いて解析し、topTagsで大きなnを指定して全遺伝子の結果を含むテーブルを作成し、それをwrite.tableでファイルに書き出したものである。これは、Rでread.tableでそのまま読み込むことができる。読み込んだら、headで中身を確認しよう。
> setwd("~/Desktop/gitc/data/IU/ex802")
> 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では、発現変動の有意性を評価するFDR値を topTagsのp.valueオプションで取り出せるが、それだけだとup regulateとdown regulateが区別できない。これは、logFCの符号で区別できる。ただし符号が＋のときが、up/downのどちらに対応するかは、元のカウントマトリクス上の並びに依存するので、元データに戻って確認する必要がある。
ここではedgeRの過程を省略したのでそのデータはないが、実際にedgeRで解析する際に読み込んでいるはずであり、仮にそのときのカウントマトリクスの変数をxとすると、上記１行目の遺伝子（logFCがマイナス）に対応する生カウントは以下のようにして調べられる。
# 以下はカウントマトリクスを 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あり、後半３つがなしに相当しているので、logFCがマイナスの時にはagonistありの条件で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)) # キーとして指定可能な遺伝子ID
[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") # keytypeをENSEMBLに変更して、対応するGO termを取得
'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:0098813 GO:0098813 nuclear chromosome segregation 45/497 261/11478
GO:0007059 GO:0007059 chromosome segregation 53/497 359/11478
GO:0000070 GO:0000070 mitotic sister chromatid segregation 31/497 169/11478
GO:0000280 GO:0000280 nuclear division 45/497 335/11478
GO:0000819 GO:0000819 sister chromatid segregation 33/497 202/11478
GO:0140014 GO:0140014 mitotic nuclear division 35/497 228/11478
結果を可視化するには、enrichplotをロードし、さまざまな可視化関数にenrichment解析の結果オブジェクトを与える形で実行する。
> library(enrichplot)
> barplot(cprof.up.bp.fisher, showCategory=10) # 棒グラフ
> goplot(cprof.up.bp.fisher, showCategory=10) # GO階層(DAG)表示
エンリッチしているGO termにそれぞれどういう遺伝子が含まれているかという情報も、結果オブジェクトに含まれており、可視化関数 cnetplot を使って可視化できる。ただし、記録されているのは入力に使った遺伝子リスト中のID（ここではENSEMBL ID)なので、わかりやすい名前(SYMBOL)で表示するために、setReadable関数を使う。また、upsetplotはどのGO term群に含まれるかによって遺伝子を分類して頻度をグラフ表示したものである。これらの表示により、得られた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)
ID Description setSize
GO:1905818 GO:1905818 regulation of chromosome separation 66
GO:0000070 GO:0000070 mitotic sister chromatid segregation 169
GO:0098813 GO:0098813 nuclear chromosome segregation 261
GO:0000819 GO:0000819 sister chromatid segregation 202
GO:0051304 GO:0051304 chromosome separation 75
GO:0007059 GO:0007059 chromosome segregation 359
可視化関数として、ridgeplotとgseplotがある。ridgeplotはshowCategoryで与えた数の上位GO termについて、スコアの分布と一様分布との差分を表示する。gseplotは、geneSetIDで指定された順位のGO termについて、スコアの分布と一様分布との差分を表示する。
> ridgeplot(cprof.bp.gsea, showCategory=30)
> gseaplot(cprof.bp.gsea, geneSetID=1)
（追加課題）オントロジーの種類がMF, CCの場合についても解析せよ。
Step 6: (Optional) TopGOによるエンリッチメント解析の実行
TopGOを使う際も、必要な準備はclusterProfilerとほぼ同様である。Step 0, 2, 4 はこのまま実行すること。TopGOでは、Fisher' exact 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:0008608 attachment of spindle microtubules to ki... 49 22 3.94 17 5.7e-11 4.0e-12 1.00
2 GO:0000819 sister chromatid segregation 202 60 16.24 3 1.3e-10 8.0e-20 1.00
3 GO:0007094 mitotic spindle assembly checkpoint sign... 44 19 3.54 28 2.7e-10 2.7e-10 0.98
4 GO:0051301 cell division 519 87 41.74 23 5.7e-08 2.0e-11 0.95
5 GO:0051984 positive regulation of chromosome segreg... 29 13 2.33 62 1.1e-07 1.1e-07 0.99
可視化のための関数としてはshowSigOfNodesがあり、有意性ありとして抽出されたGO termを、GO グラフ階層上に表示することができる。この表示から、algorithm="classic" の出力が冗長であり (同じGO termの周辺に集中している)、algorithm="weight"を用いることで冗長性が除かれ、より広範なGO termが抽出されていることが確認できる。
> showSigOfNodes(topgoData, score(resultFisher), firstSigNodes=6, useInfo='pval')
> showSigOfNodes(topgoData, score(resultWeightFisher), firstSigNodes=6, useInfo='pval')